/* ch 7 one obs per cell -- orthogonal contrasts */
/* normal prob plot */
options ls=72;
data beads;
  INPUT sizelet $ sizecat finglet $ finglnth fingcat color $ colcat
   time a1 a2 b1 b2 c1;
  htime=log10(tan((time*3.14159/180)/2)-.3);
  code=sizecat*100+fingcat*10+colcat;
  letcode=trim(sizelet)||trim(finglet)||trim(color);
   a1a1=a1*a1;
   a2a2=a2*a2;
   b1b1=b1*b1;
   b2b2=b2*b2;
   c1c1=c1*c1;
   a1b1=a1*b1*a1*b1;
   a1b2=a1*b2*a1*b2;
   a2b1=a2*b1*a2*b1;
   a2b2=a2*b2*a2*b2;
   a1c1=a1*c1*a1*c1;
   a2c1=a2*c1*a2*c1;
   b1c1=b1*c1*b1*c1;
   b2c1=b2*c1*b2*c1;
   a1b1c1=a1*b1*c1*a1*b1*c1;
   a1b2c1=a1*b2*c1*a1*b2*c1;
   a2b1c1=a2*b1*c1*a2*b1*c1;
   a2b2c1=a2*b2*c1*a2*b2*c1;
  ha1a1=a1*htime;
  ha2a2=a2*htime;
  hb1b1=b1*htime;
  hb2b2=b2*htime;
  hc1c1=c1*htime;
  ha1b1=a1*b1*htime;
  ha1b2=a1*b2*htime;
  ha2b1=a2*b1*htime;
  ha2b2=a2*b2*htime;
  ha1c1=a1*c1*htime;
  ha2c1=a2*c1*htime;
  hb1c1=b1*c1*htime;
  hb2c1=b2*c1*htime;
  ha1b1c1=a1*b1*c1*htime;
  ha1b2c1=a1*b2*c1*htime;
  ha2b1c1=a2*b1*c1*htime;
  ha2b2c1=a2*b2*c1*htime;
  LINES;
   s 1 t 64 1 o 1 100 -1 1 -1 1    1
   s 1 a 72 2 o 1 140  -1 1 0 -2   1
   s 1 b 76 3 o 1 93  -1 1 1 1     1
   m 2 t 70 1 o 1 85 0 -2 -1 1     1
   m 2 a 71 2 o 1 45  0 -2 0 -2    1
   m 2 b 79 3 o 1 105  0 -2 1 1    1
   l 3 t 70 1 o 1 40  1 1 -1 1     1
   l 3 a 75 2 o 1 37  1 1 0 -2     1
   l 3 b 77 3 o 1 35  1 1 1 1      1
   s 1 t 63 1 b 2 119 -1 1  -1  1  -1
   s 1 a 72 2 b 2 127 -1 1   0  -2  -1
   s 1 b 79 3 b 2  70 -1 1   1   1  -1
   m 2 t 69 1 b 2  49  0 -2  -1  1  -1
   m 2 a 71 2 b 2 133  0 -2   0  -2 -1
   m 2 b 80 3 b 2  86  0 -2  1   1  -1
   l 3 t 65 1 b 2  44  1 1  -1   1  -1
   l 3 a 75 2 b 2  35  1 1   0   -2 -1
   l 3 b 77 3 b 2  53  1 1   1   1  -1
  ;

proc means noprint sum;
 var a1a1 a1b1 a1b2 a1c1 a2a2 a2b1 a2b2 a2c1 b1b1 b1c1 b2b2 c1c1
     b2c1 a1b1c1 a2b1c1 a1b2c1 a2b2c1
    ha1a1 ha1b1 ha1b2 ha1c1 ha2a2 ha2b1 ha2b2 ha2c1 hb1b1 hb1c1
    hb2b2 hc1c1  hb2c1 ha1b1c1 ha2b1c1 ha1b2c1 ha2b2c1;
output out=intmdest
 sum= a1a1 a1b1 a1b2 a1c1 a2a2 a2b1 a2b2 a2c1 b1b1 b1c1 b2b2 c1c1
     b2c1 a1b1c1 a2b1c1 a1b2c1 a2b2c1
    ha1a1 ha1b1 ha1b2 ha1c1 ha2a2 ha2b1 ha2b2 ha2c1 hb1b1 hb1c1
    hb2b2 hc1c1  hb2c1 ha1b1c1 ha2b1c1 ha1b2c1 ha2b2c1;
data estimats;
set intmdest;
vla1=ha1a1/sqrt(a1a1);
vla2=ha2a2/sqrt(a2a2);
vlb1=hb1b1/sqrt(b1b1);
vlb2=hb2b2/sqrt(b2b2);
vlc1=hc1c1/sqrt(c1c1);
vla1b1=ha1b1/sqrt(a1b1);
vla1b2=ha1b2/sqrt(a1b2);
vla1c1=ha1c1/sqrt(a1c1);
vla2b1=ha2b1/sqrt(a2b1);
vla2b2=ha2b2/sqrt(a2b2);
vla2c1=ha2c1/sqrt(a2c1);
vlb1c1=hb1c1/sqrt(b1c1);
vlb2c1=hb2c1/sqrt(b2c1);
vla1b1c1=ha1b1c1/sqrt(a1b1c1);
vla2b1c1=ha2b1c1/sqrt(a2b1c1);
vla1b2c1=ha1b2c1/sqrt(a1b2c1);
vla2b2c1=ha2b2c1/sqrt(a2b2c1);
data newest; set estimats;
keep
  vla1 vla1b1 vla1b2 vla1c1 vla2
  vla2b1 vla2b2 vla2c1 vlb1 vlb1c1 vlb2 vlc1
     vlb2c1 vla1b1c1 vla2b1c1 vla1b2c1 vla2b2c1;
proc transpose prefix=est out=ests;
data ests; set ests;
proc print;
proc rank normal=blom out=plt;
var est1;
 ranks nscore;
proc print;
 proc plot;
  plot est1*nscore/vpos=19 hpos=50;
run; quit;