## Simple Poisson regression
set.seed(100)
counts <- c(rpois(5, 2), rpois(5, 4), rpois(5, 6), rpois(5, 8))
treatment <- gl(4, 5)
fit <- glm(counts ~ treatment, family=poisson)
anova(fit, test="Chisq")
## half-normal plot
hnp(fit)
## or save it in an object and then use the plot method
my.hnp <- hnp(fit, print.on=TRUE, plot=FALSE)
plot(my.hnp)
## changing graphical parameters
plot(my.hnp, lty=2, pch=4, cex=1.2)
plot(my.hnp, lty=c(2,3,2), pch=4, cex=1.2, col=c(2,2,2,1))
plot(my.hnp, main="Half-normal plot", xlab="Half-normal scores",
ylab="Deviance residuals", legpos="bottomright")
## Using a numeric vector
my.vec <- rnorm(20, 4, 4)
hnp(my.vec) # using N(0,1)
hnp(my.vec, scale=TRUE) # using N(mu, sigma^2)
## Implementing new classes
## Users provide three functions - diagfun, simfun and fitfun,
## in the following way:
##
## diagfun <- function(obj) {
## userfunction(obj, other_argumens)
## # e.g., resid(obj, type="pearson")
## }
##
## simfun <- function(n, obj) {
## userfunction(n, other_arguments) # e.g., rpois(n, fitted(obj))
## }
##
## when response is categorical (n>=3), "userfunction"
## (e.g., rmultinom) must not depend on n and n must be
## equal to 1 (unless "userfunction" works differently
## than rmultinom in a sense that it produces less than
## an entire simulated data set when n=1)
##
## fitfun <- function(data) {
## userfunction(y. ~ linear_predictor, other_arguments, data=data)
## # e.g., glm(y. ~ block + factor1 * factor2, family=poisson,
## # data=data)
## }
##
## the explanatory variables used in the linear predictor must
## be present in the data.frame passed to the "data" argument
## in hnp.
##
## when response is binary:
## fitfun <- function(data) {
## userfunction(cbind(y., m-y.) ~ linear_predictor,
## other_arguments, data=data)
## #e.g., glm(cbind(y., m-y.) ~ treatment - 1,
## # family=binomial, data=data)
## }
## Example no. 1: Using Cook's distance as a diagnostic measure
y <- rpois(30, lambda=rep(c(.5, 1.5, 5), each=10))
tr <- gl(3, 10)
fit1 <- glm(y ~ tr, family=poisson)
# diagfun
d.fun <- function(obj) cooks.distance(obj)
# simfun
s.fun <- function(n, obj) {
lam <- fitted(obj)
rpois(n, lambda=lam)
}
# fitfun
f.fun <- function(data) glm(y. ~ tr, family=poisson, data=data)
# hnp call
hnp(fit1, newclass=TRUE, diagfun=d.fun, simfun=s.fun,
fitfun=f.fun, data=data.frame(y, tr))
## Example no. 2: Implementing gamma model using package gamlss
# load package
require(gamlss)
# model fitting
y <- rGA(30, mu=rep(c(.5, 1.5, 5), each=10), sigma=.5)
tr <- gl(3, 10)
fit2 <- gamlss(y ~ tr, family=GA)
# diagfun
d.fun <- function(obj) resid(obj) # this is the default if no
# diagfun is provided
# simfun
s.fun <- function(n, obj) {
mu <- obj$mu.fv
sig <- obj$sigma.fv
rGA(n, mu=mu, sigma=sig)
}
# fitfun
f.fun <- function(data) gamlss(y. ~ tr, family=GA, data=data)
# hnp call
hnp(fit2, newclass=TRUE, diagfun=d.fun, simfun=s.fun,
fitfun=f.fun, data=data.frame(y, tr))
## Example no. 3: Implementing binomial model in gamlss
# model fitting
y <- rBI(30, bd=50, mu=rep(c(.2, .5, .9), each=10))
m <- 50
tr <- gl(3, 10)
fit3 <- gamlss(cbind(y, m-y) ~ tr, family=BI)
# diagfun
d.fun <- function(obj) resid(obj)
# simfun
s.fun <- function(n, obj) {
mu <- obj$mu.fv
bd <- obj$bd
rBI(n, bd=bd, mu=mu)
}
# fitfun
f.fun <- function(data) gamlss(cbind(y., m-y.) ~ tr,
family=BI, data=data)
# hnp call
hnp(fit3, newclass=TRUE, diagfun=d.fun, simfun=s.fun,
fitfun=f.fun, data=data.frame(y, tr, m))Run the code above in your browser using DataLab