*APPENDIX A. SAS Code for the example data analysis used
in this paper is available in ;
*electronic form. From Sturdivant et al. A smoothed residual based
goodness-of-fit statistic for ;
*nest-survival models. Studies in Avian Biology. ;
*
MACRO used to produce the USS Kernel Smoothed Statistic ;
read all var {ifate}
into yvec ; * RESPONSE VARIABLE
NAME HERE ;
read all var {pred} into
pihat ;
CLOSE piest ;
USE west ;
read all var {pred} into
wvec ;
CLOSE west ;
what = diag(wvec) ;
winv = inv(what) ;
n = nrow(pihat) ;
getwt = ceil(0.5*sqrt(n))+1 ; * SELECT THE BANDWIDTH HERE
;
* KERNEL SMOOTH ROUTINE ;
wtmat = J(n,n) ;
rx=J(n,1) ;
do
i=1 to n ;
x
= abs(pihat[i] - pihat);
rx[rank(x)]=x;
bw =
rx[getwt] ;
wtmat[,i] =
x / bw ;
end ;
* Get Kernel density values and
weights;
* UNIFORM (-a,a) ;
ukern = t(wtmat<1) ;
icolsum = 1/ukern[,+] ;
uwt
= ukern # icolsum ;
*
CUBIC ;
ctemp
= 1 - (t(wtmat))##3 ;
ckern
= ukern # ctemp ;
icolsum
= 1/ckern[,+] ;
cwt
= ckern # icolsum ;
*
NORMAL ;
nkern
= pdf('norm',t(2*wtmat)) ;
icolsum
= 1/nkern[,+] ;
nwt
= nkern # icolsum ;
*
MOMENTS and TEST STATISTICS;
USE mall ; * NAMES OF FIXED DESIGN MATRIX HERE and data set
;
read all var{lv3 sage
PctGr4} into x ; * Note: here lv3 is
all ones so used as int ;
read all var {site} into
groups ; * NAME OF LEVEL2
VARIABLE HERE ;
CLOSE mall ;
zmat = design(groups) ;
Q = x||zmat ;
USE betahat ;
read all var {Estimate}
into betahat ;
CLOSE betahat ;
USE Randeff ;
read all var {estimate}
into muhat ;
CLOSE Randeff ;
USE Sigmahat ;
read all var {estimate}
into cov2 ;
CLOSE Sigmahat ;
icov2 = 1/cov2 ;
icov2d=diag(icov2) ;
icov2a = BLOCK(icov2d,icov2d,icov2d)
;
icovbl2 =
BLOCK(icov2a,icov2a,icov2a,icov2a,icov2a,icov2a);
* BLOCKS SAME NUMBER AS GROUPS ;
faketop =
j(ncol(x),ncol(x)+ncol(zmat),0) ;
fakeleft = j(ncol(zmat),ncol(x),0)
;
comb1 = fakeleft||icovbl2 ;
R = faketop//comb1 ;
mymat = inv( t(Q)*what*Q + R ) ;
g = what * Q * mymat * R * dhat ;
M = what * Q * mymat * t(Q) ;
* CALCULATE TEST STATISTICS;
im = I(nrow(M))-M ;
midunif = t(uwt)*uwt ;
midcube = t(cwt)*cwt ;
midnorm = t(nwt)*nwt ;
aunif = t(im)*midunif*im ;
acube
= t(im)*midcube*im ;
anorm
= t(im)*midnorm*im ;
bunif
= 2*t(im)*midunif*g ;
bcube
= 2*t(im)*midcube*g ;
bnorm
= 2*t(im)*midnorm*g ;
Tuni
= t(ehat)*midunif*ehat ;
Tc
= t(ehat)*midcube*ehat ;
Tn
= t(ehat)*midnorm*ehat ;
*
CALCULATE EXPECTED VALUES ;
eunif = trace( aunif*what) +
t(g)*midunif*g ;
ecube = trace( acube*what) +
t(g)*midcube*g ;
enorm = trace( anorm*what) +
t(g)*midnorm*g ;
* CALCULATE VARIANCE ;
temp1 = wvec#(1-6*wvec)
;
temp3 = pihat#(1-pihat)#(1-2*pihat)
;
tempu = (vecdiag(aunif))##2 ;
tempc = (vecdiag(acube))##2 ;
tempn = (vecdiag(anorm))##2 ;
v1unif
= sum(tempu#temp1) ;
v2unif
= 2* trace(aunif*what*aunif*what) ;
v3unif = t(bunif)*what*bunif ;
v4unif = 2*sum(
(vecdiag(aunif))#bunif#temp3 ) ;
v1cube = sum(tempc#temp1) ;
v2cube = 2* trace(acube*what*acube*what)
;
v3cube = t(bcube)*what*bcube ;
v4cube = 2*sum( (vecdiag(acube))#bcube#temp3 )
;
v1norm
= sum(tempn#temp1) ;
v2norm
= 2* trace(anorm*what*anorm*what) ;
v3norm = t(bnorm)*what*bnorm ;
v4norm = 2*sum(
(vecdiag(anorm))#bnorm#temp3 ) ;
vunif = v1unif + v2unif + v3unif +
v4unif ;
vcube = v1cube + v2cube + v3cube + v4cube ;
vnorm = v1norm + v2norm + v3norm + v4norm ;
cubestat
= (Tc-ecube)/sqrt(vcube) ;
normstat
= (Tn-enorm)/sqrt(vnorm) ;
unifstat
= (Tuni-eunif)/sqrt(vunif) ;
punif
= 2*(1-probnorm(abs(unifstat))) ;
pcube
= 2*(1-probnorm(abs(cubestat))) ;
pnorm
= 2*(1-probnorm(abs(normstat))) ;
print
Tc ecube vcube cubestat pcube ;
print
Tn enorm vnorm normstat pnorm ;
print
Tuni eunif vunif unifstat punif ;
quit
; run ;
%MEND ;
*
NEST DATA ;
data
Mall; set Nests.mall2000nd;
LV3 =1 ;
* ADD A COLUMN OF
ONES FOR INTERCEPT ;
if ifate=0 then ness+1;
else if ifate=1 then ness+t;
/*
create indicator variables for different nesting habitats */
if hab=1 then NatGr=1; else
NatGr=0;
/* Native Grassland */
if hab=2 or hab=3 or hab=9
then CRP=1; else CRP=0; /* CRP & similar */
if hab=7 or hab=22 then Wetl=1;
else Wetl=0;
/* Wetland
sites */
if hab=20 then Road=1; else
Road=0;
/* Roadside sites */
run;
Proc Sort
data=Mall;
by
site; run;
* FIT MODEL USING PROC NLMIXED;
PROC NLMIXED
DATA=Mall tech=quanew method=gauss maxiter=1000;
parms
B0=2.42, B2=0.019, B4=0.38, s2u=0.1;
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~binary(p);
random
u~normal(0,s2u) subject=site out=randeff;
predict
p*(1-p) out=west;
predict
p out=piest ;
ods
output ParameterEstimates=betahat
(where=(Parameter=:"B")) ;
ods
output ParameterEstimates=sigmahat
(where=(Parameter=:"s2")) ;
ods
output ParameterEstimates=B0Hat
(where=(Parameter='B0') rename=(Estimate=Est_B0));
ods
output ParameterEstimates=B1Hat
(where=(Parameter='B1') rename=(Estimate=Est_B1));
ods
output ParameterEstimates=B2Hat
(where=(Parameter='B2') rename=(Estimate=Est_B2));
ods
output ParameterEstimates=s2uhat
(where=(Parameter='s2u') rename=(Estimate=Est_s2u));
run ;
%u1kern1
; * CALL KERNEL SMOOTHED STATISTIC MACRO
;