*
APPENDIX 1. 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;
*
* 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;
*
* grassland on the study site, and an observer
effect on dsr for
* day 1 of each interval (done with a dummy
variable called ‘
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
else
logit=(B0+(B8*
p=p*(exp(logit)/(1+exp(logit)));
end;
model
ifate~binomial(1,p);
RUN;
*
* 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
else
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;
*
* 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
else
logit=(B0+u+B8*
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;