VGAM (version 1.0-4)

calibrate.rrvglm: Calibration for CLO models (RR-VGLMs)

Description

Performs maximum likelihood calibration for constrained and unconstrained quadratic and additive ordination models (CQO and CAO models are better known as QRR-VGLMs and RR-VGAMs respectively).

Usage

calibrate.rrvglm(object, newdata = NULL,
        type = c("latvar", "predictors", "response", "vcov", "all3or4"),
        initial.vals = NULL, ...)

Arguments

object

The fitted rrvglm model.

newdata

A data frame with new response data (e.g., new species data). The default is to use the original data used to fit the model; however, the calibration may take a long time to compute because the computations are expensive.

type

Same as calibrate.qrrvglm. The "all3or4" is for all of them, i.e., all types. For CLO models, "vcov" is unavailable, so all 3 are returned.

initial.vals

Same as calibrate.qrrvglm. The default is a grid defined by arguments in calibrate.rrvglm.control.

…

Arguments that are fed into calibrate.rrvglm.control.

Value

See calibrate.qrrvglm.

Warning

This function assumes that the prior weights are all unity; see weightsvglm.

This function is computationally expensive for Rank >= 1, and since it uses a for() loop several times it can be slow. Setting trace = TRUE to get a running log is a good idea.

Details

Given a fitted regression CLO model, maximum likelihood calibration is theoretically easy and elegant. However, the method assumes that all responses are independent. More details and references are given in Yee (2012).

The function optim is used to search for the maximum likelihood solution. Good initial values are needed, and calibrate.rrvglm.control allows the user some control over the choice of these.

It is found empirically that the stereotype model (also known as a reduced-rank multinomial logit model) calibrates well only with grouped data. When the response vector is all 0s except for one 1, then the MLE will probably be at -Inf or +Inf.

See Also

calibrate.qrrvglm, calibrate, rrvglm, weightsvglm.

Examples

# NOT RUN {
nona.xs.nz <- na.omit(xs.nz)  # Overkill!!
nona.xs.nz$dmd       <- with(nona.xs.nz, round(drinkmaxday))
nona.xs.nz$feethr    <- with(nona.xs.nz, round(feethour))
nona.xs.nz$sleephr   <- with(nona.xs.nz, round(sleep))
nona.xs.nz$beats     <- with(nona.xs.nz, round(pulse))

p2 <- rrvglm(cbind(dmd, feethr, sleephr, beats) ~ age + smokenow +
  depressed + embarrassed + fedup + hurt + miserable +  # 11 psychological
  nofriend + moody + nervous + tense + worry + worrier, # variables
  noRRR = ~ age + smokenow, trace = FALSE, poissonff, data = nona.xs.nz,
         Rank = 2)
cp2 <- calibrate(p2, new = head(nona.xs.nz, 9), type = "all", trace = TRUE)
cp2

two.cases <- nona.xs.nz[1:2, ]  # Another example
two.cases$dmd       <- c(4, 10)
two.cases$feethr    <- c(4, 7)
two.cases$sleephr   <- c(7, 8)
two.cases$beats     <- c(62, 71)
cp2b <- calibrate(p2, new = two.cases, type = "all")
cp2b
# }