/* This example simulates data from a one-way 
breeding experiment like the on around p. 175 of
the lecture packet.  The model involves a random effect
for sire, T_i.  The magnitude of the variance component
may be set as the &sigmaT macro variable.  Fixed effects
and random effects models are fit and confidence intervals
for the population mean mu=80 are obtained.  Note how many
times the confidence interval from the fixed effects
simulations "miss" the target of mu=80*/

options ls=75 nodate;

%let sigmaT=10;
%let sigma=20;

*%let seed=0;
%let seed=1234;

/* use seed=0 to generate new data each time */
/* use nonzero seed to regenerate the same random data */

/* The next 40 or so lines of code constitute a SAS MACRO
that is used to generate the data.  You can regenerate
new data with a single line to call this macro.  The
macro code ends at the %MEND statement. */

%macro datagen(t=,nreps=,nsims=);

data one; 
   do simulation=1 to &nsims;
   do trt=1 to &t;
      Ti=normal(&seed)*&sigmaT;
      do rep=1 to &nreps;
          Eij=normal(&seed)*σ
          y=80+Ti+Eij;
          output;
      end;
   end;
   end;
run;

ods trace on;
proc mixed data=one;
   by simulation;
   class trt;
   model y=/cl;
   random trt;
   estimate "mean" intercept 1/cl;
   ods output estimates=ests covparms=covparms;
run;


proc mixed data=one;
   by simulation;
   class trt;
   model y=trt/cl;
   *random trt;
   estimate "mean" intercept 1/cl;
   ods output estimates=estsfixed;
run;

%mend datagen;


ods listing close;   /* Stops sending output to the output window */

%datagen(t=5,nreps=10,nsims=200);

data ests;
   set ests(rename=(estimate=ybar));
   tstat=abs(ybar-80)/Stderr;
   pvalue=2*(1-probt(tstat,DF));
   sim=_n_;
run;
data estsfixed;
   set estsfixed(rename=(estimate=ybar));
   tstat=abs(ybar-80)/Stderr;
   pvalue=2*(1-probt(tstat,DF));
   sim=_n_;
run;

ods listing ;   /* Restarts sending output to the output window */

*proc print data=covparms (obs=20);
proc print data=covparms; 
   title "output from PROC MIXED macro - estimated variance components";
run;

proc sort data=covparms;
   by CovParm;
run;

proc gchart;
   title "estimated variance components, remember, sigma=&sigma,sigmaT=&sigmaT";
   by covparm;
   vbar estimate;
run;

proc print data=ests;
   title "output from PROC MIXED macro - random effects";
   var Label Estimate StdErr DF lower upper tstat pvalue;
run;
proc print data=estsfixed;
   title "output from PROC MIXED  macro - fixed effects";
   var Label Estimate StdErr DF lower upper tstat pvalue;
run;

symbol value=dot;
proc gplot data=ests;
   title "confidence intervals for mu using random effects model";
   plot upper*sim=1 lower*sim=2/overlay vref=80;
run;

proc gplot data=estsfixed;
   title "confidence intervals for mu using fixed effects model";
   plot upper*sim lower*sim/overlay vref=80;
run;