/* sandm9.sas                     */
 /*  analysis of sandm7.r output   */
 /*                                */
 /* three methods, two param ests, */
 /*  conv codes, sigma-sq-a        */
  options ls=80 ;
  data a ;
    infile 'estsP4.dat' ;          * read data ;
    keep sqere sqera method siga batch ;
    input a1 a2 m1 m2 r1 r2 cm cr siga ;
    rep = mod(_n_-1,1000) ;
    batch = int(rep/20) ;           * 50 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 siga ;
    var sqere sqera ;
    output out=b mean= ;             * output ds with sqere sqera ;
 run ;
 /* print a few                      */
  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 batch ;
    model sqera = method siga method*siga batch(siga) ;
    random batch(siga) / test ;
    lsmeans method siga ;
    title 'Analysis of MSE' ;
 run ;