Code for analyzing fixed- and random-effects 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. ;
* This file:
* 1. inputs a dataset containing nest survival information,
* 2. calculates the effective sample size for the dataset,
* 3. runs a variety of models of NEST SURVIVAL in NLMIXED,
* 4. creates an AICc table for model selection, &
* 5. outputs the AICc table to HTML, RTF, and pdf files.
* Step 1: calculate the effective sample size according to the methods
* of Dinsmore et al. (2002). Here, n-ess is incremented by 1 for
* each day a nest was under observation and survived and by 1 for
* each interval for which a nest was under observation and failed.;
* This step:
* 1. calculates the contribution to n-ess for each observation
* interval and adds that contribution to the sum of ness, i.e.,
* ness column is a running total
* 2. creates dummy/indicator variables for each of 4 habitat types;
data Mall; set Sasuser.mall2000nd;
if ifate=0 then ness+1;
else if ifate=1 then ness+t;
/* create indicator variables for different nesting habitats */
/* Native Grassland */
if hab=1 then NatGr=1; else NatGr=0;
/* CRP & similar */
if hab=2 or hab=3 or hab=9 then CRP=1; else CRP=0;
/* Wetland sites */
if hab=7 or hab=22 then Wetl=1; else Wetl=0;
/* Roadside sites */
if hab=20 or hab=8 then Road=1; else Road=0;
run;
* This step finds the actual n-ess for the dataset,
* which is the maximum value in the ness column.;
Proc Univariate data=Mall;
var ness;
output out=ness max=ness;
run;
* This step sorts the data by site, which is used as a random factor in
* some models below. PROC NLMIXED assumes that a new realization
* occurs whenever the SUBJECT= variable changes from the previous
* observation, so your input data set should be clustered according
* to this variable. You can accomplish this by running PROC SORT
* prior to calling PROC NLMIXED using the SUBJECT=variable as the
* BY variable. ;
Proc Sort data=Mall;
by site; run;
* This step reformats the Fit Statistics table of NLMIXED
* so it displays more decimal places in the created tables;
proc template;
define table Stat.Nlm.FitStatistics;
notes "Fit statistics";
column Descr Value;
header H1;
define H1;
text "Fit Statistics";
space = 1;
end;
define Descr;
header = "Description";
width = 30;
print_headers = OFF;
flow;
end;
define Value;
header = "Value";
format = 12.4;
print_headers = OFF;
end;
end;
run;
* Run a fixed-effects model including nest age and the extent of
* grassland on the study site;
Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=0, B1=0, B2=0;
p=1;
do i=0 TO t-1;
logit=B0+B1*(SAge+i)+B2*PpnGr;
p=p*(exp(logit)/(1+exp(logit)));
end;
model ifate~binomial(1,p);
RUN;
* Run a fixed-effects model including nest age, the extent of
* grassland on the study site, and an observer effect on dsr for
* day 1 of each interval (done with a dummy variable called ‘Ob’);
Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=0, B2=0, B4=0, B8=0;
p=1;
do i=0 TO t-1;
if i=0 then Ob=1;
else Ob=0;
logit=(B0+(B8*Ob))+B2*(SAge+i)+B4*PpnGr;
p=p*(exp(logit)/(1+exp(logit)));
end;
model ifate~binomial(1,p);
RUN;
* Run a mixed model including nest age (fixed effect), the extent of
* grassland on the study site (fixed effect), and
* a random effect of site (random effect influences the intercept).;
Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=2.42, B2=0.019, B4=0.38, vsite=0.5;
p=1;
do i=0 TO t-1;
if i=0 then Ob=1;
else Ob=0;
logit=(B0+u)+B2*(sage+i)+B4*PctGr4;
p=p*(exp(logit)/(1+exp(logit)));
end;
model ifate~binomial(1,p);
random u~normal(0,vsite) subject=site;
RUN;
* Run a mixed model including nest age (fixed effect), the extent of
* grassland (fixed effect) on the study site, an observer effect on
* dsr for day 1 (fixed effect), and a random effect of site.;
Proc Nlmixed data=Mall tech=quanew method=gauss maxiter=1000;
parms B0=2.42, B2=0.019, B4=0.38, B8=-1, vsite=0.5;
p=1;
do i=0 TO t-1;
if i=0 then Ob=1;
else Ob=0;
logit=(B0+u+B8*Ob)+B2*(sage+i)+B4*PctGr4;
p=p*(exp(logit)/(1+exp(logit)));
end;
model ifate~binomial(1,p);
random u~normal(0,vsite) subject=site;
* The following lines of code estimate the daily survival rate
* (DSR) for a 15-day old nest on a site with a grassland
* proportion of 0.5 on (1) the day of a nest visit and
* (2) a day without a nest visit;
estimate 'dsr-visited'
exp(B0+B2*(15)+B4*0.5+B8*1)/(1+exp(B0+B2*(15)+B4*0.5+B8*1));
estimate 'dsr-not visited'
exp(B0+B2*(15)+B4*0.5+B8*0)/(1+exp(B0+B2*(15)+B4*0.5+B8*0));
RUN;
