G1=[...
   1  1  4.3  3.3  3.0  2.6  2.2  2.5  3.4  4.4
   2  1  3.7  2.6  2.6  1.9  2.9  3.2  3.1  3.9
   3  1  4.0  4.1  3.1  2.3  2.9  3.1  3.9  4.0
   4  1  3.6  3.0  2.2  2.8  2.9  3.9  3.8  4.0
   5  1  4.1  3.8  2.1  3.0  3.6  3.4  3.6  3.7
   6  1  3.8  2.2  2.0  2.6  3.8  3.6  3.0  3.5
   7  1  3.8  3.0  2.4  2.5  3.1  3.4  3.5  3.7
   8  1  4.4  3.9  2.8  2.1  3.6  3.8  4.0  3.9
   9  1  5.0  4.0  3.4  3.4  3.3  3.6  4.0  4.3
  10  1  3.7  3.1  2.9  2.2  1.5  2.3  2.7  2.8
  11  1  3.7  2.6  2.6  2.3  2.9  2.2  3.1  3.9
  12  1  4.4  3.7  3.1  3.2  3.7  4.3  3.9  4.8
  13  1  4.7  3.1  3.2  3.3  3.2  4.2  3.7  4.3];
G2=[...
   1  2  4.3  3.3  3.0  2.6  2.2  2.5  2.4  3.4
   2  2  5.0  4.9  4.1  3.7  3.7  4.1  4.7  4.9
   3  2  4.6  4.4  3.9  3.9  3.7  4.2  4.8  5.0
   4  2  4.3  3.9  3.1  3.1  3.1  3.1  3.6  4.0
   5  2  3.1  3.1  3.3  2.6  2.6  1.9  2.3  2.7
   6  2  4.8  5.0  2.9  2.8  2.2  3.1  3.5  3.6
   7  2  3.7  3.1  3.3  2.8  2.9  3.6  4.3  4.4
   8  2  5.4  4.7  3.9  4.1  2.8  3.7  3.5  3.7
   9  2  3.0  2.5  2.3  2.2  2.1  2.6  3.2  3.5
  10  2  4.9  5.0  4.1  3.7  3.7  4.1  4.7  4.9
  11  2  4.8  4.3  4.7  4.6  4.7  3.7  3.6  3.9
  12  2  4.4  4.2  4.2  3.4  3.5  3.4  3.9  4.0
  13  2  4.9  4.3  4.0  4.0  3.3  4.1  4.2  4.3
  14  2  5.1  4.1  4.6  4.1  3.4  4.2  4.4  4.9
  15  2  4.8  4.6  4.6  4.4  4.1  4.0  3.8  3.8
  16  2  4.2  3.5  3.8  3.6  3.3  3.1  3.5  3.9
  17  2  6.6  6.1  5.2  4.1  4.3  3.8  4.2  4.8
  18  2  3.6  3.4  3.1  2.8  2.1  2.4  2.5  3.5
  19  2  4.5  4.0  3.7  3.3  2.4  2.3  3.1  3.3
  20  2  4.6  4.4  3.8  3.8  3.8  3.6  3.8  3.8];
G1=G1(:,3:10);
G2=G2(:,3:10);
[n1,p]=size(G1);
[n2,p]=size(G2);
I1=eye(n1);
j1=ones(n1,1);
gbar1=(1/n1)*G1'*j1;
J1=ones(size(I1));
S1=(1/(n1-1))*G1'*(I1-(1/n1)*J1)*G1;
I2=eye(n2);
j2=ones(n2,1);
gbar2=(1/n2)*G2'*j2;
J2=ones(size(I2));
S2=(1/(n2-1))*G2'*(I2-(1/n2)*J2)*G2;
gbar=(n1*gbar1+n2*gbar2)/(n1+n2);
W=(n1-1)*S1+(n2-1)*S2;
B=n1*(gbar1-gbar)*(gbar1-gbar)'+n2*(gbar2-gbar)*(gbar2-gbar)';
SX=[G1;G2];
IS=eye(size(SX,1));
JS=ones(size(IS));
S2=(1/(size(SX,1)-1))*SX'*(IS-(1/size(SX,1))*JS)*SX;
%
tMAdj=([0 .5 1 1.5 2 3 4 5]'-(ones(8,1))*mean([0 .5 1 1.5 2 3 4 5]))
A=[ones(8,1) tMAdj tMAdj.^2 tMAdj.^3 tMAdj.^4 tMAdj.^5 tMAdj.^6]
betahat=inv(A'*inv(S2)*A)*A'*inv(S2)*gbar
T21=33*(gbar-A*betahat)'*inv(S2)*(gbar-A*betahat)
T22=33*(gbar'*inv(S2)*gbar-gbar'*inv(S2)*A*betahat)
%
hours=[0 .5 1 1.5 2 3 4 5];
axis([0 6 2.5 5])
hold on
plot(hours,gbar1,':')
text('Position',[hours(8:8),gbar1(8:8)],'String','Control')
plot(hours,gbar2,'-.')
text('Position',[hours(8:8),gbar2(8:8)],'String','Obese')
plot(hours,gbar,'-')
text('Position',[hours(8:8),gbar(8:8)],'String','All groups')
title ('Profiles and mean profile for Two Patient groups')
ylabel('Plasma Inorganic Phosphate')
xlabel('Hours')
hold off


