data norep;
  input a b y;
  y=y/10;
  cards;
1 1 21
1 2 68
1 3 84
1 4 14
1 5 146
1 6 79
1 7 48
2 1 17
2 2 81
2 3 84
2 4 14
2 5 120
2 6 37
2 7 45
3 1 144
3 2 148
3 3 270
3 4 309
3 5 365
3 6 364
3 7 314
4 1 574
4 2 624
4 3 374
4 4 633
4 5 655
4 6 656
4 7 598
5 1 662
5 2 817
5 3 533
5 4 807
5 5 797
5 6 808
5 7 824
6 1 752
6 2 940
6 3 743
6 4 879
6 5 933
6 6 878
6 7 805
7 1 41
7 2 102
7 3 107
7 4 55
7 5 181
7 6 114
7 7 61
;

proc means noprint;
  var y;
  output out=overall mean=mean;
run;

proc sort data=norep out=norep;
  by a;
run;

proc means noprint data=norep;
  var y;
  by a;
  output out=ameans mean=ybari;
run;

proc sort data=norep out=norep;
  by b;
run;

proc means noprint data=norep;
  var y;
  by b;
  output out=bmeans mean=ybarj;
run;

data norep;
  merge norep overall;
  drop _TYPE_ _FREQ_ mean;
  if _N_=1 then ybar=mean;
    else ybar=ybar; 
  retain ybar;

data norep;
  merge norep bmeans;
  by b;

proc sort data=norep out=norep;
  by a;
run;

data norep;
  merge norep ameans;
  by a;
  drop _TYPE_ _FREQ_;
  ybari_ybar=ybari-ybar;
  ybarj_ybar=ybarj-ybar;
  Pij=y*ybari_ybar*ybarj_ybar;
  isq=ybari_ybar**2;
  jsq=ybarj_ybar**2;

proc means data=norep noprint;
  var Pij;
  output out=sumPij sum=P;
run;

data justb; set norep;
  if a^=1 then delete;

proc means data=justb noprint;
  var jsq;
  output out=sumjsq sum=sumjsq;
run; 

data justa; set norep;
  if b^=1 then delete;

proc means data=justa noprint;
  var isq;
  output out=sumisq sum=sumisq;
run; 

data nonadd;
  merge sumPij sumisq sumjsq;
  SSnonadd=P**2/(sumisq*sumjsq);
  keep SSnonadd P sumisq sumjsq;


proc glm data=norep outstat=glmout;
  class a b;
  model y=a b;
run;

data test;
  merge glmout nonadd;
  SSE=SS-SSnonadd;
  DFE=DF-1;
  MSE=SSE/DFE;
  F=SSnonadd/MSE;
  pvalue=1-probf(F,1,DFE);
  if _N_>1 then delete;
  keep SSnonadd SSE MSE DFE F pvalue;

proc print;
run;