Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- library(gridBase)
- library(grid)
- library(mgcv)
- library(statmod)
- ## nicer version of gam.check/rqgam.chack
- # - deviance residuals for the Q-Q and histogram
- # - RQR for resids vs linear pred
- # - response vs fitted
- # hist.p gives the quantiles of the residuals to show
- better_check <- function(b, k.sample = 5000, k.rep = 200, rep = 0, level = 0.9,
- rl.col = 2, rep.col = "gray80", hist.p=c(0.01, 0.099), ...){
- # layout stuff
- opar <- par(mfrow=c(2,2))
- # grab the randomised quantile residuals
- # requires statmod package
- # need to do the right thing for mgcv's Tweedie
- if(grepl("^Tweedie",b$family$family)){
- if(is.null(environment(b$family$variance)$p)){
- p.val <- b$family$getTheta(TRUE)
- environment(b$family$variance)$p <- p.val
- }
- qres <- qres.tweedie(b)
- # and for negbin
- }else if(grepl("^Negative Binomial",b$family$family)){
- # need to set $theta
- if("extended.family" %in% class(b$family)){
- # for SNW's extended family, need to set TRUE in getTheta as theta
- # is on the wrong scale
- b$theta <- b$family$getTheta(TRUE)
- }else{
- b$theta <- b$family$getTheta()
- }
- qres <- qres.nbinom(b)
- }else{
- stop("Only negative binomial and Tweedie response distributions are supported.")
- }
- # values of the linear predictor
- linpred <- napredict(b$na.action, b$linear.predictors)
- ## get the deviance residuals
- resid <- residuals(b, type = "deviance")
- ## Q-Q plot of deviance resids
- qq.gam(b, rep = rep, level = level, type = "deviance", rl.col = rl.col,
- rep.col = rep.col, main="Quantile-quatile plot", ...)
- ## Q resids vs. linear pred
- plot(linpred, qres, main="Resids vs. linear pred.",
- xlab="linear predictor", ylab="Randomised quantile residuals",
- pch=19, cex=0.5, col=rgb(0, 0, 0, 0.6), ...)
- lines(lowess(linpred, qres), col = 2)
- ## histogram of deviance resids
- vps <- baseViewports()
- pushViewport(vps$figure)
- # partial histogram
- hist(resid[resid >= quantile(resid, hist.p[1]) &
- resid <= quantile(resid, hist.p)[2]],
- xlab = "deviance residuals", main = "Histogram of residuals",...)
- box()
- # setup the smaller plot
- pushViewport(viewport(x=-.4, y=.45 ,width=.3, height=.3))
- oopar <- par(plt = gridFIG(), new=TRUE)
- # histogram of full data
- hist(resid, xlab = "", main="", axes=FALSE, ylab="", ...)
- axis(1, cex.axis=0.5, mgp=c(3, 0.1, 0), tck=-0.02)
- # reset grid options for next plot
- popViewport(2)
- par(oopar)
- ## Response vs. Fitted Values
- plot(fitted(b), b$y,
- main="Response vs. Fitted Values",
- xlab="Fitted Values", ylab="Response",
- pch=19, cex=0.5, col=rgb(0, 0, 0, 0.6), ...)
- lines(lowess(b$fitted.values, b$y), col = 2)
- par(opar)
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement