load schizo.dat;
grp=schizo(:,1);
Y=sqrt(schizo(:,2:5));

%
% Construct Plot
%

[N,d]=size(Y);
X=[ones(N,1) dummyvar(grp)];
B_tilde=pinv(X'*X)*X'*Y
L=[1 1 0;
   1 0 1];
ybar=(L*B_tilde)'
clear L
minutes=[0 15 30 45]';
plot(minutes,ybar(:,1),'-',minutes,ybar(:,2),'--')
set(gca,'FontSize',16)
xlabel('Time after Insulin')
ylabel('Root fatty Acid')
set(gca,'LineWidth',2)
legend('Normal','Schizophrenic','Location','SouthWest')
hold on
plot(minutes,ybar(:,1),'*',minutes,ybar(:,2),'*','MarkerSize',10)
hold off
print -deps profile2m.eps

%
% Construct Coefficients for Orthogonal Polynomial Trend Components
%

lin=minutes;
Poly=[lin lin.^2 lin.^3];
Poly_Orth=[];
K=ones(d,1);
for i=1:3
  m=(eye(d)-ppo(K))*Poly(:,i);
  Poly_Orth=[Poly_Orth m/sqrt(m'*m)];
  K=[K Poly(:,i)];
end

M{1} = ones(d,1)/4;
M{2} = null(ones(d,1)');
L{1} = [2 1 1]/2;
L{2} = [0 1 -1];
MM{1}=M{1};
MM{2}=[Poly_Orth [1 1 -1 -1]'/2];

%
% Test Ho: L{j}B M{i} = 0
% Construct Simultaneous CIs for L{j}*B*MM{i}, Col(MM{i}) in Col(M{i})
%

r=rank(X);
S=Y'*(eye(N)-ppo(X))*Y/(N-r);
tit=[' i  j      T^2       F    df1  df2  p-value  crit'];
tit_ci=['    Psi    SE    Lower  Upper'];
form=[' %d  %d  %8.3f  %8.3f  %d   %d   %5.4f %6.3f '];
form_ci=[' %6.3f %6.3f  %6.3f %6.3f '];
for i=1:2
  for j=1:2
    Psi=L{j}*B_tilde*M{i};
    V=M{i}'*S*M{i};
    k_N=inv(L{j}*pinv(X'*X)*L{j}');
    T2=k_N*Psi*inv(V)*Psi';
    q=size(M{i},2);
    F=T2*(N-r-q+1)/(q*(N-r));
    pval=1-fcdf(F,q,N-r-q+1);
    crit=finv(.95,q,N-r-q+1);
    crit=crit*q*(N-r)/(N-r-q+1);
    M_opt=M{i}*inv(V)*Psi';
    ii=find(abs(M_opt)==max(max(abs(M_opt))));
    M_opt=sign(M_opt(ii(1)))*M_opt/abs(M_opt(ii(1)));
    disp(['L, M and MM'])
    L_j=L{j}
    M_i=M{i}
    MM_i=[MM{i} M_opt]
    disp(tit)
    [Str, errmsg]=sprintf(form, [i j T2 F q N-r-q+1 pval crit]);
    disp(Str)
    Psi=L{j}*B_tilde*MM_i;
    V=MM_i'*S*MM_i;
    SE=sqrt(diag(V)/k_N);
    CI=[Psi' SE Psi'-SE*sqrt(crit) Psi'+SE*sqrt(crit)];
    disp(tit_ci)
    for jj=1:size(CI,1)
      [Str, errmsg]=sprintf(form_ci,CI(jj,:));
      disp(Str)
    end
  end
end