% You need SPM on your search path

% Get data for one voxel
load voxel_15_52
X=large_X;


% Which models to you want to compare ?
%comparison='AllwTemp_vs_mean';
%comparison='AllwTemp_vs_All';
comparison='All_vs_U1';
%comparison='U1_vs_mean';
%comparison='All_vs_mean';

switch comparison
    case 'U1_vs_mean',
        Xfull=[X(:,1),X(:,9)];
        Xred=X(:,9);
    case 'All_vs_U1',
        Xfull=[X(:,1),X(:,3),X(:,5),X(:,7),X(:,9)];
        Xred=[X(:,1),X(:,9)];
    case 'All_vs_mean',
        Xfull=[X(:,1),X(:,3),X(:,5),X(:,7),X(:,9)];
        Xred=[X(:,9)];
    case 'AllwTemp_vs_All',
        Xfull=X;
        Xred=[X(:,1),X(:,3),X(:,5),X(:,7),X(:,9)];
    case 'AllwTemp_vs_mean',
        Xfull=X;
        Xred=X(:,9);
end
       

% Number of parameters in designs
pfull=size(Xfull,2);
pred=size(Xred,2);

% Projection matrices
Pfull=Xfull*pinv(Xfull'*Xfull)*Xfull';
Pred=Xred*pinv(Xred'*Xred)*Xred';

% Projections (ie. model predictions)
yfull=Pfull*y;
yred=Pred*y;

% Residual forming matrices
N=size(y,1);
Rfull=eye(N)-Pfull;
Rred=eye(N)-Pred;

% Sums of squares
RSS_red=y'*Rred*y
RSS_full=y'*Rfull*y
ESS=RSS_red-RSS_full

% Alternative formula is 
% ESS=trace(y'*(Rred-Rfull)*y)

% Ranks
rESS=rank(Rred-Rfull)
rRSS=rank(Rfull)

% F statistic
F=(ESS/RSS_full)*(rRSS/rESS)
p=1-spm_Fcdf(F,rESS,rRSS)

disp('Square of (partial) correlation coefficient:');
disp(' Of the variability unexplained by the reduced model');
disp(' the proportion explained by the full model is');
R=(RSS_red-RSS_full)/RSS_red
disp('Correlation')
r=sqrt(R)
F=(R/(1-R))*(rRSS/rESS)

% Plot models
figure
subplot(2,2,2);
title('Full model');
imagesc(Xfull);
ylabel('Scans');
set(gca,'XTicklabel',[]);
%xlabel('Regressors');
title('X');

colormap('gray');
subplot(2,2,4);
title('Reduced model');
imagesc(Xred);
title('X0');
colormap('gray');
ylabel('Scans');
set(gca,'XTicklabel',[]);
xlabel('Regressors');

% Plot time series and full model fit
subplot(2,2,1);
plot(y);
hold on
plot(yfull,'r');
%title(sprintf('Time series at voxel x=%d, y=%d',xx,yy));
title(sprintf('R^2=%1.2f',R));
%legend('Data','Full model');
ylabel('fMRI');
%xlabel('Scans');

% Plot time series and reduced model fit
subplot(2,2,3);
plot(y);
hold on
plot(yred,'r');
%title(sprintf('Time series at voxel x=%d, y=%d',xx,yy));
%legend('Data','Reduced model');
ylabel('fMRI');
xlabel('Scans');