/*data from John and Williams, 1995, p.146)*/
data a;
input
rep   block     gen     y;
datalines;
 1       1       11    4.1172
 1       1        4    4.4461
 1       1        5    5.8757
 1       1       22    4.5784
 1       2       21    4.6540
 1       2       10    4.1736
 1       2       20    4.0141
 1       2        2    4.3350
 1       3       23    4.2323
 1       3       14    4.7572
 1       3       16    4.4906
 1       3       18    3.9737
 1       4       13    4.2530
 1       4        3    3.3420
 1       4       19    4.7269
 1       4        8    4.9989
 1       5       17    4.7876
 1       5       15    5.0902
 1       5        7    4.1505
 1       5        1    5.1202
 1       6        6    4.7085
 1       6       12    5.2560
 1       6       24    4.9577
 1       6        9    3.3986
 2       1        8    3.9926
 2       1       20    3.6056
 2       1       14    4.5294
 2       1        4    4.3599
 2       2       24    3.9039
 2       2       15    4.9114
 2       2        3    3.7999
 2       2       23    4.3042
 2       3       12    5.3127
 2       3       11    5.1163
 2       3       21    5.3802
 2       3       17    5.0744
 2       4        5    5.1202
 2       4        9    4.2955
 2       4       10    4.9057
 2       4        1    5.7161
 2       5        2    5.1566
 2       5       18    5.0988
 2       5       13    5.4840
 2       5       22    5.0969
 2       6       19    5.3148
 2       6        7    4.6297
 2       6        6    5.1751
 2       6       16    5.3024
 3       1       11    3.9205
 3       1        1    4.6512
 3       1       14    4.3887
 3       1       19    4.5552
 3       2        2    4.0510
 3       2       15    4.6783
 3       2        9    3.1407
 3       2        8    3.9821
 3       3       17    4.3234
 3       3       18    4.2486
 3       3        4    4.3960
 3       3        6    4.2474
 3       4       12    4.1746
 3       4       13    4.7512
 3       4       10    4.0875
 3       4       23    3.8721
 3       5       21    4.4130
 3       5       22    4.2397
 3       5       16    4.3852
 3       5       24    3.5655
 3       6        3    2.8873
 3       6        5    4.1972
 3       6       20    3.7349
 3       6        7    3.6096
;

/*ad hoc h2*/
ods output diffs=diffs;
proc mixed data=a lognote;
class rep block gen;
model y=gen rep;
random rep*block;
lsmeans gen/pdiff;
run;

data diffs;
set diffs;
half_var_e=0.5*stderr**2;
run;

proc means data=diffs mean;
var half_var_e;
output out=half mean=;
run;

proc mixed data=a lognote mmeqsol;
class rep block gen;
model y=rep;
random gen rep*block/g;
ods output covparms=cp mmeqsol=mmeqsol g=g;
run;

data cp;
set cp;
if covparm='gen';
run;

data cp;
merge cp half;
h2=estimate/(estimate+half_var_e);

proc print data=cp;
run;

/*simulated*/

proc iml;

/*read from G matrix*/
use g;
read all var _NUM_ into G;
G=G[,5:ncol(G)];
from_g=1; to_g=24;             /*manually identify rows pertaining to genotypic effects*/
D=G[from_g:to_g,from_g:to_g];
F=G[from_g:to_g,];
*print g d f;

/*read coefficients from MME*/
use mmeqsol;
read all var _NUM_ into C22;
c22=c22[,5:ncol(c22)-1];
from_C22=5; to_C22=46;             /*manually identify rows pertaining to genotypic effects*/
c22=c22[from_c22:to_c22,from_c22:to_c22];

M=G-C22;
inv_G=inv(G);
Q=F*inv_G*M*inv_G*t(F);
omega=(D||q)//(q||q);
root_omega=root(omega);
gamma=t(root_omega);  
n=nrow(gamma);
z=j(n,1,0);
r2=0;
sim_max=10000;
a=0; b=0; c=0;
do sim=1 to sim_max;
  do i=1 to n;
    z[i]=normal(1);
  end;
  w=gamma*z;
  g=w[1:n/2];       g=g-g[:];
  g_hat=w[n/2+1:n]; g_hat=g_hat-g_hat[:];
  r2=r2+ ( t(g_hat)*g )**2/( t(g)*g*t(g_hat)*g_hat );
  a=a+t(g_hat)*g;
  b=b+t(g)*g;
  c=c+t(g_hat)*g_hat;
end;
r2=r2/sim_max; 
r2_=a*a/(b*c);
print r2 r2_; 

quit;
run;