1 %% Introduction to signal analysis in Matlab: GLM I
   2 % CBU workshop, Ed Sohoglu 2018
   3 
   4 %% Regression
   5 clearvars
   6 close all
   7 rng(1) % for reproducability- code below will generate 'random' numbers in a predictable way
   8 
   9 nReg = 1; % number of regressors
  10 nSubj = 10; % number of subjects
  11 
  12 % make design matrix
  13 age = randi([20 60],nSubj,1);
  14 X = age; % X = design matrix
  15 X = [X ones(size(X,1),1)]; % add constant
  16 
  17 % make data using GLM equation Y = B1*X1 + B2*X2 + B3*X3 etc. + Constant + Error 
  18 % here beta 1 represents influence of age covariate and beta 2 is the constant
  19 Y = 3*X(:,1) + 500 + 20*randn(nReg*nSubj,1);
  20 
  21 % plot data and design matrix
  22 figure;
  23 subplot(1,2,1); scatter(age,Y);
  24 subplot(1,2,2); imagesc(X); colormap(gray); colorbar;
  25 
  26 % compute betas
  27 betas = pinv(X)*Y;
  28 
  29 % error sum of squares
  30 prediction = X*betas;
  31 res = Y-prediction;
  32 SS_error = sum(res.^2);
  33 subplot(1,2,1); hold on; plot(age,prediction); hold off % plot data and prediction
  34 
  35 % effect sum of squares
  36 prediction = X(:,1:nReg)*betas(1:nReg); % Only use regressor/group regressors (ignore covariates and constant terms) 
  37 SS_effect = sum((prediction-mean(prediction)).^2);
  38 
  39 % compute mean squares
  40 df_effect = nReg;
  41 df_error = size(Y,1)-df_effect-1;
  42 MS_effect = SS_effect./df_effect;
  43 MS_error = SS_error/df_error;
  44 
  45 % compute F statistic, p-value and R2
  46 F = MS_effect./MS_error;
  47 p = fcdf(F,df_effect,df_error,'upper'); % requires statistics tooolbox
  48 
  49 SS_total = sum((Y-mean(Y)).^2);
  50 R2 = SS_effect/SS_total; % R2 is the proportion of variance explained (take square root to compute correlation coefficient)
  51 
  52 %% One-way ANOVA
  53 clearvars
  54 close all
  55 rng(1) % for reproducability- code below will generate 'random' numbers in a predictable way
  56 
  57 nReg = 3; % number of regressors/groups
  58 nSubj = 10; % number of subjects per group
  59 
  60 % make design matrix
  61 X = kron(eye(nReg),ones(nSubj,1)); % X = design matrix
  62 X = [X ones(size(X,1),1)]; % add constant
  63 
  64 % make data using GLM equation Y = B1*X1 + B2*X2 + B3*X3 etc. + Constant + Error  
  65 % here betas 1 to 3 represent means for the three groups and beta 4 is the constant
  66 Y = 300*X(:,1) + 200*X(:,2) + 400*X(:,3) + 500 + 20*randn(nReg*nSubj,1);
  67 
  68 % plot data and design matrix
  69 figure;
  70 subplot(1,2,1); scatter(1:nSubj*nReg,Y);
  71 subplot(1,2,2); imagesc(X); colormap(gray); colorbar;
  72 
  73 % compute betas
  74 betas = pinv(X)*Y;
  75 
  76 % error sum of squares
  77 prediction = X*betas;
  78 res = Y-prediction;
  79 SS_error = sum(res.^2);
  80 subplot(1,2,1); hold on; plot(1:nSubj*nReg,prediction); hold off % plot data and prediction
  81 
  82 % effect sum of squares
  83 prediction = X(:,1:nReg)*betas(1:nReg); % ignore constant terms 
  84 SS_effect = sum((prediction-mean(prediction)).^2);
  85 
  86 % compute mean squares
  87 df_effect = nReg-1;
  88 df_error = size(Y,1)-df_effect-1;
  89 MS_effect = SS_effect./df_effect;
  90 MS_error = SS_error/df_error;
  91 
  92 % compute F statistic, p-value and R2
  93 F = MS_effect./MS_error;
  94 p = fcdf(F,df_effect,df_error,'upper'); % requires statistics tooolbox
  95 
  96 SS_total = sum((Y-mean(Y)).^2);
  97 R2 = SS_effect/SS_total; % R2 is the proportion of variance explained (take square root to compute correlation coefficient)
  98 
  99 %% One-way repeated measures ANOVA
 100 clearvars
 101 close all
 102 rng(1) % for reproducability- code below will generate 'random' numbers in a predictable way
 103 
 104 nReg = 3; % number of regressors/conditions
 105 nSubj = 10; % number of subjects per condition
 106 
 107 % make design matrix
 108 X = kron(eye(nReg),ones(nSubj,1)); % X = design matrix
 109 X = [X kron(ones(nReg,1),eye(nSubj))]; % add constant for each subject
 110 
 111 % make data using GLM equation Y = B1*X1 + B2*X2 + B3*X3 etc. + Constant + Error   
 112 % here betas 1 to 3 represent means for the three groups, beta 4 is the age
 113 % covariate and beta 5 is the constant
 114 subject_means = 200*randn(nSubj,1)+500;
 115 subject_means = kron(ones(nReg,1),subject_means);
 116 Y = 300*X(:,1) + 200*X(:,2) + 400*X(:,3) + subject_means + 20*randn(nReg*nSubj,1);
 117 
 118 % plot data and design matrix
 119 Y2plot = reshape(Y,nSubj,nReg);
 120 figure; subplot(1,2,1); hold on
 121 for i=1:nReg
 122     scatter(1:nSubj,Y2plot(:,i));
 123 end
 124 hold off
 125 subplot(1,2,2); imagesc(X); colormap(gray); colorbar;
 126 
 127 % compute betas
 128 betas = pinv(X)*Y;
 129 
 130 % error sum of squares
 131 prediction = X*betas;
 132 res = Y-prediction;
 133 SS_error = sum(res.^2);
 134 prediction2plot = reshape(prediction,nSubj,nReg);
 135 subplot(1,2,1); hold on; ax = gca; ax.ColorOrderIndex = 1; plot(1:nSubj,prediction2plot); hold off % plot data and prediction
 136 
 137 % effect sum of squares
 138 prediction = X(:,1:nReg)*betas(1:nReg); % ignore constant terms
 139 SS_effect = sum((prediction-mean(prediction)).^2);
 140 
 141 % compute mean squares
 142 df_effect = nReg-1;
 143 df_error = size(Y,1)-df_effect-nSubj;
 144 MS_effect = SS_effect./df_effect;
 145 MS_error = SS_error/df_error;
 146 
 147 % compute F statistic, p-value and R2
 148 F = MS_effect./MS_error;
 149 p = fcdf(F,df_effect,df_error,'upper'); % requires statistics tooolbox
 150 
 151 SS_total = sum((Y-mean(Y)).^2);
 152 R2 = SS_effect/SS_total; % R2 is the proportion of variance explained (take square root to compute correlation coefficient)
 153 
 154 %% ANCOVA
 155 
 156 % In this exercise, you will combine what you learned from the Regression
 157 % and One-way ANOVA examples to build an ANCOVA model.
 158 % The data are RTs from three groups of subjects (N = 10 per group). You
 159 % will use ANCOVA to test for the group effect, after controlling for the
 160 % influence of age.
 161 
 162 nReg = 3;
 163 nSubj = 10; % per group
 164 rng(1) % for reproducability- code below will generate 'random' numbers in a predictable way
 165 
 166 % make design matrix
 167 X = kron(eye(nReg),ones(nSubj,1)); % X = design matrix
 168 age = linspace(20,70,nReg*nSubj)';
 169 X = [X age]; % Add age covariate
 170 X = [X ones(size(X,1),1)]; % add constant
 171 figure; imagesc(X); colorbar;
 172 
 173 % make data using GLM equation Y = B1*X1 + B2*X2 + B3*X3 etc. + Constant + Error   
 174 % here betas 1 to 3 represent means for the three groups, beta 4 is the age
 175 % covariate and beta 5 is the constant
 176 Y = 300*X(:,1) + 200*X(:,2) + 400*X(:,3) + 3*X(:,4) + 500 + 20*randn(nReg*nSubj,1);
 177 
 178 % plot data and design matrix
 179 figure;
 180 subplot(1,2,1); scatter(1:nSubj*nReg,Y);
 181 subplot(1,2,2); imagesc(X); colormap(gray); colorbar;
 182 
 183 % compute betas
 184 betas = pinv(X)*Y;
 185 
 186 % error sum of squares
 187 prediction = X*betas;
 188 res = Y-prediction;
 189 SS_error = sum(res.^2);
 190 subplot(1,2,1); hold on; plot(1:nSubj*nReg,prediction); hold off % plot data and prediction
 191 
 192 % effect sum of squares
 193 prediction = X(:,1:nReg)*betas(1:nReg); % ignore constant terms (and covariate)
 194 SS_effect = sum((prediction-mean(prediction)).^2);
 195 
 196 % compute mean squares
 197 df_effect = nReg-1;
 198 df_error = size(Y,1)-df_effect-nSubj;
 199 MS_effect = SS_effect./df_effect;
 200 MS_error = SS_error/df_error;
 201 
 202 % compute F statistic, p-value and R2
 203 F = MS_effect./MS_error;
 204 p = fcdf(F,df_effect,df_error,'upper'); % requires statistics tooolbox
 205 
 206 %% One-way repeated measures ANOVA
 207 clearvars
 208 close all
 209 rng(1) % for reproducability- code below will generate 'random' numbers in a predictable way
 210 
 211 nReg = 4; % number of regressors/conditions
 212 nSubj = 10; % number of subjects per condition
 213 
 214 % make design matrix
 215 A = [repmat(1,nSubj,1); repmat(-1,nSubj,1);];
 216 A = kron(A,[1 1]');
 217 B = [repmat(1,nSubj,1); repmat(-1,nSubj,1);];
 218 B = kron([1 1]',B);
 219 AB = A.*B;
 220 X = [A B AB ones(size(A,1),1)];
 221 X = [X kron([[1 1 -1 -1]' [1 -1 1 -1]' [1 -1 -1 1]' [1 1 1 1]' ],eye(nSubj))]; % add constant for each subject
 222 
 223 % % make data using GLM equation Y = B1*X1 + B2*X2 + B3*X3 etc. + Constant + Error   
 224 % % here betas 1 to 3 represent means for the three groups, beta 4 is the age
 225 % % covariate and beta 5 is the constant
 226 % subject_means = 200*randn(nSubj,1)+500;
 227 % subject_means = kron(ones(nReg,1),subject_means);
 228 % Y = 300*X(:,1) + 200*X(:,2) + 400*X(:,3) + subject_means + 100*randn(nReg*nSubj,1);
 229 Y = [1+randn(10,1); 2+randn(10,1); 3+randn(10,1); 4+randn(10,1)];
 230 
 231 % plot data and design matrix
 232 Y2plot = reshape(Y,nSubj,nReg);
 233 figure; subplot(1,2,1); hold on
 234 for i=1:nReg
 235     scatter(1:nSubj,Y2plot(:,i));
 236 end
 237 hold off
 238 subplot(1,2,2); imagesc(X); colormap(gray); colorbar;
 239 
 240 % compute betas
 241 betas = pinv(X)*Y;
 242 
 243 % error sum of squares
 244 prediction = X*betas;
 245 res = Y-prediction;
 246 SS_error = sum(res.^2);
 247 prediction2plot = reshape(prediction,nSubj,nReg);
 248 subplot(1,2,1); hold on; ax = gca; ax.ColorOrderIndex = 1; plot(1:nSubj,prediction2plot); hold off % plot data and prediction
 249 
 250 % effect sum of squares
 251 prediction = X(:,[1 5 6:end])*betas([1 5 6:end]); % ignore constant terms
 252 SS_effect = sum((prediction-mean(prediction)).^2);
 253 
 254 % compute mean squares
 255 df_effect = 1;
 256 df_error = 9;
 257 MS_effect = SS_effect./df_effect;
 258 MS_error = SS_error/df_error;
 259 
 260 % compute F statistic, p-value and R2
 261 F = MS_effect./MS_error;
 262 p = fcdf(F,df_effect,df_error,'upper'); % requires statistics tooolbox
 263 
 264 SS_total = sum((Y-mean(Y)).^2);
 265 R2 = SS_effect/SS_total; % R2 is the proportion of variance explained (take square root to compute correlation coefficient)