newhnp {hnp}

Method for non-implemented model classes

Description

Uses user defined functions to produce the (half-)normal plot with simulated envelope.

Usage

newhnp(object, sim=99, conf=.95, halfnormal=T, plot.sim=T,
       verb.sim=F, how.many.out=F, print.on=F, paint.out=F,
       col.paint.out, diagfun, simfun, fitfun, ...)

Arguments

object

fitted model object or numeric vector.

sim

number of simulations used to compute envelope. Default is 99.

conf

confidence level of the simulated envelope. Default is 0.95.

halfnormal

logical. If TRUE, a half-normal plot is produced. If FALSE, a normal plot is produced. Default is TRUE.

plot.sim

logical. Should the (half-)normal plot be plotted? Default is TRUE.

verb.sim

logical. If TRUE, prints each step of the simulation procedure. Default is FALSE.

how.many.out

logical. If TRUE, the number of points out of the envelope is printed. Default is FALSE.

print.on

logical. If TRUE, the number of points out of the envelope is printed on the plot. Default is FALSE.

paint.out

logical. If TRUE, points out of the simulation envelope are plotted in a different color. Default is FALSE.

col.paint.out

If paint.out=TRUE, sets the color of points out of the envelope. Default is "red".

diagfun

user-defined function used to obtain the diagnostic measures from the fitted model object (only used when newclass=TRUE). Default is resid.

simfun

user-defined function used to simulate a random sample from the model estimated parameters (only used when newclass=TRUE).

fitfun

user-defined function used to re-fit the model to simulated data (only used when newclass=TRUE).

...

extra graphical arguments passed to plot.hnp.

Details

By providing three user-defined functions, newhnp produces the half-normal plot with simulated envelope for a model whose class is not yet implemented in the package.

The first function, diagfun, must extract the desired model diagnostics from a model fit object. The second function, simfun, must return the response variable (numeric vector or matrix), simulated using the same error distributions and estimated parameters from the fitted model. The third and final function, fitfun, must return a fitted model object. See the Examples section.

Author(s)

Rafael A. Moral <[e-mail]>, John Hinde and Clarice G. B. Demétrio

Examples

## 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))
##   }
##
## fitfun <- function(y.) {
##  userfunction(y. ~ linear_predictor, other_arguments, data=data)
##    # e.g., glm(y. ~ block + factor1 * factor2, family=poisson,
##    #           data=mydata)
##  }
##
## when response is binary:
## fitfun <- function(y.) {
##  userfunction(cbind(y., m-y.) ~ linear_predictor,
##               other_arguments, data=data)
##    #e.g., glm(cbind(y., m-y.) ~ treatment - 1,
##    #          family=binomial, data=data)
##  }

## Not run: 
## 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
my.data <- data.frame(y, tr)
f.fun <- function(y.) glm(y. ~ tr, family=poisson, data=my.data)

# hnp call
hnp(fit1, newclass=TRUE, diagfun=d.fun, simfun=s.fun, fitfun=f.fun)

## 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
my.data <- data.frame(y, tr)
f.fun <- function(y.) gamlss(y. ~ tr, family=GA, data=my.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
my.data <- data.frame(y, tr, m)
f.fun <- function(y.) gamlss(cbind(y., m-y.) ~ tr,
                               family=BI, data=my.data)

# hnp call
hnp(fit3, newclass=TRUE, diagfun=d.fun, simfun=s.fun, fitfun=f.fun)

## End(Not run)