format short
format compact
clc
echo on
% Artificial data for two groups
NN = [10 15]; % number per group
NCases = sum(NN); NTreatments = 2;
MM = [5 5]; % means -- no difference!
SD = [1 1]; % standard deviation
Sample = [MM(1)+SD(1)*randn(NN(1),1); MM(2)+SD(2)*randn(NN(2),1)];
%
pause
clc
% Now calculate the various sums of squares:
%
% Individual -- 
II = Sample'*Sample
%
% Treatment --
T1 = Sample(1:NN(1)); T2 = Sample((NN(1)+1):NCases);
TT = (sum(T1)^2)/NN(1) + (sum(T2)^2)/NN(2)
%
% Combination --
CC = (sum(Sample)^2)/length(Sample)
%
pause
clc
% Now calculate the partition of variance...
% Total square deviations: II - CC with NCases-1 degrees of freedom
Total = II - CC   % df = 34
%
% Treatment: TT - CC with NTreatments-1 degrees of freedom
Treatment = TT - CC   % df = 1
%
% Residual: II - TT with  NCases - NTreatments degrees of freedom
Residual = II - TT   % df = NCases - NTreatments = 33
% Do they add up?
Treatment + Residual
pause
clc
% Now calculate the F ratio:
MSTR = Treatment/(NTreatments-1)  % df =1
MSE = Residual/(NCases-NTreatments) % df = 34
F = MSTR/MSE  % Should be small, since there is really no difference
pause
clc
% Now compare against expectations for F(1,34)
F
finv(1-[.05 .01 .005 .001],1,34) % significance thresholds for .05, .01, etc.
pause
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OK, now do it again with a difference!
%
% Artificial data for two groups
NN = [10 15]; % number per group
NCases = sum(NN); NTreatments = 2;
MM = [5 7]; % means -- one sd's apart
SD = [2 2]; % standard deviation
Sample = [MM(1)+SD(1)*randn(NN(1),1); MM(2)+SD(2)*randn(NN(2),1)];
%
pause
clc
% Now calculate the various sums of squares:
%
% Individual -- 
II = Sample'*Sample
%
% Treatment --
T1 = Sample(1:NN(1)); T2 = Sample((NN(1)+1):NCases);
TT = (sum(T1)^2)/NN(1) + (sum(T2)^2)/NN(2)
%
% Combination --
CC = (sum(Sample)^2)/length(Sample)
%
pause
clc
% Now calculate the partition of variance...
% Total square deviations: II - CC with NCases-1 degrees of freedom
Total = II - CC   % df = 34
%
% Treatment: TT - CC with NTreatments-1 degrees of freedom
Treatment = TT - CC   % df = 1
%
% Residual: II - TT with  NCases - NTreatments degrees of freedom
Residual = II - TT   % df = NCases - NTreatments = 33
% Do they add up?
Treatment + Residual
pause
clc
% Now calculate the F ratio:
MSTR = Treatment/(NTreatments-1)  % df =1
MSE = Residual/(NCases-NTreatments) % df = 34
F = MSTR/MSE  % Should be larger, since this time there was a difference!
pause
clc
% Now compare against expectations for F(1,34)
F
finv(1-[.05 .01 .005 .001],1,34) % significance thresholds for .05, .01, etc.
pause
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OK, now try with a larger number of groups
%
% Artificial data for five groups
NN = [10 15 9 6 8]; % number per group
NCases = sum(NN); NTreatments = 5;
MM = [5 5 5 6 5]; % means -- a bit of a difference!
SD = [1 1 1 1 1]; % standard deviation
Sample = [MM(1)+SD(1)*randn(NN(1),1); MM(2)+SD(2)*randn(NN(2),1); ...
    MM(3)+SD(3)*randn(NN(3),1); MM(4)+SD(4)*randn(NN(4),1); ...
    MM(5)+SD(5)*randn(NN(5),1) ];
%
pause
clc
% Now calculate the various sums of squares:
%
% Individual -- 
II = Sample'*Sample
%
% Treatment --
TT = 0;
bounds = [1 cumsum(NN)+1];
for n = 1:NTreatments
    b1=bounds(n); b2=bounds(n+1)-1; TT = TT + (sum(Sample(b1:b2))^2)/NN(n);
end
TT
%
% Combination --
CC = (sum(Sample)^2)/NCases
%
pause
clc
% Now calculate the partition of variance...
% Total square deviations: II - CC with NCases-1 degrees of freedom
Total = II - CC   % df = 47
%
% Treatment: TT - CC with NTreatments-1 degrees of freedom
Treatment = TT - CC   % df = 4
%
% Residual: II - TT with  NCases - NTreatments degrees of freedom
Residual = II - TT   % df = NCases - NTreatments = 43
% Do they add up?
Treatment + Residual
pause
clc
% Now calculate the F ratio:
MSTR = Treatment/(NTreatments-1)  % df = 4
MSE = Residual/(NCases-NTreatments) % df = 43
F = MSTR/MSE  % Should be large enough?
pause
clc
% Now compare against expectations for F(4,43)
F
% significance thresholds for .05, .01, etc.:
finv(1-[.05 .01 .005 .001],NTreatments-1,NCases-NTreatments)
pause





%
