attachment:GLM_workshop.m of SignalAnalysisMatlabSchedule - Methods
location: attachment:GLM_workshop.m of SignalAnalysisMatlabSchedule

Attachment 'GLM_workshop.m'

Download

   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)

Attached Files

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.
  • [get | view] (2018-11-15 14:50:04, 2556.1 KB) [[attachment:ArcherBoyd+Guerit_vocoder_Nov18.zip]]
  • [get | view] (2019-11-21 17:21:33, 1661.8 KB) [[attachment:CBU_FunctionsVocoder_22Nov19.zip]]
  • [get | view] (2019-11-19 16:53:57, 4.6 KB) [[attachment:Code_Optimisation_20Nov19.zip]]
  • [get | view] (2016-03-09 10:59:16, 259.0 KB) [[attachment:DP_Optimisation_2016.zip]]
  • [get | view] (2015-01-20 17:22:38, 614.9 KB) [[attachment:EMEGdata.mat]]
  • [get | view] (2016-02-03 12:32:14, 26.8 KB) [[attachment:ExampleScripts.zip]]
  • [get | view] (2018-01-12 13:39:53, 3150.4 KB) [[attachment:Filtering&Oscillations_ForUpload_TallieJan18.pdf]]
  • [get | view] (2018-07-19 14:49:35, 2.5 KB) [[attachment:Functions and Calculus.m]]
  • [get | view] (2018-07-19 14:49:30, 433.8 KB) [[attachment:Functions and Calculus.pdf]]
  • [get | view] (2015-02-02 09:27:07, 522.0 KB) [[attachment:GLM1.pdf]]
  • [get | view] (2019-11-19 13:06:58, 790.4 KB) [[attachment:GLM1_19Nov19.pdf]]
  • [get | view] (2019-11-19 13:07:11, 2.5 KB) [[attachment:GLM1_19Nov19.zip]]
  • [get | view] (2018-01-12 13:00:28, 9.0 KB) [[attachment:GLM1_Jan2018.m]]
  • [get | view] (2018-01-12 13:02:14, 354.1 KB) [[attachment:GLM1_Jan2018.pdf]]
  • [get | view] (2018-11-13 11:34:38, 11.5 KB) [[attachment:GLM1_Nov2018.m]]
  • [get | view] (2018-11-13 11:34:31, 1417.1 KB) [[attachment:GLM1_Nov2018.pptx]]
  • [get | view] (2018-11-13 12:07:58, 11.5 KB) [[attachment:GLM1_Nov2018_ES18.m]]
  • [get | view] (2018-11-13 12:07:52, 271.1 KB) [[attachment:GLM1_Nov2018_ES18.pdf]]
  • [get | view] (2015-02-02 09:27:27, 2.1 KB) [[attachment:GLM1_script.m]]
  • [get | view] (2016-02-11 18:26:35, 560.4 KB) [[attachment:GLM_part1.pdf]]
  • [get | view] (2018-01-08 09:49:09, 8.5 KB) [[attachment:GLM_workshop.m]]
  • [get | view] (2015-02-18 12:51:09, 580.2 KB) [[attachment:IntroCalculus 18Feb15.pdf]]
  • [get | view] (2015-03-02 18:21:14, 4.3 KB) [[attachment:IntroFilteringOscillations_25Feb15.zip]]
  • [get | view] (2016-02-18 15:52:26, 237.9 KB) [[attachment:IntroGLM2_17Feb16.pdf]]
  • [get | view] (2015-02-04 11:58:10, 125.6 KB) [[attachment:IntroGLM2_4Feb15.pdf]]
  • [get | view] (2016-02-04 17:19:34, 445.3 KB) [[attachment:IntroMatrixAlgebra.pdf]]
  • [get | view] (2015-03-04 12:00:10, 1.3 KB) [[attachment:IntroOptimisation_4Mar15.zip]]
  • [get | view] (2018-01-09 12:44:37, 19.3 KB) [[attachment:Matlab_Files.zip]]
  • [get | view] (2015-01-20 17:29:01, 6.1 KB) [[attachment:MatrixAlgebra.m]]
  • [get | view] (2018-11-13 12:08:14, 917.4 KB) [[attachment:Matrix_Algebra_AT18.pdf]]
  • [get | view] (2015-03-02 18:19:54, 260.9 KB) [[attachment:Olaf FilteringOscillations 25Feb15.pdf]]
  • [get | view] (2015-02-18 12:51:00, 1.8 KB) [[attachment:Olaf IntroCalculus 18Feb15 Matlab.zip]]
  • [get | view] (2016-02-18 15:52:39, 2.9 KB) [[attachment:Olaf IntroGLM2 17Feb2016.m]]
  • [get | view] (2015-02-04 11:58:19, 2.9 KB) [[attachment:Olaf IntroGLM2 4Feb15.m]]
  • [get | view] (2015-03-04 12:00:04, 113.0 KB) [[attachment:Olaf Optimisation 4Mar15.pdf]]
  • [get | view] (2015-01-16 13:56:32, 3.2 KB) [[attachment:Sampling.m]]
  • [get | view] (2015-01-16 13:56:14, 754.0 KB) [[attachment:Sampling.pdf]]
  • [get | view] (2016-02-04 17:19:23, 1107.2 KB) [[attachment:SignalSamplingNoise.pdf]]
  • [get | view] (2018-11-13 12:08:19, 1479.3 KB) [[attachment:Signal_sampling_noise_AT18.pdf]]
  • [get | view] (2018-01-11 09:31:52, 259.2 KB) [[attachment:optimisation.zip]]
 All files | Selected Files: delete move to page copy to page

You are not allowed to attach a file to this page.