% lecture 01
% Moo K. Chung mkchung@wisc.edu
% 
% CORTICAL THICKNESS MODELING


% loading *.mat file
load autism_cortical_surface.mat

% 1st subject, autism, inner surface
% need to reformat the data structure for visualization
surf.vertices=squeeze(autism_inner(1,:,:))';
surf.faces=tri;
figure_wire(surf, 'yellow');

figure
surf.vertices=squeeze(autism_outer(1,:,:))';
surf.faces=tri;
figure_wire(surf, 'blue');

% visualizing sclar data on surface mesh


figure_trimesh(surf,surf.vertices(:,2))

% cortical thickness for autism subject 1 
thickness=squeeze(L2norm(autism_outer(1,:,:) - autism_inner(1,:,:)));
figure_trimesh(surf,thickness)

% cortical thickness at vertex 10 for all subject
thickness1=L2norm(autism_outer(:,:,10) - autism_inner(:,:,10))
thickness2=L2norm(control_outer(:,:,10) - control_inner(:,:,10))

figure;
plot(thickness1,'or');
hold on;
plot(thickness2,'ob');

% box plots

help boxplot

figure;
subplot(1,2,1)
boxplot(thickness1)
xlabel('autism')
ylabel('thickness (mm)')
subplot(1,2,2)
boxplot(thickness2)
xlabel('control')
ylabel('thickness (mm)')
% two sample t-test

help ttest2

[h, p, ci, stats]=ttest2(thickness1, thickness2, 0.05,'both','unequal')

% tstat 0.7947

%--------------------------------------------------------
% general linar model 1:
% thickness = lambda1 + beta1*group
% we test for the significance of the coefficient beta1

thickness=[thickness1; thickness2]
group=[zeros(16,1); ones(11,1)]


%H_0: reduced model 
%     thickness = lambda_1 
n=length(thickness);
Z=ones(n,1)

lambda1 = inv(Z'*Z)*Z'*thickness

% for this model, the least squares estimation is the sample mean
mean(thickness)

SSE0 = sum((thickness - (Z*lambda1)).^2)


%H_1: full model
%     thickness = lambda1 + beta1*group

Z=[ones(n,1) group]

gamma = inv(Z'*Z)*Z'*thickness  % gamma = (lambda1, beta1)
% lambda1=3.81, beta1=-0.41

SSE1 = sum((thickness - (Z*gamma)).^2)

F=(SSE0-SSE1)./(SSE0/(n-1-1))

%F= 0.6384
% This is not exactly equal to Tstat^2. Slight difference is 
% expected since the T-statistic is approximately computed.

%pvalue

pvalue= 1-fcdf(F,1,n-1-1)
%pvalue =0.4318

%-----------------------------------------------------------
% general linar model 2:

% thickness = lambda1 + lambda2*age + beta1*group
% we test for the significance of the coefficient beta1
% while accountin for age effect.

age_control = [15 18 18 15 13 18 15 21 17 16 23]'
age_autism = [15 20 17 13 12 15 25 14 15 14 24 18 10 12 22 12]'
age=[age_control;age_autism]



%H_0: reduced model 
%     thickness = lambda_1 + lambda2*age

Z=[ones(n,1) age]
lambda = inv(Z'*Z)*Z'*thickness
% lambda1 =2.2046, lambda2=0.0872

SSE0 = sum((thickness - (Z*lambda)).^2)

%H_1: full model
%     thickness = lambda1 + lambda2*age + beta1*group

Z=[ones(n,1) age group]

gamma = inv(Z'*Z)*Z'*thickness
% lambda 1= 2.3830, lambda2=0.0861, beta1 = -0.3933

SSE1 = sum((thickness - (Z*gamma)).^2)

F=(SSE0-SSE1)./(SSE0/(n-2-1))

%F= 0.6384
% This is not exactly equal to Tstat^2. Slight difference is 
% expected since the T-statistic is approximately computed.

%pvalue

1-fcdf(F,1,n-2-1)
% pavlue = 0.4414

% HW3: Using the code developed in HW1, compute the surface area of both
% inner and couter cortical surfaces for all subjects. Then determine the
% significance of the group effect on cortical thickness while accounting 
% for total cortical surface area and age for all mesh vertices. What is 
% the maximum F statistic value and the corresponding p-value? 
% Provide MATLAB codes.


for i=1:40962
   % put some whatever codes here...
end;

figure_trimesh(surf,p)

figure_trimesh(surf,-log10(p))

figure_trimesh(surf,F)

% HW4 (ANCOVA): Set up the general linear model testing for interaction between the
% group and age for all mesh vertices. That is, you need to test the 
% significance of the coefficient beta1 in the model: 
%
% thickness = lambda1 + lambda2*age + lambda3*group + beta1*age*group
% 
% What vertex gives you the most significant interaction? How will you
% interprete your result?
%
%

%reduced model

thickness = lambda1 + lambda2*age + lambda3*group

% NOTE: Only registered students submit all 8 homework problems on October 9
%      1:00pm. 
