% Statistical and Computational Methods in Brain Image Analysis
% Chapter 18 Mixed Effect Models
%
% (C) Moo K. Chung 
% University of Wisconsin-Madison
% mkchung@wisc.edu


load hippocampus.mat

% The dataset is obtained from 206 MRI.
% The first 124 entries are cross sectional. 
% The remaining 82 entries are repeat scans. 

% group
% high income family: 0    low income family: 1


%-----------------------------------------------------------
% Laplace-Beltrami eigenfunction based heat kernel smoothing
[A, C] =FEM(hippoleft);
[Vl,Dl] = eigs(C,A, 500, 'sa');

% sigma=0.5 and k=500 degree used for hippocampus
hippolefts = lb_smooth([],hippoleft,0.5, 500, Vl, Dl);

for i=1:4
    figure;
    figure_trimesh(hippolefts,Vl(:,i),'rwb')
    colormap('hot')
    view([0 90]); camlight
    caxis([-0.04 0.04]);
end;


%---------------------------------------------------------
% FIXED EFFECT MODEL ON 124 SUBJECTS

Age=term(age);
Gender=term(gender);
Group=term(group);

% Change surface mesh to SURF-STAT format

leftsurf.coord= hippolefts.vertices';
leftsurf.tri=hippolefts.faces';

% length of displacement
left_disp2 = sqrt(sum(left_disp.^2,3));


% testing group difference
slm0 = SurfStatLinMod( left_disp2, 1+ Gender + Age + Group, leftsurf);
slm = SurfStatT(slm0, group); 
figure_origami(hippolefts,slm.t)
colormap('hot')

pvalue= [0.01 0.05 0.1]
threshold=randomfield_threshold(slm, pvalue)

%-----------------------------------
% MIXED EFFECT MODELS

Age=term(age);
Gender=term(gender);
Group=term(group);
Subject = term(var2fac(subject));


M1 = 1+ Gender + Age + Group + random(Subject) + I;
slm = SurfStatLinMod(left_disp2, M1, leftsurf);
slm = SurfStatT( slm, group);
figure_origami(hippolefts,slm.t)
colormap('hot')

pvalue= [0.0001 0.001 0.01 0.05 0.1]
threshold=randomfield_threshold(slm, pvalue)



