ver Hoef and Boveng (2007) plots

verhoef.R

# 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]]))
}