options ls=75 nodate;

data latinsq;
   /* ARRAY statements enable several variables to be processed as a group */
   array sqware{3} (3,0,-3);  
   array slight1{3} (1,0,-1);
   array slight2{3} (1,0,-1);
   array slight3{3} (1,0,-1);
   array iheight{3} (2,0,-2); /* this has been fixed.  */
   array treatment{3} (-1,0,1);
   keep square row row9 col trt growth; 
   input square row col trt; 
   row9=3*(square-1)+row;  /* here we define a new row variable  taking value 1 through 9 */
   sigma=1;  /* try various values of sigma to generate the data */
   if square=1 then do;
      growth = 10 + sqware{square} + slight1{row} + iheight{col}+ treatment{trt} + sigma*rannor(123);
   end;
   /* the argument of the RANNOR function is a seed, for random number generation */
   else if square=2 then do;
      growth = 10 + sqware{square} + slight2{row} + iheight{col}+ treatment{trt} + sigma*rannor(123);
   end;
   else do;
      growth = 10 + sqware{square} + slight3{row} + iheight{col}+ treatment{trt} + sigma*rannor(123);
   end;
   cards;
   1  1  1  3  
   1  1  2  1 
   1  1  3  2
   1  2  1  1  
   1  2  2  2 
   1  2  3  3
   1  3  1  2  
   1  3  2  3 
   1  3  3  1
   2  1  1  2  
   2  1  2  3 
   2  1  3  1
   2  2  1  1  
   2  2  2  2 
   2  2  3  3
   2  3  1  3  
   2  3  2  1
   2  3  3  2
   3  1  1  2  
   3  1  2  1
   3  1  3  3
   3  2  1  3  
   3  2  2  2 
   3  2  3  1
   3  3  1  1  
   3  3  2  3
   3  3  3  2
 ;

data latinsq;
   set latinsq;
   fert=trt;
   iheight=col;
   slight=row;
   growth=round(growth,0.1);
run;

data one;  /* this is a subset of the data corresponding to the first three rows */
   set latinsq;
   if square=1;
run;
   
proc print;
   var square row col trt growth;

proc glm data=one;
   class square row col trt;
   model growth = row col trt;
   lsmeans row col trt;
run;
proc glm data=latinsq;
   class square row row9 col trt;
   *model growth = square row(square) col trt;
   model growth = row9 col trt;
   *lsmeans square row(square) col trt;
   lsmeans row9 col trt;
run;

options ps=25;

proc plot;
   plot growth*trt;
   *plot growth*trt=iheight;
   *plot growth*trt=slight;
   *plot growth*trt=square;
   plot growth*trt=col;
   plot growth*trt=row;
   plot growth*trt=square;
run;