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