* APPENDIX 2.  Simulation code for analyzing performance of generalized linear mixed models of nest-;

*  survival data from periodic nest visits with Proc NLMIXED in SAS (SAS Institute, 2000).;

* From Rotella et al. Extending methods for modeling heterogeneity in nest-survival        ;

* data using generalized mixed models.  Studies in Avian Biology.                                    ;

    * -- to avoid the problem of SAS Log Window becomes full;

*PROC PRINTTO LOG='C:\Windows\Temp\LOGFILE.Txt';

RUN;

* beginning of the macro program 'RNDM';

%MACRO RPT (reps=,numyrs=,numsites=,numnests=,trueB0=,trueB1=,s2site=);    

* The macro parameters are:;

*    # of reps, # of yrs, # of sites, # of nests, true Beta values, and true random effects (s**2);

%Do L=1 %to &reps;      * how many replications are to be done;

DATA nests;

*ods listing close;

/***  set parameters  ***/

retain seed_site (2); * these can be changed to get different simulations;

 ssite=sqrt(&s2site);  *calculate the std dev for the random effect of site;

 site=0;                        * initializes the counter for site number;

        /***  end of parameters ***/

DO yr=1 TO &numyrs;    * number of years per dataset;

      DO s=1 TO &numsites;          * number of sites per dataset;

        site=site+1;                * because of this, a different set of sites are simulated each year;

                                          * This may not be desired.;

        call rannor(seed_site,x2);

        siterndm=0+ssite*x2;

            DO n=1 TO &numnests;    * number of nests per site;

              nestid=n;

                  init=round((ranuni(0)*59)+1);

                  agestart=round((ranuni(0)*33)+1);

                  ppngrass=ranuni(0);

                  output;

end; end; end;

run;

data intervals;

 set nests;

do j=1 to &numnests;

      interval=0;

      datestart=init;

      fate=1;

      do until (check=0);

                  interval=interval+1;

                  a=ranuni(0);

                  if a <0.25 then intleng=5;

                  if 0.25 < a < 0.5 then intleng=6;

                  if 0.5 < a < 0.75 then intleng=7;

                  if a > 0.75 then intleng=8;

*                       drop a;

                  sum=agestart+intleng;

                  if sum>35 then intleng=(35-agestart);

                        do k=0 TO intleng-1;

*                             if i=0 then Ob=1;

*                             else Ob=0;

                              logit=(&trueB0+&trueB1*ppngrass+siterndm);

                              dsr=exp(logit)/(1+exp(logit));

                              fate=fate*ranbin(0,1,dsr);

                        end;

                  ageend=agestart+intleng;

                  dateend=datestart+intleng;

                  drop j k sum ;

                  output ;

                  agestart=ageend;

                  datestart=dateend;

                  if agestart>=35 then x=0; else x=1;

                  if fate=0 then z=0; else z=1;

                  check=x*z; if check=0 then return;

      end;

end;

run;

Proc Nlmixed tech=quanew method=gauss maxiter=1000;

parms B0=0, B1=0;

        p=1;

           do i=0 TO intleng-1;

           if i=0 then Ob=1;

           else Ob=0;

                   logit=B0+B1*ppngrass;

              p=p*(exp(logit)/(1+exp(logit)));

           end;

model fate~binomial(1,p);

ods output ParameterEstimates=B0Hat

         (where=(Parameter='B0') rename=(Estimate=Est_B0 StandardError=SE_B0 Lower=LCI_B0 Upper=UCI_B0));

ods output ParameterEstimates=B1Hat

         (where=(Parameter='B1') rename=(Estimate=Est_B1 StandardError=SE_B1 Lower=LCI_B1 Upper=UCI_B1));

RUN;

data ests;

  merge B0Hat B1Hat;

    B0True=&TrueB0; B1True=&TrueB1;

      if LCI_B0 <= &TrueB0 <= UCI_B0 then CovB0=1; else covB0=0;

    if LCI_B1 <= &TrueB1 <= UCI_B1 then CovB1=1; else covB1=0;

      label SE_B0="SE_B0"; label LCI_B0="LCI_B0"; label UCI_B0="UCI_B0";

      label SE_B1="SE_B1"; label LCI_B1="LCI_B1"; label UCI_B0="UCI_B1";

      keep B0True Est_B0 SE_B0 LCI_B0 UCI_B0 CovB0 B1True Est_B1 SE_B1 LCI_B1 UCI_B1 CovB1;

run;

data _Null_;

 set ests;

 proc append base=work.NullModel;

run;

* Run the random-effects of SITE model;

Proc sort data=intervals;           * data must be sorted by the random effect variable;

 by site;

run;

Proc Nlmixed data=intervals tech=quanew method=gauss maxiter=1000;

parms B0=4, B1=.25, vsite=.4;

        p=1;

           do i=0 TO intleng-1;

           if i=0 then Ob=1;

           else Ob=0;

                   logit=B0+B1*ppngrass+u;

              p=p*(exp(logit)/(1+exp(logit)));

           end;

model fate~binomial(1,p);

random u~normal(0,vsite) subject=site;

ods output ParameterEstimates=B0Hat2

         (where=(Parameter='B0') rename=(Estimate=Est_B0 StandardError=SE_B0 Lower=LCI_B0 Upper=UCI_B0));

ods output ParameterEstimates=B1Hat2

         (where=(Parameter='B1') rename=(Estimate=Est_B1 StandardError=SE_B1 Lower=LCI_B1 Upper=UCI_B1));

ods output ParameterEstimates=vsiteHat

         (where=(Parameter='vsite') rename=(Estimate=Est_vsite StandardError=SE_vsite Lower=LCI_vsite Upper=UCI_vsite));

RUN;

data ests2;

  merge B0Hat2 B1Hat2 vsiteHat;

    B0True=&TrueB0; B1True=&TrueB1; s2siteTrue=&s2site;

      if LCI_B0 <= &TrueB0 <= UCI_B0 then CovB0=1; else CovB0=0;

    if LCI_B1 <= &TrueB1 <= UCI_B1 then CovB1=1; else covB1=0;

      if LCI_vyr <= s2siteTrue <= UCI_vsite then Covvsite=1; else Covvsite=0;

      label SE_B0="SE_B0"; label LCI_B0="LCI_B0"; label UCI_B0="UCI_B0";

      label SE_B1="SE_B1"; label LCI_B1="LCI_B1"; label UCI_B0="UCI_B1";

      label SE_vsite="SE_vsite"; label LCI_vsite="LCI_vsite"; label UCI_vsite="UCI_vsite";

      keep B0True     Est_B0    SE_B0    LCI_B0    UCI_B0    CovB0

             B1True     Est_B1    SE_B1    LCI_B1    UCI_B1    CovB1

             s2siteTrue Est_vsite SE_vsite LCI_vsite UCI_vsite Covvsite;

run;

data _Null_;

 set ests2;

 proc append base=work.v_Site_Model;

run;

%End;

%MEND RPT ;

%RPT (reps=100,numyrs=1,numsites=15,numnests=25,trueB0=2,trueB1=1.75,s2site=.25);   Run;