/* A N A L Y S E  R E A L  D A T A S E T */

/*Estimate variance components*/

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 univariate data=diffs;
var stderr;
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;

/*end of simulation for empirical datset*/

/*predict response to selection by simulation*/

proc iml;

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;

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];
*print c22;

M=G-C22;
inv_G=inv(G);
Q=F*inv_G*M*inv_G*t(F);
omega=(D||q)//(q||q);  
call svd(u,lambda,v, omega); 
gamma=u*diag(sqrt(lambda)); 
n=nrow(gamma);
z=j(n,1,0);
r2=0;
sim_max=100000;
a=0; b=0; c=0;
sel_mean=j(n/2,1,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_hat=w[n/2+1:n];
  r_g_hat=rank(g_hat);
  v=g_hat||r_g_hat||g; 
  call sort(v, {1}, {1}); 
  do j=1 to n/2;
    sel_mean[j]=sel_mean[j]+sum(v[1:j,3])/j;
  end;
  g=g-g[:];
  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_; 
create r2 var {r2 r2_};
append;

sel_mean=sel_mean/sim_max;

use cp;
read all var _num_ into var_g; print var_g;
h2=var_g[5];
var_g=var_g[1];
i=n/2;
one=j(i,1,1);
c22=c22[1:24, 1:24];
a=trace(c22)-t(one)*c22*one/i; a=2*a/(i-1);
h2_c=1-a/2/var_g;
print h2_c;

create sel_mean var {sel_mean};
append;

create h2 var {h2 h2_c var_g};
append;

print sel_mean;

quit;

proc print data=sel_mean;
run;

/* ad hoc approach*/
proc iml;
use h2;
read all var _num_ into y;
print y;
h2=y[1];
h2_c=y[2];
var_g=y[3];
var_p=var_g/h2_c;
var_e=var_p-var_g;
sim_max=100000;
n=24;
g_hat=j(n,1,0);
g=g_hat;
sel_mean_naive=g;
do sim=1 to sim_max;
  do i=1 to n;
    g[i]=sqrt(var_g)*normal(1);
	g_hat[i]=g[i]+sqrt(var_e)*normal(1);
  end;
  r_g_hat=rank(g_hat);
  v=g_hat||r_g_hat||g; 
  call sort(v, {1}, {1}); 
  do i=1 to n;
    sel_mean_naive[i]=sel_mean_naive[i]+sum(v[1:i,3])/i;
  end;
end;

sel_mean_naive=sel_mean_naive/sim_max;

print sel_mean_naive;

create sel_mean_naive var {sel_mean_naive};
append;

quit;
run;

proc print data=sel_mean_naive;
run;

/*end of: predict response to selection by simulation*/

/*B O O T S T R A P*/

/* this bootstrap code essentially repeats the code for the analysis of the dataset in 
   bootstrap a do-loop*/

%macro boot(b=100);
  
/**************turn printing off ************************/

proc printto log='d:\hpp\log.tmp' new; run;
options nospool;
ods results=off;
ods listing close;

/*first get V*/
proc mixed data=a;
class rep block gen;
model y=rep;
random gen rep*block/vc;
ods output cholv=cholv;
run;

data r2_sim2;
run;

/*now bootstrap*/
%do i=1 %to &b;

  data b;
    array col col1-col72;
	array z z1-z72;
	do i=1 to 72;
	  z[i]=normal(-1);
	end;
	do j=1 to 72;                         /* y is a vector of length 72 in this example*/
      set cholv nobs=size point=j;
	  y_sim=0;
	  do i=1 to 72;
	    y_sim=y_sim+col[i]*z[i];
	  end;
	  output;
    end;
	stop;
  run;

  data b;
  set b;
  keep y_sim;

  data b;
  merge a b;

  proc mixed data=b lognote mmeqsol;
  class rep block gen;
  model y_sim=rep;
  random gen rep*block/g;
  ods output mmeqsol=mmeqsol g=g;
  run;

  /*simulated*/

  proc iml;

  use g;
  read all var _NUM_ into G;
  G=G[,5:ncol(G)];
  from_g=1; to_g=24;             /*check manually which rows of G contain the genetic effects*/
  D=G[from_g:to_g,from_g:to_g];
  F=G[from_g:to_g,];

  use mmeqsol;
  read all var _NUM_ into C22;
  c22=c22[,5:ncol(c22)-1];
  from_C22=5; to_C22=46;             /*check manually which rows of G contain the genetic 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=100000;
  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_sim=r2/sim_max; 
  r2_sim_=a*a/(b*c);
  *print r2 r2_; 

  create r2_sim var {r2_sim r2_sim_};
  append;
  quit;

  data r2_sim2;
  set r2_sim2 r2_sim;
  run;

  %put &i;

%end;

proc means data=r2_sim2 noprint;
var r2_sim r2_sim_;
output out=result mean=mean1 mean2 var=var1 var2 std=std1 std2;
run;

/************turn printing back on***********************************/

proc printto log=log;
ods listing;

proc print data=result;
run;

%mend;

%boot(b=1000);

data result2;
merge result r2;
bias1=mean1-r2;
bias2=mean2-r2_;
mse1=bias1**2+var1;
mse2=bias2**2+var2;
run;

proc print data=result2;
run;