## Contents

```% This matlab file is a companion to the
% 2 lectures on the General Linear Model
%
% Cyril Pernet January/February 2011
```

## 1. Linearity

Definition

1. the quality or state of being linear 2. Electronics: the extent to which any signal modification process, as detection, is accomplished without amplitude distortion 3. Physics: the extent to which any effect is exactly proportional to its cause

TWO KEY NOTIONS

additivity (or superposition): the total response to a set of inputs is the sum of individual inputs (y = x1+x2)

scaling (or homogeneity of degree 1): the magnitude of the system output is proportional to the system input (y = ax)

```% ---------------------
% Linear correlations
% ---------------------

clear all; clc

% create 2 variables and correlate
x1 = [-15:14]' + randn(30,1);
x2 = 3*x1;
r = corr(x1,x2,'type','Pearson');
figure('Name','Correlations')
subplot(3,2,1); title('Linear correlation');
plot([x1 x2],'LineWidth',3); grid on;
title(['Scalling effect: r = ' num2str(r)]);

% Superposition
x3 = x1 + x2;
r = corr(x1,x3,'type','Pearson');
subplot(3,2,2);
plot([x1 x3],'LineWidth',3); grid on;
title(['Superposition effect: r = ' num2str(r)]);

% Scaling and Superposition
x4 = x2 + 2*x3 + 5; % = 11*x1 + 5
r = corr(x1,x4,'type','Pearson');
subplot(3,2,3);
plot([x1 x4],'LineWidth',3); grid on;
title(['Scaling and Superposition effect: r = ' num2str(r)]);

% Non linear ^2
x5 = x1.^2;
r = corr(x1,x5,'type','Pearson');
subplot(3,2,4);
plot([x1 x5],'LineWidth',3); grid on;
title(['Squaring effect: r = ' num2str(r)]);

% non linear abs
x6 = abs(x1);
r = corr(x1,x6,'type','Pearson');
subplot(3,2,5);
plot([x1 x6],'LineWidth',3); grid on;
title(['Absolute value effect: r = ' num2str(r)]);

% non linear but ..
x7 = 3*x1 + x1.^2 + x1.^3;
r = corr(x1,x7,'type','Pearson');
subplot(3,2,6);
plot([x1 x7],'LineWidth',3); grid on;
title(['polynomial effect: r = ' num2str(r)]);  % yes polynomes are special cases

% exercise: try a different argument 'type' for the corr function
```

## 2. Linear algebra

with linear algebra one uses matrices to write equations as we shall see this turns out to be very usefull as it can solve very quickly complex linear systems

a simple linear system with two variables can be 2x - y = 0 -x + 2y = 3

in matrix notation this is

```clc; clear; close

X = [2 -1; -1 2];
Y = [0 3]';

% we can know apply the same principle as solving a single equation
% for instance y = bx is solved by finding b = 1/x * y
% with matrices we also use the inverse, except that the inverse of
% a matrix is a bit more complicated

I = eye(2);  % this is called the identity matrix because X*I = I*X = X (just as x*1 = 1*x = x)
Xinv = inv(X);  % this is the inverse of A defined as X*Xinv = I jus as 1/x * x = 1
test = (Xinv*X == I);

% the solution of our linear system is therefore
B = inv(X)*Y; % just as with numbers b = 1/x * y

% which we can check
XB = X*B; % = Y

% in geometrical terms, each column of a matrix describes a vector - the
% solution is to find how much (B) of each vector (X(:,i)) we need to arrive to the point Y
```

## 3. Simple linear regression

In a regression we have some observations Y that we want to explain by a variable X by contrast with a linear system as describe above, there is usually no perfect solution the model is y = b1x + b2 + e with b1 how much of x we need to get close to y, b2 the intercept or offset and e the error that is the part that is not explained by x. However, we can see that a regression is a linear model since y is expressed as a function of x using additivity and scaling. The statisitcs is simply to evaluate how much we explain by x given the error.

```clc; clear

% y = b1x + b2 + e
% let generate some data
x = rand(20,1)*5; % x is the VI we manipulate
Y = 4*x + 3 + randn(20,1); % Y is the VD we measure, we run here a forward model + add noise

% get a solution, ie find b1=4 and b2=3
X = [x ones(20,1)]; % ones are for b2 which is not a function of x but a constant to add for each x

% X is not squared so inv(X) is not possible
% a simple solution is to multiply by the transpose
% making the matrix X'X square
% Y = XB --> X'Y = X'XB --> inv(X'X)X'Y = inv(X'X)X'XB --> but inv(X'X)X'X = I thus

B = inv(X'*X)*X'*Y;

% check what it looks like
Yhat = X*B; % Yhat is the model
Res  = Y - Yhat; % Residual (error)

% get some statistics
% = estimate if the model Yhat explain above the error Res

SStotal  = norm(Y-mean(Y)).^2; % what is the norm of a vector? = sum((Y-mean(Y)).^2)
SSeffect = norm(Yhat-mean(Yhat)).^2; % just like in books, sum of the difference xi - mean(x)
SSerror  = norm(Res-mean(Res)).^2;
df       = rank(X)-1; % what is the rank ? see also http://en.wikipedia.org/wiki/Linear_independence
dferror  = length(Y) - df - 1;

R2 = SSeffect / SStotal; % how much we explain of the total
F  = (SSeffect / df) / (SSerror / dferror); % how much is explain given the error, this is the same as R2*dfe / (1-R2)*df
p  = 1 - fcdf(F,df,dferror);

% make a figure
figure('Name','Simple Regression');
subplot(1,2,1); plot(Y,'x','LineWidth',5);  hold on
plot(Yhat,'r','LineWidth',2); grid on;
title(['Data and model F=' num2str(F) ' p=' num2str(p)])
subplot(1,2,2); normplot(Res);
```

## 4. Multiple linear regression

This is the same as above except that we have more variables - however, using linear alebra we can find as easily the regression coefficients. In addition to testing the model, we usually also test for the weight of each regressor. The relative contribution of 1 regressor to the data given the model is the partial correlation coefficient (how much x1 explains of y when we control for the influence of all other xs on both y and x1) and the relative contribution of 1 regressor to the total variance is the semi-partial correlation coefficient (how much x1 explains of y when we control for the influence of all other xs on x1 only)

```clc; clear; close

% y = b1x1 + b2x2 + b3x3 + b4x4 + b5
% let's take data stored in matlab exemples
Y=hald(:,5) ;
X=hald(:,1:4);

% get a solution,
X = [X ones(length(X),1)];
B    = inv(X'*X)*X'*Y; % just as above
Yhat = X*B;
Res  = Y - Yhat;

SStotal = norm(Y - mean(Y)).^2;
SSeffect = norm(Yhat - mean(Yhat)).^2;
SSerror  = norm(Res-mean(Res)).^2;

df      = rank(X)-1;
dferror = length(Y) - df - 1;
R2      = SSeffect / SStotal;
F       = (SSeffect / df) / (SSerror / dferror);
p       = 1 - fcdf(F,df,dferror);

% make a figure
figure('Name','Multiple Regression');
subplot(1,2,1); plot(Y,'x','LineWidth',5);  hold on
plot(Yhat,'r','LineWidth',2); grid on;
title(['R^2 ' num2str(R2) ' Data and model F=' num2str(F) ' p=' num2str(p)],'FontSize',14)
subplot(1,2,2); normplot(Res);

% ------------------------------------
% semi-partial correlation coefficient
% ------------------------------------
% let's think of the model and what it means it terms of geometry.
% the data Y can be described as a vector and a point in R_20
% a space with 20 dimensions. We then establish a model X with 5
% regressors; that is we look for a combination of these 5 vectors which
% will get as close as possible to Y. To find the actual contribution of x1
% to the data for this model one needs to look at how much x1 explains
% to the total variance, ie we want to compare the R2 between the full and
% a reduced model without x1 - the difference will be how much x1 explains in Y.

Xreduced    = X(:,2:end); % reduced model all minus 1st regressor
Breduced    = inv(Xreduced'*Xreduced)*Xreduced'*Y;
Yhatreduced = Xreduced*Breduced;
Resreduced  = Y - Yhatreduced;

dfreduced       = rank(Xreduced) -1 ;
dferrorreduced = length(Y) - dfreduced - 1;
SSeffectreduced = norm(Yhatreduced-mean(Yhatreduced)).^2;
SSerrorreduced  = norm(Resreduced-mean(Resreduced)).^2;
R2reduced       = SSeffectreduced / SStotal;
Freduced       = (SSeffectreduced / dfreduced) / (SSerrorreduced / dferrorreduced);
preduced       = 1 - fcdf(Freduced,dfreduced,dferrorreduced);

Semi_Partial_corr_coef = R2 - R2reduced;
dfe_semi_partial_coef  = df - dfreduced;
F_semi_partail_coef    = (Semi_Partial_corr_coef*dferror) / ...  % variance explained by x1
((1-R2)*dfe_semi_partial_coef); % unexplained variance overall
p_semi_partial_coef    = 1 - fcdf(Semi_Partial_corr_coef, df, dfe_semi_partial_coef); % note df is from the full model

% make a figure
figure('Name','Multiple Regression - Full versus reduced model');
subplot(2,2,1); plot(Y,'x','LineWidth',5);  hold on
plot(Yhat,'r','LineWidth',2); grid on;
title(['Data and model F=' num2str(F) ' p=' num2str(p)])
subplot(2,2,2); plot(Y,'x','LineWidth',5);  hold on
plot(Yhatreduced,'r','LineWidth',2); grid on;
title(['Data and reduced model F=' num2str(Freduced) ' p=' num2str(preduced)])
subplot(2,2,3); plot(Yhat,'b','LineWidth',2); grid on; hold on
plot(Yhatreduced,'r','LineWidth',2);
title('Modelled data for both models')
subplot(2,2,4); plot((Res-Resreduced).^2,'k','LineWidth',2); grid on; hold on
title('Difference of residuals squared')

% --------------------------------
% partial correlation coefficient
% --------------------------------
% As we shall see below, this is easily obtained using projections - but
% for now let just think about what we want to measure. We are interested
% in knowing the correlation between y and x1 controlling for the effect of
% the other xs, that is removing the effect of other xs. Compred to
% semi-parital coef we also want to remove the effect of xs on y or if you
% prefer we want to compute how much of x1 we need to get to the point
% defined by the xs which is the closest to Y

% above we removed X(:,2:end) from Y and got Resreduced
% we need to do the same for x1 so that we can correlate x1 and y witout xs
x = X(:,1);
B = inv(Xreduced'*Xreduced)*Xreduced'*x;
xhat = Xreduced*B;
Resx = x - xhat;

% the correlation between Resreduced and Resx is the partial coef
Partial_coef = corr(Resx,Resreduced);

% --------------------------------------------------
```

## 5. Matlab tools for regressions (from stat toolbox)

-------------------------------------------------

```clc; close
[B,Bint,r,rint,stats] = regress(Y,X); % Removes NaN data

% B are the beta parameters exactly what we have abbove with inv(X'*X)*X'*Y
% B int are the 95% confidence intervals for the coefficient
% r are the residuals - same as above
% r int are the 95% confidence intervals of residuals - note it should
% include 0 otherwise would suggest that a point is not well fitted (outlier?)
% this can be checked with
rcoplot(r,rint);
% stats: R2 statistic, the F statistic and its p value, and an estimate of the error variance.

regstats(Y,X(:,2:end)); % call a GUI - WARNING NO INTERCEPT SHOULD BE USED
% stats = regstats(y,X,model,whichstats) does not call the gui and returns
% in stats the paramters specified in whichstats

[br,statsr] = robustfit(X(:,2:end),Y); % robust estimate using iterative weighted least square

% GLMFIT works for regression and ANOVAs
[b,dev,stats] = glmfit(X(:,2:end),Y,'normal'); % pretty cool feature is that you can have other types of distributions

% shame none of the above do partial coef and semi-partial coef.

% ------------------
```
```Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.071144e-17.
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.071144e-17.
Warning: X is rank deficient, rank = 4
Warning: X is ill conditioned, or the model is
overparameterized, and
some coefficients are not identifiable.  You should use
caution
in making predictions.
```

## 6. One-way ANOVA

-----------------

We now describe a one-way ANOVA as a linear model, just as we did with regression - the standard model is y = u + xi + e that is the data are equal to the mean + effect of treatment + error. The effect of the treatment for each condition is suppose to be the same for each y in a given condition - this can be thus be modelled using 1s.

```clc; clear; close all

% Y = u1 + u2 + u3

u1 = rand(10,1) + 11.5;
u2 = rand(10,1) + 7.2;
u3 = rand(10,1) + 5;
Y = [u1; u2; u3];

% describe X as 3 conditions
x1 =[ones(10,1); zeros(20,1)]; % condition 1 as a given value for y(1:10)
x2 =[zeros(10,1); ones(10,1); zeros(10,1)]; % condition 2 as a given value for y(11:20)
x3 =[zeros(20,1); ones(10,1)]; % condition 3 as a given value for y(21:30)
X =[x1 x2 x3 ones(30,1)]; % again we add the intercept here the grand mean

% now we can do the almost exactly same analysis as above
B = pinv(X)*Y; % almost as above

% here we use a pseudoinverse:
% X is rank deficient, i.e. regressors are not independent, since any linear combination of
% 3 columns can give us the 4th one, thus X'X is also rank deficient and singular ie inv(X'*X)
% doesn't exist -- there is no matrix A such as A(X'X) = I - however pinv gives a unique
% solution that minimizes the square distance to the data.
% see
% http://mathworld.wolfram.com/Moore-PenroseMatrixInverse.html
% http://en.wikipedia.org/wiki/Moore_Penrose_pseudoinverse

Yhat = X*B;
Res = Y - Yhat;
SStotal = norm(Y-mean(Y)).^2;
SSeffect = norm(Yhat-mean(Yhat)).^2;
SSerror  = norm(Res-mean(Res)).^2;
df = rank(X)-1;
dferror = length(Y) - df - 1;
R2 = SSeffect / SStotal;
F = (SSeffect / df) / (SSerror / dferror);
p = 1 - fcdf(F,df,dferror);

% make a figure
figure('Name','One way ANOVA');
subplot(1,3,1); plot(Y,'x','LineWidth',5);  hold on
plot(Yhat,'r','LineWidth',2); grid on;
title(['Data and model F=' num2str(F) ' p=' num2str(p)])
subplot(1,3,2); imagesc(X); colormap('gray'); title('Design matrix');
subplot(1,3,3); normplot(Res);

% for comparison try [P,ANOVATAB,STATS] = ANOVA1(Y,g,'off')
% with g = [ones(30,1)]; g(11:20)=2; g(21:30)=3;

% ------------------------
```

## 7. Projection matrices

------------------------

```% see  http://en.wikipedia.org/wiki/Projection_(linear_algebra)

% P is the orthogonal projection of Y onto X sch as Y^=P*Y
% the projection is orthoganal which means this minimize (in the least
% square sense) the distance between Y and X
% P = Xpinv(X) and Y^ = Xpinv(X)Y
% which is the same as above with B = pinv(X'X)X'Y and Y^ = XB

% -------------------
```

## 8. way ANOVA again

-------------------

```% following the results obtained above
P = X*pinv(X); % orthogonal projection of Y onto X
R = eye(size(Y,1)) - P; % projection onto the error space null(X)
E = Y'*R*Y; % SS error
C = diag([1 1 1 0]); % test effect of treatement (thus remove last column of X)
C0 = eye(size(X,2)) - C*pinv(C);
X0 = X*C0; % Reduced model (i.e. only intercept)
R0 = eye(size(Y,1)) - (X0*pinv(X0)); % projection onto null(C(X))
M = R0 - R; % M is the projection matrix onto Xc, what we want
H = (B'*X'*M*X*B); % SSeffect
F_Rsquare = (H./(rank(X)-1)) ./ (E/(size(Y,1)-rank(X)));
p_Rsquare = 1 - fcdf(F_Rsquare, (rank(X)-1), (size(Y,1)-rank(X)));

% ex. use a difference C to compute gp 1 vs gp 2
% C = diag([1 -1 0 0])

% ------------------------
```

## 9. repeated measure ANOVA

------------------------ Y = XB + E with X = [Factors Interactions / Subject]

```clc; clear; close all

% Y = u1 + u2 + u3

u1 = rand(10,1) + 11.5;
u2 = rand(10,1) + 7.2;
u3 = rand(10,1) + 5;
Y = [u1; u2; u3]; % this time this not 3 gp but 10 subjects with 3 measures

% create the design matrix for the different factors
nb_subjects =10;
nb_conditions =3;
Subjects = repmat(eye(nb_subjects),nb_conditions,1); % error
x = kron(eye(nb_conditions),ones(nb_subjects,1));  % effect
X = [x Subjects]; % no more ones for the grand mean but a subject specific mean
figure; imagesc(X); colormap('gray'); title('Repearted measure design','Fontsize',14)

% Compute as usual
df  = nb_conditions -1;
dfe = size(Y,1)  - nb_subjects - df;

P     = X*pinv(X'*X)*X'; % our projection matrix
R     = eye(size(Y,1)) - P; % projection on error space
SSe   = diag(Y'*R*Y); % Y projected onto the error
Betas = pinv(x)*Y;  % compute without cst/subjects
yhat  = x*Betas; % yhat computed based on the treatment with subject effect - we use little x
SS    = norm(yhat-mean(yhat)).^2;
F_values = (SS/df) ./ (SSe/dfe);
p_values = 1 - fcdf(F_values, df, dfe);

% ------------------------
```

## 10. multivariate regression

------------------------

```clc; clear; close all

% y = b1x1 + b5
Y=[Acceleration Displacement];% now we have 2 Ys

% get B
X = [Cylinders Weight ones(length(Weight),1)];
B    = inv(X'*X)*X'*Y; % just as before but now we have 2 columns, 1 per Y :-)

% use projection
T        = (Y-repmat(mean(Y),size(Y,1),1))'*(Y-repmat(mean(Y),size(Y,1),1));  % SSCP Total
P        = X*pinv(X'*X)*X';
R        = eye(size(Y,1)) - P;
E        = (Y'*R*Y); % SSCP Error

% stats for the model
C = eye(size(X,2));
C(:,size(X,2)) = 0;
C0   = eye(size(X,2)) - C*pinv(C);
X0   = X*C0;
R0   = eye(size(Y,1)) - (X0*pinv(X0));
M    = R0 - R;
H    = (B'*X'*M*X*B);

% if we were doing univariate stats the SS are on the diagonal of the SSCP
% matrices so simply use diag
Rsquare   = diag(H)./diag(T);
F_Rsquare = (diag(H)./(rank(X)-1)) ./ (diag(E)/(size(Y,1)-rank(X)));
p_Rsquare = 1 - fcdf(F_Rsquare, (rank(X)-1), (size(Y,1)-rank(X)));

% what about multivariate regression then
% we decompose the product of the SSCP effect * inv(SSCP error)
% an eigen value decomposition is like a PCA, in matrix term we
% look for a new set of vectors orthogonal to each other but their
% combination gives the same result (in space) than the original
% matrix see http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix

eigen_values = eig(inv(E)*H);  % for real life processing better to use SVD to get positive results only

p = size(Y,2); % = number of variables (dimension)
q = size(X,2)-1; % = number of regressors (df)
n = size(Y,1); % nb of observations (dfe)
s = min(p,q);
m = (abs(q-p)-1)/2;
N = (n-q-p-2)/2;
d = max(p,q);

theta = max(eigen_values) / (1+max(eigen_values)); % Roy
V = sum(eigen_values ./ (1+eigen_values)); % Pillai

R2_Roy_value = theta; % = 1st canonical correlation
R2_Roy_F     = ((n-d-1)*max(eigen_values))/d;
R2_Roy_p     = 1-fcdf(R2_Roy_F, d, (n-d-1));

R2_Pillai_value = V / s; % average of canonical correlations
R2_Pillai_F     = ((2*N+s+1)*V) / ((2*m+s+1)*(s-V)');
R2_Pillai_p     = 1-fcdf(R2_Pillai_F,(s*(2*m+s+1)),(s*(2*N+s+1)));

% what about the effect of 1 regressor in particular? (partial coef)
% let's look at the 2nd regressor
C    = diag([0 1 0]);
C0   = eye(size(X,2)) - C*pinv(C);
X0   = X*C0;
R0   = eye(size(Y,1)) - (X0*pinv(X0));
M    = R0 - R;
H    = B'*X'*M*X*B;
eigen_values = eig(inv(E)*H);

df_multivariate = size(Y,2); % nb of DV
dfe_multivariate = size(Y,1)-size(Y,2)-(size(X,2)-1); % N-df-nb of covariates

% Roy's test
theta = max(eigen_values) / (1+max(eigen_values));
% Pillai
V = sum(eigen_values ./ (1+eigen_values));
% Lawley-Hotelling's generalized T2
U_continuous = sum(eigen_values);

F_continuous= (dfe_multivariate*max(eigen_values))/df_multivariate;
pval = 1-fcdf(F_continuous, df_multivariate, dfe_multivariate);
```