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 ;
%MACRO u1kern1 ;
PROC IML ;
USE piest ;
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 ;
ehat = yvec-pihat ;
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 ;
dhat = betahat//muhat ;
* CREATE g vector and M matrix ;
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 Ob=1;
else Ob=0;
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 ;
