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 ;