Skip to content

Instantly share code, notes, and snippets.

@dill
Created August 28, 2015 03:18
Show Gist options
  • Save dill/c375e99cdbbd33819073 to your computer and use it in GitHub Desktop.
Save dill/c375e99cdbbd33819073 to your computer and use it in GitHub Desktop.
ver Hoef and Boveng (2007) plots
# makes plots (somewhat) like those in ver Hoef and Boveng (2007)
# Jay M. Ver Hoef and Peter L. Boveng 2007. QUASI-POISSON VS. NEGATIVE BINOMIAL REGRESSION: HOW SHOULD WE MODEL OVERDISPERSED COUNT DATA? Ecology 88:2766–2772. http://dx.doi.org/10.1890/07-0043.1
# http://www.utstat.utoronto.ca/reid/sta2201s/QUASI-POISSON.pdf
# calling with something like:
# par(mfrow=c(1,2))
# # define the breaks
# br <- c(seq(0,10,by=1/2))
# verhoef(mod0nbxy, br, "Negative binomial")
# verhoef(mod0qpxy, br, "Quasi-poisson")
# will give you what you want
verhoef <- function(model, breaks, main="", rtype="pearson"){
# calculate (y -mu)^2
av_sq_resids <- (model$y - residuals(model, type=rtype))^2
cat("av sq residual range=",range( av_sq_resids ),"\n")
# get linear predictor
lin_pred <- fitted(model)
cat("linear predictor range=",range(lin_pred),"\n")
# function make the binned data
mk_pts <- function(x, y){
# cut the x's
cx <- cut(x, breaks=breaks, include.lowest=TRUE)
# aggregate y's by the x's
cy <- aggregate(y, list(cx), FUN=mean)
# get the grouping names and make them into numbers
cx <- gsub("[\\[\\]\\(\\)]", "", as.character(cx),perl=TRUE)
cx <- unlist(lapply(strsplit(cx, ","), function(x) mean(as.numeric(x))))
# return the table of the x's and the y's
return(list(table(cx),cy))
}
# get the points to plot
pts <- mk_pts(av_sq_resids, lin_pred)
# use the model variance function to get the line
mod_var <- model$family$variance(fitted(model))
# actually make a plot
plot(sort(lin_pred), sort(mod_var), type="l",
xlab="Linear pred", ylab="Variance", main=main)
points(as.numeric(names(pts[[1]])), pts[[2]]$x, cex=as.numeric(pts[[1]]))
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment