/* This example simulates data from a one-way experiment
with a random treatment effect.  The magnitude of the
variance component may be set in the %let sigmaT=; statement
below.  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 seed=0;
%let sigmaT=20;
%let sigma=20;

%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;

proc print data=estsfixed (obs=20);
run;

%mend datagen;


ods listing close;
%datagen(t=5,nreps=4,nsims=200);

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

ods listing;

proc print data=covparms (obs=20);
   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=twofixed;
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 "output from PROC MIXED  macro - random effects";
   plot upper*sim=1 lower*sim=2/overlay vref=80;
run;

proc gplot data=estsfixed;
   title "output from PROC MIXED  macro - fixed effects";
   plot upper*sim lower*sim/overlay vref=80;
run;