newhnp: Method for non-implemented model classes
In hnp: Half-Normal Plots with Simulation Envelopes

Description Usage Arguments Details Author(s) Examples

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

1
2
3

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, ...)

`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`.

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.

Rafael A. Moral <rafael_moral@yahoo.com.br>, John Hinde and Clarice G. B. Demétrio

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