%Choosing working directory
cd ~/surfstat

clear all; close all; clc

surf = SurfStatReadSurf( 'av.obj' );

%Adding the midline mask
mask = SurfStatMaskCut( surf );

%Adding the brainstem mask
maskb = (mask & SurfStatROI([0; -16; -8], 20, surf ) == 0);

cd ~/surfstat/glimfiles

%Creating the matrix of data.
[SubjectNum, gender, age, ThickFileLeft, ThickFileRight]=textread('icbm_glim.csv', ' %s %s %f %s %s');

%To read and assign the actual thickness data
T1 = SurfStatReadData( [ThickFileLeft] );
T2 = SurfStatReadData( [ThickFileRight] );
T3 = horzcat(T1, T2);

%Create terms
Gender = term (gender);
Age = term (age);
%age_normalized = (age-mean(age))/std(age);  %if you want to normalize age
%Age = term(age_normalized);  %if you want to normalize age

%% 1 You are now ready to create your model and to estimate it

Y = 1 + Age + Gender;
figure(1)
image(Y);
slm = SurfStatLinMod( T3, Y, surf );

% MAIN EFFECT OF GENDER, DIRECTION 1 (MALE > FEMALE)

%To get your t statistic for group
contrast_direction1 = Gender.male - Gender.female;
slm_direction1 = SurfStatT ( slm, contrast_direction1);
figure(2)
SurfStatView ( slm_direction1.t.*maskb, surf, 'tmap males > females, removing age' );
set(gcf,'numbertitle','off','name', 'tmap males > females, removing age')
%SurfStatColLim([-12,12]);

%To get thresholded p values using Random Field Theory
[ pval, peak, clus ] = SurfStatP( slm_direction1, maskb );
figure(3)
SurfStatView( pval, surf, 'RFT map males > females, removing age');
set(gcf,'numbertitle','off','name', 'RFT map males > females, removing age')

%To get thresholded p values using False Discovery Rate
qval = SurfStatQ( slm_direction1, maskb );
figure(4)
SurfStatView( qval, surf, 'FDR map males > females, removing age');
set(gcf,'numbertitle','off','name', 'FDR map males > females, removing age')

% MAIN EFFECT OF GENDER, DIRECTION 2 (FEMALE > MALE)

%To get your t statistic for group
contrast_direction2 = Gender.female - Gender.male;
slm_direction2 = SurfStatT ( slm, contrast_direction2);
figure(5)
SurfStatView ( slm_direction2.t.*maskb, surf, 'tmap females > males, removing age' );
set(gcf,'numbertitle','off','name', 'tmap females > males, removing age')
%SurfStatColLim([-12,12]);

%To get thresholded p values using Random Field Theory
[ pval, peak, clus ] = SurfStatP( slm_direction2, maskb );
figure(6)
SurfStatView( pval, surf, 'RFT map females > males, removing age');
set(gcf,'numbertitle','off','name', 'RFT map females > males, removing age')

%To get thresholded p values using False Discovery Rate
qval = SurfStatQ( slm_direction2, maskb );
figure(7)
SurfStatView( qval, surf, 'FDR map females > males, removing age');
set(gcf,'numbertitle','off','name', 'FDR map females > males, removing age')

% MAIN (NEGATIVE) EFFECT OF AGE

slm = SurfStatT ( slm, -age);
figure(8)
SurfStatView ( slm.t.*maskb, surf, 'tmap for -age, removing gender' );
set(gcf,'numbertitle','off','name', 'tmap for -age, removing gender')
%SurfStatColLim([-12,12]);

%To get thresholded p values using Random Field Theory
[ pval, peak, clus ] = SurfStatP( slm, maskb );
figure(9)
SurfStatView( pval, surf, 'RFT map for -age, removing gender');
set(gcf,'numbertitle','off','name', 'RFT map for -age, removing gender')

%To get thresholded p values using False Discovery Rate
qval = SurfStatQ( slm, maskb );
figure(10)
SurfStatView( qval, surf, 'FDR map for -age, removing gender');
set(gcf,'numbertitle','off','name', 'FDR map for -age, removing gender')

%Mean cortical thickness
figure(11)
SurfStatView( mean(T3), surf, 'mean cortical thickness');
set(gcf,'numbertitle','off','name','mean cortical thickness')
SurfStatColLim([0,5.5]);

%Standard deviation of cortical thickness
figure(12)
SurfStatView( std(T3), surf, 'std cortical thickness');
set(gcf,'numbertitle','off','name','std cortical thickness')
SurfStatColLim([0,.7]);

Back to Statistical Analyses Using SurfStat

CIVET Home