clear all
close all
% The data showing the score of normalized damage events are shown
% in the table below. There were n1 = 13 donors classified as young (?45
% years old) and n2 = 17 classified as old (>45 years old). To calculate the
% microdamage progression score, the counts of damage events (extensions,
% surface originations, widenings, and combinations) are normalized to the
% bone area and summed up.
% Test the equality of mean scores for the two groups (young/old).
microdamage = [... 
   1     0.1582         0    0.3164    0.3164; ...
   1     0.5902    0.1180    0.2361         0; ...
   1     0.1196    0.2392    0.4785    0.1196;...
   1     0.2527    0.1264    0.6318         0;...
   1     0.2378         0    0.2378    0.2378;...
   1     0.1276         0         0    0.1276;...
   1          0         0    0.4058         0;...
   1          0         0    0.1349         0;...
   1          0         0    0.2106    0.1053;...
   1          0         0    0.1788         0;...
   1     0.5619         0    0.2809    0.4214;...
   1     0.1566         0    0.3132    0.9395;...
   1     0.6327         0    0.3163    0.2109;...
   2     0.2746         0    0.8239    0.2746;...
   2     0.3609         0    0.1203    0.1203;...
   2          0         0    0.6856    0.3428;...
   2     0.4740    0.1580    0.4740    0.1580;...
   2     0.1690         0    0.3379    0.6758;...
   2     0.4638    0.1546    0.7730    0.4638;...
   2     0.6332         0    0.7914    0.4749;...
   2     0.1622         0    0.3244         0;...
   2     0.2033         0    0.2033    0.4067;...
   2     0.4101         0    0.2734    0.1367;...
   2     0.2948         0    0.4421    0.5895;...
   2     0.6621    0.1655    0.1655    0.3311;...
   2     0.5029         0    1.0059    0.5029;...
   2     0.2053         0    0.5131    0.3079;...
   2     0.6780         0    0.3390    0.1130;...
   2          0    0.1208    0.2415    0.2415;...
   2     0.1738         0    0.2607    0.4345];
microdamage = round(1000 * microdamage)/1000;
agegroup    = microdamage(:,1);
extension   = microdamage(:,2);
surforigin  = microdamage(:,3);
widening    = microdamage(:,4);
combination = microdamage(:,5);
score = extension + surforigin +  widening + combination;

x1 = score(agegroup==1);
x2 = score(agegroup==2);
m1 = mean(x1); n1 = length(x1); s1=std(x1);
m2 = mean(x2); n2 = length(x2); s2=std(x2);

m1  %0.7341
m2  %1.1541
n1   %13
n2   %17
s1  %0.4351
s2  %0.4541
Fstat = s1^2/s2^2 
     % Fstat =  0.9181
%For two sided alternative the p-value is 
pval = 2 * min(fcdf(Fstat, n1-1, n2-1), 1-fcdf(Fstat, n1-1, n2-1)) %0.8972

sp = sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2)/(n1 + n2 - 2)) %0.4460

t1 = (m1 - m2)/(sp * sqrt(1/n1 + 1/n2)) %-2.5560
t2 = (m1 - m2)/sqrt(s1^2/n1 + s2^2/n2)  %-2.5711

df1=n1+n2-2 %28
df2 = (s1^2/n1 + s2^2/n2)^2 /((s1^2/n1)^2/(n1-1) + (s2^2/n2)^2/(n2-1)) %26.5179

pval1 = tcdf(t1, df1) %0.0082
pval2 = tcdf(t2, df2) %0.0080

 
     
