% Statistical and Computational Methods in Brain Image Analysis
% Chapter 7. General Multivariate Linear Models
%
% (C) Moo K. Chung
% University of Wisconsin-Madison
% mkchung@wisc.edu
%
% May 1, 2019 updated

%-----------------------------------------------------
% Simulating multivariate normal distribution

V=[2 1
     1 2];

Htr=chol(V)

Htr'*Htr


z1 = normrnd(0,1,1000,1); 
z2=normrnd(0,1,1000,1); 
w=Htr'*[z1 z2]';

figure; plot(w(1,:), w(2,:),'.k')


hatV = cov(w')
hatH=chol(hatV)'
newz=inv(hatH)*w;

figure; plot(newz(1,:), newz(2,:) ,'.k')

figure_bg('w')
print_pdf('pp')

figure; qqplot(newz(1,:));
figure; qqplot(newz(2,:));

%----------------------------------------
% Covariance matrix factorization

V=[2 1
     1 2]
 
[Q D]=eig(V)

Q*D.^(1/2)*Q'

V^(1/2)

w=V^(1/2)*[z1 z2]'


%------------------------------  
% Visualizing amygdala and hippocampus template for MIDUS data set, which
% was published in
%
% [1] Chung, M.K., Schaefer, S.M., van Reekum, C. M., Schmitz, L.P., Sutterer, M., 
% Davidson, R.J. 2014. A unified kernel regression for diffusion wavelets on 
% manifolds detects aging-related changes in the amygdala and hippocampus, 
% International Conference on Medical Image Computing and Computer Assisted Intervention 
% (MICCAI). Lecture Notes in Computer Science(LNCS) 8674:791-298
% http://www.stat.wisc.edu/~chung/papers/chung.2014.MICCAI.pdf
%
% [2] Chung, M.K., Kim, S.-G., Schaefer, S.M., van Reekum, C. M., Schmitz, L.P., 
% Sutterer, M., Davidson, R.J. 2014. Improved statistical power with a sparse 
% shape model in detecting an aging effect in the hippocampus and amygdala, 
% Proceedings of SPIE Medical Imaging 9034, 90340Y
%
% If you are using this dataset or codes, please refernce [1] or [2] above.

load midus.mat

amygl=template{1,1}
amygr=template{1,2}
hippol=template{2,1}
hippor=template{2,2}

figure; figure_wire(amygl,'black','red');
hold on; figure_wire(amygr,'black','blue');
hold on; figure_wire(hippol,'black','yellow');
hold on; figure_wire(hippor,'black','green');
view([-120 15])

% Visualizing displacement vector field for the 10th subject

dispal= displacement{10,1,1};
dispar= displacement{10,1,2};
disphl= displacement{10,2,1};
disphr= displacement{10,2,2};

hold on;
hold on; figure_quiver2(amygl, dispal, 10, 'red')
hold on; figure_quiver2(amygr, dispar, 10, 'blue')
hold on; figure_quiver2(hippol, disphl, 10, 'yellow')
hold on; figure_quiver2(hippor, disphr, 10, 'green')

legend('Left amygdala', 'Right amygala', 'Left hippocampus', 'Right hippocampus')
legend('boxoff')
view([-120 15])
axis off


%==========================
% Keith Worsley's SusfStat Package. You must d/l and install the package
% from http://www.math.mcgill.ca/keith/surfstat (Septembrer 26, 2008 version)
% http://galton.uchicago.edu/faculty/InMemoriam/worsley/research/surfstat (Febuary 21, 2009 version)


% Data reformatting
% It is necessary to reformat the data into a 3D array of subjects x mesh vertices x coordinates
% to be used in SurfStat

dispal=[];
for i=1:69
    dispal(i,:,:)=displacement{i,1,1} ;
end;

dispar=[];
for i=1:69
    dispar(i,:,:)=displacement{i,1,2};
end;

disphl=[];
for i=1:69
    disphl(i,:,:)=displacement{i,2,1};
end;

disphr=[];
for i=1:69
    disphr(i,:,:)=displacement{i,2,2};
end;
   

% T-test in SurfStat

dispal2 = sqrt(sum(dispal.^2,3));
dispar2 = sqrt(sum(dispar.^2,3));
disphl2 = sqrt(sum(disphl.^2,3));
disphr2 = sqrt(sum(disphr.^2,3));

Gender=term(gender);
Age=term(age);
Brain=term(brain);


slm1 = SurfStatLinMod( dispal2, 1+Brain + Gender+Age);
slm = SurfStatT( slm1, gender);
figure; figure_trimesh( amygl ,slm.t)

slm1 = SurfStatLinMod( dispar2, 1+Brain + Gender+Age);
slm = SurfStatT( slm1, gender);
hold on; figure_trimesh( amygr ,slm.t)

slm1 = SurfStatLinMod( disphl2, 1+Brain + Gender+Age);
slm = SurfStatT( slm1, gender);
hold on; figure_trimesh( hippol ,slm.t)

slm1 = SurfStatLinMod( disphr2, 1+Brain + Gender+Age);
slm = SurfStatT( slm1, gender);
hold on; figure_trimesh( hippor ,slm.t)
view([-120 15])

%print -dpdf -r300 midus.tstat.pdf
% F-test in SurfStat

slm0 = SurfStatLinMod( dispal2, 1+Brain + Age);
slm1 = SurfStatLinMod( dispal2, 1+Brain + Gender+Age);
slm = SurfStatF(slm1, slm0); 
figure; figure_trimesh( amygl ,slm.t)

slm0 = SurfStatLinMod( dispar2, 1+Brain +Age);
slm1 = SurfStatLinMod( dispar2, 1+Brain + Gender+Age);
slm = SurfStatF( slm1, slm0);
hold on; figure_trimesh( amygr ,slm.t)

slm0 = SurfStatLinMod( disphl2, 1+Brain +Age);
slm1 = SurfStatLinMod( disphl2, 1+Brain + Gender+Age);
slm = SurfStatF( slm1, slm0);
hold on; figure_trimesh( hippol ,slm.t)

slm0 = SurfStatLinMod( disphr2, 1+Brain +Age);
slm1 = SurfStatLinMod( disphr2, 1+Brain + Gender+Age);
slm = SurfStatF( slm1, slm0);
hold on; figure_trimesh( hippor ,slm.t)
view([-120 15])

%print -dpdf -r300 midus.Fstat.pdf


% Hotelling's T2 statistic
% testing the effect of gender
% Compare this to Hotelling's T2 routine given in
% http://pages.stat.wisc.edu/~mchung/research/amygdala/

dispalm = dispal(find(~gender),:,:);
dispalf = dispal(find(gender),:,:);
h =SurfStat2HotelT2(dispalm, dispalf);
figure; figure_trimesh( amygl ,h.t)

disparm = dispar(find(~gender),:,:);
disparf = dispar(find(gender),:,:);
h =SurfStat2HotelT2(disparm, disparf);
hold on; figure_trimesh( amygr ,h.t)

disphlm = disphl(find(~gender),:,:);
disphlf = disphl(find(gender),:,:);
h =SurfStat2HotelT2(disphlm, disphlf);
hold on; figure_trimesh( hippol ,h.t)

disphrm = disphr(find(~gender),:,:);
disphrf = disphr(find(gender),:,:);
h =SurfStat2HotelT2(disphrm, disphrf);
hold on; figure_trimesh( hippor ,h.t)
view([-120 15])


% MGLM (T-stat) in SurfStat
% testing the effect of gender while accouting for age and brain size
% differences.

slm0 = SurfStatLinMod( dispal, Age + Brain + Gender);
slm = SurfStatT( slm0, gender);
figure; figure_trimesh( amygl ,slm.t)

slm0 = SurfStatLinMod( dispar, Age + Brain + Gender);
slm = SurfStatT( slm0, gender);
hold on; figure_trimesh( amygr ,slm.t)

slm0 = SurfStatLinMod( disphl, Age + Brain + Gender);
slm = SurfStatT( slm0, gender);
hold on; figure_trimesh( hippol ,slm.t)

slm0 = SurfStatLinMod( disphr, Age + Brain + Gender);
slm = SurfStatT( slm0, gender);
hold on; figure_trimesh( hippor ,slm.t)
view([-120 15])



