format compact
format long
echo on

% See Raz J.
% "A Statistical Critique of SPM
%    with Suggestions for Alternative Methods for fMRI Analysis" 
% available from http://www.sph.umich.edu/~jonraz/papers.html
%
% Data from "oddball" design fMRI experiment by Turetsky (Penn Psychiatry)
% Tones presented every two seconds
% Image taken every two seconds in between tone presentations
% Random presentation of "target" and "standard" tone
% Target probability 15%
% 210 trials per experiment; fMRI images in trials 11 through 210.

% conditionvec: an indicator (0 or 1) for standard or target condition
% for the current trial (col 1) and each of 9 previous trials (cols 2-10)

% Time series for significant voxels: voxel11_34_18  voxel17_31_19
% for non-significant voxels:         voxel14_40_18  voxel29_38_18 

% (Significance determined by SPM?)
pause
load('example.mat');
load('example2.mat');

size(conditionvec)
conditionvec(:,1)'
sum(conditionvec(:,1))
pause

size(voxel11_34_18)

pause

figure(1); newplot; hold off;
plot(1:200,voxel11_34_18,'b');
hold on;
plot(1:200,voxel17_31_19,'b');
plot(1:200,voxel14_40_18,'r');
plot(1:200,voxel29_38_18,'r');
pause

ones200=ones(200,1);
X=[ones200 conditionvec];


[B,BINT,R,RINT,STATS]=regress(voxel11_34_18,X);

% REGRESS Multiple linear regression using least squares.
%    b = REGRESS(y,X) returns the vector of regression coefficients, B.
%    Given the linear model: y = Xb, 
%    (X is an nxp matrix, y is the nx1 vector of observations.) 
%    [B,BINT,R,RINT,STATS] = REGRESS(y,X,alpha) uses the input, ALPHA
%    to calculate 100(1 - ALPHA) confidence intervals for B and the 
%    residual vector, R, in BINT and RINT respectively. 
%    The vector STATS contains the R-square statistic along with the F 
%    and p values for the regression.

% Note: here y = Xb is voxel11_34_18 = X*B
%       so B = X\voxel11_34_18

NB = X\voxel11_34_18;
[B NB]'

pause

figure(2); hold off; newplot; plot(1:10,B(2:11),'b+-')

STATS

pause

[B,BINT,R,RINT,STATS]=regress(voxel17_31_19,X);
hold on; plot(1:10,B(2:11),'g+-')
STATS

pause

[B,BINT,R,RINT,STATS]=regress(voxel14_40_18,X);
plot(1:10, B(2:11),'ro--')
STATS

pause

[B,BINT,R,RINT,STATS]=regress(voxel29_38_18,X);
plot(1:10,B(2:11),'ko--')

STATS