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(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. + Error  
  65 % here betas 1 to 3 represent means for the three groups and beta 4 is the constant
  66 Y = 500*X(:,1) + 700*X(:,2) + 600*X(:,3) + 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. + Error   
 112 % here betas 1 to 3 represent means for the three conditions, betas 4 to 13
 113 % are the subject-specific constants
 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 % You will need to generate the RT data, arranged as three groups of subjects (N = 10 per group, similar to one-way ANOVA example above.
 159 % You will also need to generate the age covariate (similar to Regression example above)
 160 % Use ANCOVA to test for the group effect, after controlling for the influence of age.
 161 
 162 nReg = 3; % number of (group) regressors
 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 
 172 % make data using GLM equation Y = B1*X1 + B2*X2 + B3*X3 etc. + Error   
 173 % here betas 1 to 3 represent means for the three groups, beta 4 is the age
 174 % covariate and beta 5 is the constant
 175 Y = 500*X(:,1) + 700*X(:,2) + 600*X(:,3) + 3*X(:,4) + 20*randn(nReg*nSubj,1);
 176 
 177 % plot data and design matrix
 178 figure;
 179 subplot(1,2,1); scatter(1:nSubj*nReg,Y);
 180 subplot(1,2,2); imagesc(X); colormap(gray); colorbar;
 181 
 182 % compute betas
 183 betas = pinv(X)*Y;
 184 
 185 % error sum of squares
 186 prediction = X*betas;
 187 res = Y-prediction;
 188 SS_error = sum(res.^2);
 189 subplot(1,2,1); hold on; plot(1:nSubj*nReg,prediction); hold off % plot data and prediction
 190 
 191 % effect sum of squares
 192 prediction = X(:,1:nReg)*betas(1:nReg); % ignore constant terms (and covariate)
 193 SS_effect = sum((prediction-mean(prediction)).^2);
 194 
 195 % compute mean squares
 196 df_effect = nReg-1;
 197 df_error = size(Y,1)-df_effect-nSubj;
 198 MS_effect = SS_effect./df_effect;
 199 MS_error = SS_error/df_error;
 200 
 201 % compute F statistic, p-value and R2
 202 F = MS_effect./MS_error;
 203 p = fcdf(F,df_effect,df_error,'upper'); % requires statistics tooolbox
 204 
 205 %% Two-way (2x2) ANOVA
 206 
 207 % In this exercise, you can try a more complex two-way (2x2) anova.
 208 % Try to assess each of the three effects in this design
 209 % (main effect A, main effect B and AxB interaction)
 210 % Tip: each factor (column) in the design matrix will need to code for 
 211 % differences between the levels of that factor (e.g. for factor A: kron([1 1 -1 -1]',ones(nSubj,1))').
 212 % The column for the interaction can be formed by element-wise
 213 % multiplication of factors A and B.
 214 
 215 clearvars
 216 close all
 217 rng(1) % for reproducability- code below will generate 'random' numbers in a predictable way
 218 
 219 nReg = 4; % number of regressors (conditions)
 220 nSubj = 10; % number of subjects per group
 221 chooseFactor = 1; % which factor do you want to assess? (1 for A, 2 for B, 3 for interaction)
 222 
 223 % make design matrix
 224 A = kron([1 1 -1 -1]',ones(nSubj,1)); % Factor A
 225 B = kron([1 -1 1 -1]',ones(nSubj,1)); % Factor B
 226 AB = A.*B; % AxB interaction
 227 X = [A B AB];
 228 X = [X ones(size(X,1),1)]; % add constant
 229 
 230 % make data using GLM equation Y = B1*X1 + B2*X2 + B3*X3 etc. + Error   
 231 % here betas 1 to 4 represent means for the four groups (2x2 design) and beta 5 is the constant
 232 Y = 500*X(:,1) + 700*X(:,2) + 600*X(:,3) + 20*randn(nReg*nSubj,1);
 233 
 234 % plot data and design matrix
 235 Y2plot = reshape(Y,nSubj,nReg);
 236 figure; subplot(1,2,1); hold on
 237 for i=1:nReg
 238     scatter(1:nSubj,Y2plot(:,i));
 239 end
 240 hold off
 241 subplot(1,2,2); imagesc(X); colormap(gray); colorbar;
 242 
 243 % compute betas
 244 betas = pinv(X)*Y;
 245 
 246 % error sum of squares
 247 prediction = X*betas;
 248 res = Y-prediction;
 249 SS_error = sum(res.^2);
 250 prediction2plot = reshape(prediction,nSubj,nReg);
 251 subplot(1,2,1); hold on; ax = gca; ax.ColorOrderIndex = 1; plot(1:nSubj,prediction2plot); hold off % plot data and prediction
 252 
 253 % effect sum of squares
 254 prediction = X(:,chooseFactor)*betas(chooseFactor); % only include factor you want to assess
 255 SS_effect = sum((prediction-mean(prediction)).^2);
 256 
 257 % compute mean squares
 258 df_effect = 1; % Degrees of freedom for the factor under consideration
 259 df_error = size(Y,1)-(nReg-1)-1; % But for error degrees of freedom, use full model
 260 MS_effect = SS_effect./df_effect;
 261 MS_error = SS_error/df_error;
 262 
 263 % compute F statistic, p-value and R2
 264 F = MS_effect./MS_error;
 265 p = fcdf(F,df_effect,df_error,'upper'); % requires statistics tooolbox
 266 
 267 SS_total = sum((Y-mean(Y)).^2);
 268 R2 = SS_effect/SS_total; % R2 is the proportion of variance explained (take square root to compute correlation coefficient)