/* The point of this little simulation is to investigate the
small sample properties of confidence intervals for variance
components in the two-factor nested random effects model.
The methods we consider are 
1) a chisquare approximation to  the linear combo of MS terms
required to estimate a variance component

sigmahat * df/chisq(high) < sigma < sigmahat *df/chisq(low)

2) the asymptotic normality of the maximum likelihood estimate
of the  variance component with asymptotic standard error
estimated using the information matrix:

sigmahat +- 1.96 SE(sigmahat)

In particular, several small values for the number of levels
of the nested factor are considered */

options ls=75 nodate; 

%let seed=0;
%let sigma_container=1.9;  /* container varcomp=3.5 */
%let sigma=0.5;            /* error varcomp=.25 */
%let numsims=100;


%macro datagen(numcontainers);

data one; 
   array trt t1-t3 (4.0,4.8,9.3);
   array tname $ tn1-tn3 ("12b","6b","s");
   do simindex=1 to &numsims;
   do i=1 to 3;
   do container=1 to &numcontainers;
   B=rannor(&seed)*&sigma_container;
   do subsample=1 to 3;
            E=normal(&seed)*&sigma;
            truemean=trt{i};
            y=truemean+B+E;
            *salinity=tname{i};
            salinity=i;
            output;
         end;
      end;
   end;
   end;
run;

ods trace on;
ods listing close;
proc mixed method = type3 cl data=one;
   by simindex;
   class salinity container;
   model y=salinity;
   random container(salinity)/ cl;
   ods output covparms=varcomps;
run;

proc mixed cl data=one;
   by simindex;
   class salinity container;
   model y=salinity;
   random container(salinity)/ cl;
   ods output covparms=varcomps_wald;
run;

/*proc print data=varcomps;run;*/

data varcomps_wald;  set varcomps_wald;  method="wald"; run;
data varcomps;  set varcomps;  method="symm"; run;
data all;
   set varcomps varcomps_wald;
run;


ods listing ;

data sigma;
   set all;
   if CovParm="Residual";
   cover=(&sigma**2 > Lower) & (&sigma**2 < Upper); 
run; 

data sigma_container;
   set all;
   if CovParm="container(salinity)";
   cover=(&sigma_container**2 > Lower) & (&sigma_container**2 < Upper); 
run; 

data sigma;
   set sigma sigma_container;
run;
proc sort data=sigma;   by CovParm method simindex; run;
ods listing ;
proc freq data=sigma;
   title "number of containers is &numcontainers";
   by CovParm;
   table method*cover;
run;
/*
proc print data=sigma (obs=5);
   where cover=0;
run;
*/

%mend datagen;
%datagen(2);
%datagen(3);
%datagen(4);
%datagen(5);
%datagen(10);