VGAM (version 1.1-14)

calibrate.qrrvglm: Calibration for CQO and CAO models

Description

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

Usage

calibrate.qrrvglm(object, newdata = NULL,
    type = c("latvar", "predictors", "response", "vcov", "everything"),
    lr.confint = FALSE, cf.confint = FALSE,
    level = 0.95, initial.vals = NULL, ...)

Value

Several methods are implemented to obtain confidence intervals/regions for the calibration estimates. One method is when lr.confint = TRUE, then a 4-column matrix is returned with the confidence limits being the final 2 columns (if type = "everything" then the matrix is returned in the lr.confint list component). Another similar method is when cf.confint = TRUE. There may be some redundancy in whatever is returned. Other methods are returned by using type

and they are described as follows.

The argument type determines what is returned. If type = "everything" then all the type values are returned in a list, with the following components. Each component has length nrow(newdata).

latvar

Calibrated latent variables or site scores (the default). This may have the attribute "objectiveFunction" which is usually the log-likelihood or the deviance.

predictors

linear/quadratic or additive predictors. For example, for Poisson families, this will be on a log scale, and for binomial families, this will be on a logit scale.

response

Fitted values of the response, evaluated at the calibrated latent variables.

vcov

Wald-type estimated variance-covariance matrices of the calibrated latent variables or site scores. Actually, these are stored in a 3-D array whose dimension is c(Rank(object), Rank(object), nrow(newdata)). This type has only been implemented for binomialff and poissonff models with canonical links and noRRR = ~ 1 and, for CQOs, I.tolerances = TRUE or eq.tolerances = TRUE.

Details

Given a fitted regression CQO/CAO model, maximum likelihood calibration is theoretically easy and elegant. However, the method assumes that all the responses are independent, which is often not true in practice. More details and references are given in Yee (2018) and ch.6 of Yee (2015).

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

References

Abate, J. and Whitt, W. (1992). The Fourier-series method for inverting transforms of probability distributions. Queueing Systems, 10, 5--88.

ter Braak, C. J. F. (1995). Calibration. In: Data Analysis in Community and Landscape Ecology by Jongman, R. H. G., ter Braak, C. J. F. and van Tongeren, O. F. R. (Eds.) Cambridge University Press, Cambridge.

See Also

calibrate.qrrvglm.control, calibrate.rrvglm, calibrate, cqo, cao, optim, uniroot.

Examples

if (FALSE) {
hspider[, 1:6] <- scale(hspider[, 1:6])  # Stdze environmental variables
set.seed(123)
siteNos <- c(1, 5)  # Calibrate these sites
pet1 <- cqo(cbind(Pardlugu, Pardmont, Pardnigr, Pardpull, Zoraspin) ~
        WaterCon + BareSand + FallTwig + CoveMoss + CoveHerb + ReflLux,
        trace = FALSE,
        data = hspider[-siteNos, ],  # Sites not in fitted model
        family = poissonff, I.toler = TRUE, Crow1positive = TRUE)
y0 <- hspider[siteNos, colnames(depvar(pet1))]  # Species counts
(cpet1 <- calibrate(pet1, trace = TRUE, newdata = data.frame(y0)))
(clrpet1 <- calibrate(pet1, lr.confint = TRUE, newdata = data.frame(y0)))
(ccfpet1 <- calibrate(pet1, cf.confint = TRUE, newdata = data.frame(y0)))
(cp1wald <- calibrate(pet1, newdata = y0, type = "everything"))
}

if (FALSE) {
# Graphically compare the actual site scores with their calibrated
# values. 95 percent likelihood-based confidence intervals in green.
persp(pet1, main = "Site scores: solid=actual, dashed=calibrated",
      label = TRUE, col = "gray50", las = 1)
# Actual site scores:
xvars <- rownames(concoef(pet1))  # Variables comprising the latvar
est.latvar <- as.matrix(hspider[siteNos, xvars]) %*% concoef(pet1)
abline(v = est.latvar, col = seq(siteNos))
abline(v = cpet1, lty = 2, col = seq(siteNos))  # Calibrated values
arrows(clrpet1[,  3], c(60, 60), clrpet1[,  4], c(60, 60),  # Add CIs
       length = 0.08, col = "orange", angle = 90, code = 3, lwd = 2)
arrows(ccfpet1[,  3], c(70, 70), ccfpet1[,  4], c(70, 70),  # Add CIs
       length = 0.08, col = "limegreen", angle = 90, code = 3, lwd = 2)
arrows(cp1wald$latvar - 1.96 * sqrt(cp1wald$vcov), c(65, 65),
       cp1wald$latvar + 1.96 * sqrt(cp1wald$vcov), c(65, 65),  # Wald CIs
       length = 0.08, col = "blue", angle = 90, code = 3, lwd = 2)
legend("topright", lwd = 2,
       leg = c("CF interval", "Wald  interval", "LR interval"),
       col = c("limegreen", "blue", "orange"), lty = 1)
}