overdispersion in pcount models: QAIC?

[Suresh Andrew Sethi]

Dear ‘unmarked’ gurus,

I’m attempting to model the influence of habitat variables on
occurrence and detection for juvenile salmon in lake habitats using a
Royle N-mixture model for repeated minnow trapping counts (pcount() in
unmarked). The global model suggests overdispersion from a Poisson
state process model (likely from schooling of fish), as assessed by a
parametric bootstrap goodness of fit testing procedure using a Chi-
squared fit statistic (Fiske and Chandler 2011; MacKenzie and Bailey
2004).

From my understanding, there’s a few options to deal with
overdispersion in this context which could be implemented with
‘unmarked’:
1. Accommodate overdispersion by modeling the state process as
Negative Binomial.
2. Accommodate overdispersion using ZIP state process.

And my question to you:
3. Use a Poisson state process, but adjust model inference using a
variance inflation factor (e.g. QAIC, and assuming that one believes
their model is structurally sufficient to characterize expected
abundance but that they wish to get more honest estimates of
coefficient uncertainty).

I’ve implemented options 1 and 2, but they are not satisfactory in
this case (NB model passes par. boot. g.o.f. test but very sensitive
to choice of K parameter for integrated likelihood; ZIP model
demonstrates systematic upward bias in residuals for predicted
abundance whereas models with Poisson and NB state processes show no
such bias). In particular I would appreciate your comment on the use
of QAIC for Royle N-mixture models on repeated counts (single-season,
pcount() type model).

Could one extend the thinking of MacKenzie and Bailey (2004) and White
et al. (2002) to come up with an estimate of a (constant) variance
inflation factor for use in a repeated count N-mixture model
application using a parametric bootstrap routine (e.g. code provided
in Fiske and Chandler 2011) and a Chi square fit statistic as:
c.hat = (observed.chi.sq.stat)/mean(par.boot.chi.sq.stats)?
Regarding such an estimate of c.hat, White et al. (2002, pg. 375)
summarize the spirit of this approach: "The logic is that the mean of
the simulated deviances represents the expected value of the deviance
under the null model of no violations of assumptions (i.e., perfect
fit of the model to the data). Thus, cˆ = observed deviance divided by
expected deviance provides a measure of the amount of overdispersion
in the original data."

Thanks for your comments, and thanks for writing and maintaining the
very impressive 'unmarked' package.

Cheers,

Suresh Andrew Sethi
US Fish and Wildlife Service, Alaska

[Richard Chandler]

Hi Suresh,

Good question. I don't think there is a correct answer to this, but I see no reason why you couldn't compute QAIC and use it in the way that Burnham and Anderson (2002) and others recommend. However, when you do this, it is important to remember that you have left the realm of parametric statistics (so you probably shouldn't use functions in unmarked like simulate and parboot). If you are comfortable abandoning parametric statistics, then you could also consider the non-parametric bootstrap (see ?nonparboot) as alternative to adjusting standard errors using variance inflation factors.  

Richard

_____________________________________
Richard Chandler, post-doc
USGS Patuxent Wildlife Research Center
301-497-5696

From:	Suresh Andrew Sethi <sas...@gmail.com>
To:	unmarked <unma...@googlegroups.com>
Date:	04/16/2012 02:57 PM
Subject:	[unmarked] overdispersion in pcount models: QAIC?
Sent by:	unma...@googlegroups.com

---

[Kristen Kirkby]

Hello,

I’m in a similar position trying to get at overdispersion. I’m working with a pcount Poisson model for abundance and detection probability of juvenile coho salmon in pools of Oregon streams.  Estimates from the negative binomial models wouldn’t ever settle as I increased initial K values, so I’ve stuck with Poisson.  My models don’t seem to have great fit, and I’d like to be able to consider and account for overdispersion.  To get at this, I’ve calculated the deviance residuals from observed and expected (fitted function) counts. I then took the sum of the squared deviance residuals and divided by the degrees of freedom to calculate a dispersion parameter with which to adjust my AIC and SE values.  I’ve included the code below.  The equations I took from stats courses and books.

After running this code, I do get a dispersion parameter that’s roughly similar to dispersion from a quasipoisson glm on just the abundance model, and also similar to Suresh’s Chi-squared-derived dispersion parameter.  (Unfortunately, the dispersion parameters are all disconcertingly high, ranging from around 4.5-6.5 for models in three different streams.)  In any case, everything in the code seems to function just fine, but I’d very much like a second, more experienced opinion on the appropriateness of using this method for pcount Poisson models in unmarked.

Your thoughts and advice are much appreciated!  Please let me know if you have any questions for me.

Kristen

 

Code:

full.model<-pcount(~diver.h+areap.h+depthp.h+woody.h~depth.h+area.h+ddf.h+cd.h,K=350,h.umf)

k<-10    # number of parameters

observed <- getY(full.model@data)

expected <- fitted(full.model)

df<-(length(observed)-sum(is.na(observed)))-k  #degrees of freedom

 

dev.i<-(sign(observed-expected)*sqrt(2*(observed*log(observed/expected)-(observed-expected))))  # Individual deviance residuals

dev.i.sq<-dev.i^2    

sumsqdev<-sum(dev.i.sq,na.rm=T) #sum of squared residuals

dispersion<-sumsqdev/df

c.hat<-dispersion   # c.hat / dispersion parameter 

 

q.aic.m1<-(-2*logLik(full.model)/c.hat)+(2*k*(k+1)/(n-k-1)) #Quasi-AICc

 

# Adjusted SE and P-values

se.full.model<-SE(full.model)

se.adj<-((sqrt(c.hat))*se.full.model)

z.adj<-(coef(full.model)/se.adj)

p.adj<-2*pnorm(-abs(z.adj))

[Jeffrey Royle]

hi Sarah, just a couple of comments on this:

 If your species are failing the GoF test, and this is due to overdispersion of abundance, you might try the negative binomial model for N. We had an example in the unmarked workshop I think [and doing the GoF testing of that]. (a link to the workshop material can be found off the unmarked site).

 I don't know much about this c.hat business -- that always struck me as kind of ad hoc and rule-of-thumbish.  Maybe some other users can chime in about that.

 I guess the CI's for p-hat are > 1 because the CI's are based on asymptotic standard errors (just +/- 2*SE) -- so-called "Wald-type" intervals.   In small samples (or even moderate samples) the likelihood is not very symmetric and the asymptotic SE is not very accurate. (or at least 2*SE is not an accurate interval).

 

regards

andy

 

 

 

On Tue, Oct 23, 2012 at 10:45 AM, Sarah Thompson <thom...@umn.edu> wrote:

Hi Unmarked users --

I'm using distsamp and I've had some issues with models failing the GOF test. It seems like this will be an issue for quite a few people and it would be great to have some expert opinion on how to proceed.
--
2 of 4 of my species fail the GOF test on a highly saturated model -- I thought that a way to deal with this might be to adjust AIC and SEs via c.hat (QAIC ala ~ MacKenzie and Bailey) as Suresh suggested previously.

I was a little perplexed by Richard's suggestions in response to that question, but I tried to follow them :)

1.) In particular, if I have "left the realm of parametric statistics" & parboot is no longer appropriate -- was it OK to use parboot() in the first place to generate c.hat?

2.) And also in attempting to following Richard's suggestions, I tried to use nonparboot to generate bootstrapped SEs, but found that for some reason, calling (method=nonparboot) from within linearComb did not work for me (I got an error for unused argument).

So, I returned to QAIC adjustments.
It is easy enough to estimate a c.hat from Chi-square stats generated by parboot() and to adjust AIC by this value. MacKenzie and Bailey recommend then adjusting SEs of betas by sqrt(c.hat) and that is a little bit trickier for me.
3.) This statement will fully give away how much I am *not* a statistician...
I see this c.hat adjustment to SEs as a way to show that you are being cautious with your conclusions and inflating your SEs depending on the degree of model non-fit. I personally do not know how to adjust the parameter SEs from within predict() (--which is the method I am using to generate estimated abundance and detection values)... I am tempted to simply use predict() normally, then adjust the subsequent SEs generated around the predicted abundance/detection rate estimates. I suspect that is a major violation of some statistical rule .. but it sure would be easy.

Anyways, if anyone else is dealing with this or has sage advice, I'd love to hear it!

4.) Oh, and 1 other thing. I occasionally get estimates of detection rates that are greater than 100% (particularly with the upper CI). Does that indicate that I am doing something wrong or is it pretty common?

Sarah

[Richard Chandler]

Hi Sarah,

Can you post the various error messages you are running into with parboot and nonparboot? And can you show the code you are using to get detection probs > 1? 

Richard

[Sarah Thompson]

 Thanks so much for the quick feedback!  I wish I would have asked this weeks ago...

-------------
Re Andy - I'll give the neg binomial a try. I had tried it in the past, but it used to be *much* slower in previous unmarked versions and I don't recall why I ultimately moved away from neg bin/gdistsamp. I suspect that model fit was similar.. but i will confirm that.
And I definitely agree that c.hat seems very ad hoc...just not sure what else to do.
I think that my form of overdispersion has to do with clustered point counts on sites & that I *should* go with BUGS to deal with the random effect of site.. but I am kind of set on using unmarked for this one for the sake of timeliness.
------------
Re Richard --
I hate publicly posting code because I'm sure anyone with more skill will laugh at it...
:-\

This is what generated high detection rates -- let me know if you need more info.
#top model
cc_mix1k         <- distsamp(~ sDATE + YEAR ~ sGRASS1 + sLIT + sROB + sSHRUB +sTREE + sTREE1 +sGRASS1K +sTREE1K, umfccsp, output = "abund", keyfun="halfnorm")
# generating predictions for detection rates
date105090<-(quantile(covs3$sDATE, c(1, 5, 9)/10))
(ccsp_date<- data.frame(sDATE=c(date105090),
                        YEAR=c("2009", "2009","2009","2010","2010","2010", "2011","2011","2011"),
                        sGRASS1=0, sLIT=0, sROB=0, sSHRUB=0, sTREE=0,
                        sTREE1=0, sGRASS1K=0,
                        sTREE1K=0))
(dateccsppred <-predict(cc_mix1k, type = "det", newdata = ccsp_date, appendData = TRUE))
uci1 <-(dateccsppred$Predicted)+(dateccsppred$SE)
lci1 <- (dateccsppred$Predicted)-(dateccsppred$SE)
(ccspdet_date <-cbind(c(2009,2009,2009,2010,2010,2010, 2011,2011,2011),c(2, 11, 23), lci1, dateccsppred$Predicted, uci1))

# 10, 50, 90th percentile for date by year
#       year  day    lci1           est    uci1
#[1,] 2009   2    51.89033 57.78156 63.67278
#[2,] 2009  11   60.87274 69.14673 77.42071
#[3,] 2009  23   69.45202 86.11592 102.77982 #
#[4,] 2010   2    65.73713 81.25977 96.78241
#[5,] 2010  11   77.09029 97.24291 117.39554 #
#[6,] 2010  23   89.25504 121.10715 152.95926  ## yikes
#[7,] 2011    2   48.11647 55.69513 63.27379
#[8,] 2011  11   58.15907 66.64992 75.14076
#[9,] 2011  23   68.21210 83.00637 97.80063


#-----------------nonparboot---------
#works fine without the method=nonparboot ---
b_mix500a  <- distsamp(~ sDATE ~ sGRASS1 + sGRASS1K, umfbobo, output = "abund",keyfun="halfnorm")
b_mix500a<-nonparboot(b_mix5002, B=100)
linearComb(b_mix500a, matrix(c(1, 0.5, 0.5,  1, 0, 0,  1, 0, 0.5), 3, 3, byrow=TRUE),
           type = "state", method = "nonparboot")
#error message:
#Error in .local(obj, coefficients, ...) :
# unused argument(s) (method = "nonparboot")
#-----------------------------------------------------

[Richard Chandler]

Hi Sarah,

The predictions you are getting are for the scale parameter "sigma" of the half-normal distribution. This parameter determines how quickly detection probability decreases with distance. It is not the same thing as detection probability, but you can compute detection probability once you have an estimate of sigma. The easiest way to do this is using the getP() function.

With respect to nonparboot, are you trying to get bootstrapped SEs for predictions? There might be a way of doing that by following the example on the help page for nonparboot. Another option (that is kind of ugly but makes more sense to me) is to make predictions from every bootstrap replicate and then compute the SDs of the bootstrap distribution. Something like this:  

sigma.hats <- sapply(b_mix500a@bootstrapSamples, function(x) predict(x, type = "det", newdata = ccsp_date)[,1])

apply(sigma.hats, 1, sd) # SEs for the predictions

Richard

[Sarah Thompson]

Hi Richard,

I figured that my issue with the high detection estimates was some sort of error; glad I asked!
I totally forgot that it was giving me sigmas with predict().

And yep, with nonparboot, I'm trying to get SEs for predictions... I will try your suggestions --

Thanks!!

Sarah