/* sandm11.sas                     */
 /*  analysis of sandm10.r output   */
 /*                                 */
 /* three methods, two param ests,  */
 /*  conv codes, sigma-sq-a, pattern */
  options ls=80 ;
  data a ;
    infile 'estsall.dat' ;          * read data ;
    keep sqere sqera method siga pattern batch ;
    input a1 a2 m1 m2 r1 r2 cm cr siga pattern ;
    rep = mod(_n_-1,100) ;          * only 100 this time ;
    batch = int(rep/20) ;           * 5 batches of 20 ;
 /* make 3 obs from each -- bs & i  */
    method = 'anova' ;    sqere = (a1-1)**2 ;
    sqera = (a2-siga)**2 ;    output ;
    method= 'ml' ;        sqere = (m1-1)**2 ;
    sqera = (m2-siga)**2 ;    output ;
    method = 'reml' ;     sqere = (r1-1)**2 ;
    sqera = (r2-siga)**2 ;    output ;
 run ;
 /* get batch means                  */
 /*                      nway noprint*/
  proc means data=a mean nway noprint ;
    class batch method pattern siga ;
    var sqere sqera ;
    output out=b mean= ;             * output ds with sqere sqera ;
 run ;
 /* print a few                      */
 /***************not this time*****
  proc print data=b (obs=10) ;
    title 'batch means' ;
 run ;
 /* analyze to compare mse         */
 /*  means & stdev of estimators   */
  proc glm data=b ;
    class method siga pattern batch ;
    model sqera = method siga method*siga pattern pattern*method 
                   pattern*siga pattern*batch*siga ;
    random pattern*batch*siga ;
    lsmeans method siga pattern ;
    title 'Analysis of MSE' ;
 run ;