/* This program uses features in
   SAS GLMMIX macro to fit logistic
   regression models to nest predation
   data from a split-plot experiment. 
   The code is stored in the file

          nestglmm.sas

   We currently do not have permission 
   to give you the data.          */


options linesize=64;

data set1;
  infile 'c:\mydocuments\courses\st557\sas\nestall.dat';
input wshed $ loc $ round roadside $ 
        roadtype $ adjhab $ rdepth rwidth
        foreback percip mtemp ftotal floss
        btotal bloss ctotal closs;
   if(roadtype='grav') then x1=1; else x1=0;
   if(adjhab='nrc') then x2=1; else x2=0;
   if(roadside='nfen') then x3=1; else x3=0;
   if(roadside='wd') then x4=1; else x4=0;
   x5=x1*x2; x6=x2*x4; x7=x1*x3; 
   slope=0; loss=floss; total=ftotal; output;
   slope=1; loss=bloss; total=btotal; output;
   keep wshed loc x1 x2 x3 x4 x5 x6 x7 
        slope loss total;
   run;


   data set1; set set1;
      slopec=slope; run;

/*  Fit a model use GLMMIX to include
    a random location effect  */

%inc  'c:\mydocuments\courses\st557\sas\glmm800.sas' /
            nosource;
            run;

/*  Logistic regression with random
    effects for transect and sub-transects */


%glimmix(data=set1,
   stmts=%str(
     class loc ;
     model loss/total = x1 x2 x3 x4 
         x5 x6 x7 slope  /
        solution  ddfm=kr;
     random intercept slope / subject=loc g gcorr;),
	 error=binomial, link=logit, 
     converge=1e-8, maxit=20, out=set
   )
run; 


proc print data=setp (obs=5); run;


/*  Logistic regression with heterogeneous
    extra-binomial variance adjustments
    for the foreslope and backslope*/

%glimmix(data=set1,
   stmts=%str(
     class loc ;
     model loss/total = x1 x2 x3 x4 
         x5 x6 x7 slope  / 
        solution  ddfm=kr;
     repeated / type=un subject=loc r rcorr; ),
	 error=binomial, link=logit, 
     converge=1e-8, maxit=20, out=setp    
   )
run;


proc print data=setp (obs=5); run;