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).
1 2 3 4
The fitted CQO/CAO model.
A data frame with new response data, such as 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.
Note that the creation of the model frame associated with
What type of result to be returned.
The first are the calibrated latent variables or site scores.
This is always computed.
likelihood ratio based confidence intervals?
characteristic function based confidence intervals?
Initial values for the search.
For rank-1 models, this should be a vector having length
Arguments that are fed into
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).
optim is used to search for
the maximum likelihood solution. Good initial values are
needed, and arguments in
allows the user some control over the choice of these.
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
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
and they are described as follows.
type determines what is returned.
type = "everything" then all the
are returned in a list, with the following components.
Each component has length
Calibrated latent variables or site scores
This may have the attribute
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.
Fitted values of the response, evaluated at the calibrated latent variables.
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
This function is computationally expensive.
trace = TRUE to get a running log can be a good idea.
This function has been tested but not extensively.
Despite the name of this function, CAO models are handled as well
to a certain extent.
Some combinations of parameters are not handled, e.g.,
lr.confint = TRUE only works for rank-1,
type = "vcov" only works for
models with canonical links and
noRRR = ~ 1,
and higher-order rank models need
eq.tolerances = TRUE or
I.tolerances = TRUE
For rank-1 objects,
lr.confint = TRUE is recommended
type = "vcov" in terms of accuracy and overall generality.
"qrrvglm" objects it is necessary that
all response' tolerance matrices are positive-definite
which correspond to bell-shaped response curves/surfaces.
deviance slot is used for the optimization rather than
loglikelihood slot, therefore one can calibrate using
real-valued responses. (If the
loglikelihood slot were used
then functions such as
dpois would be used
log = TRUE and then would be restricted to feed in
integer-valued response values.)
Maximum likelihood calibration for
Gaussian logit regression models may be performed by
rioja but this applies to a single environmental variable
data("SWAP", package = "rioja").
calibrate() estimates values of the
latent variable rather than individual explanatory variables,
hence the setting is more on ordination.
T. W. Yee.
Recent work on the standard errors by
David Zucker and
Sam Oman at HUJI
is gratefully acknowledged—these are returned in the
vcov component and provided inspiration for
A joint publication is being prepared on this subject.
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
## Not run: 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")) ## End(Not run) ## Not run: # 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) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.