vglm fit or a variance-covariance matrix,
and preprocesses it for
uninormal so that quasi-variances can be computed.
1 2 3
A single integer from the set
Character. Optional, for labelling the variance-covariance matrix.
Optional, passed into
Character. Label for for the reference level.
Optional numeric vector of length at least 3 specifying the indices of the factor from the variance-covariance matrix.
an optional vector of estimated coefficients
Suppose a factor with L levels is an explanatory variable in a
regression model. By default, R treats the first level as baseline so
that its coefficient is set to zero. It estimates the other L-1
coefficients, and with its associated standard errors, this is the
conventional output. From the complete variance-covariance matrix one
can compute L quasi-variances based on all pairwise difference
of the coefficients. They are based on an approximation, and can be
treated as uncorrelated. In minimizing the relative (not absolute)
errors it is not hard to see that the estimation involves a RCIM
rcim) with an exponential link function
object is a model, then at least one of
coef.indices must be non-
NULL. The value of
coef.indices, if non-
NULL, determines which rows and
columns of the model's variance-covariance matrix to use. If
coef.indices contains a zero, an extra row and column are
included at the indicated position, to represent the zero variances
and covariances associated with a reference level. If
factorname should be
the name of a factor effect in the model, and is used in order to
extract the necessary variance-covariance estimates.
Quasi-variances were first implemented in R with qvcalc. This implementation draws heavily from that.
A L by L matrix whose i-j element is the logarithm of the variance of the ith coefficient minus the jth coefficient, for all values of i and j. The diagonal elements are abitrary and are set to zero.
The matrix has an attribute that corresponds to the prior
weight matrix; it is accessed by
and replaces the usual
vglm. This weight matrix has ones on
the off-diagonals and some small positive number on
Negative quasi-variances may occur (one of them and
only one), though they are rare in practice. If
so then numerical problems may occur. See
qvcalc() for more information.
This is an adaptation of
qvcalc() in qvcalc.
It should work for all
models with one linear predictor, i.e., M = 1.
For M > 1 the factor should appear only in one of the
It is important to set
maxit to be larger than usual for
rcim since convergence is slow. Upon successful
convergence the ith row effect and the ith column effect
should be equal. A simple computation involving the fitted and
predicted values allows the quasi-variances to be extracted (see
A function to plot comparison intervals has not been written here.
T. W. Yee, based heavily on
qvcalc() in qvcalc
written by David Firth.
Firth, D. (2003). Overcoming the reference category problem in the presentation of statistical models. Sociological Methodology 33, 1–18.
Firth, D. and de Menezes, R. X. (2004). Quasi-variances. Biometrika 91, 65–80.
Yee, T. W. and Hadi, A. F. (2014). Row-column interaction models, with an R implementation. Computational Statistics, 29, 1427–1445.
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 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
# Example 1 data("ships", package = "MASS") Shipmodel <- vglm(incidents ~ type + year + period, poissonff, offset = log(service), # trace = TRUE, model = TRUE, data = ships, subset = (service > 0)) # Easiest form of input fit1 <- rcim(Qvar(Shipmodel, "type"), uninormal("explink"), maxit = 99) qvar(fit1) # Easy method to get the quasi-variances qvar(fit1, se = TRUE) # Easy method to get the quasi-standard errors (quasiVar <- exp(diag(fitted(fit1))) / 2) # Version 1 (quasiVar <- diag(predict(fit1)[, c(TRUE, FALSE)]) / 2) # Version 2 (quasiSE <- sqrt(quasiVar)) # Another form of input fit2 <- rcim(Qvar(Shipmodel, coef.ind = c(0, 2:5), reference.name = "typeA"), uninormal("explink"), maxit = 99) ## Not run: qvplot(fit2, col = "green", lwd = 3, scol = "blue", slwd = 2, las = 1) # The variance-covariance matrix is another form of input (not recommended) fit3 <- rcim(Qvar(cbind(0, rbind(0, vcov(Shipmodel)[2:5, 2:5])), labels = c("typeA", "typeB", "typeC", "typeD", "typeE"), estimates = c(typeA = 0, coef(Shipmodel)[2:5])), uninormal("explink"), maxit = 99) (QuasiVar <- exp(diag(fitted(fit3))) / 2) # Version 1 (QuasiVar <- diag(predict(fit3)[, c(TRUE, FALSE)]) / 2) # Version 2 (QuasiSE <- sqrt(quasiVar)) ## Not run: qvplot(fit3) # Example 2: a model with M > 1 linear predictors ## Not run: require("VGAMdata") xs.nz.f <- subset(xs.nz, sex == "F") xs.nz.f <- subset(xs.nz.f, !is.na(babies) & !is.na(age) & !is.na(ethnicity)) xs.nz.f <- subset(xs.nz.f, ethnicity != "Other") clist <- list("sm.bs(age, df = 4)" = rbind(1, 0), "sm.bs(age, df = 3)" = rbind(0, 1), "ethnicity" = diag(2), "(Intercept)" = diag(2)) fit1 <- vglm(babies ~ sm.bs(age, df = 4) + sm.bs(age, df = 3) + ethnicity, zipoissonff(zero = NULL), xs.nz.f, constraints = clist, trace = TRUE) Fit1 <- rcim(Qvar(fit1, "ethnicity", which.linpred = 1), uninormal("explink", imethod = 1), maxit = 99, trace = TRUE) Fit2 <- rcim(Qvar(fit1, "ethnicity", which.linpred = 2), uninormal("explink", imethod = 1), maxit = 99, trace = TRUE) ## End(Not run) ## Not run: par(mfrow = c(1, 2)) qvplot(Fit1, scol = "blue", pch = 16, main = expression(eta), slwd = 1.5, las = 1, length.arrows = 0.07) qvplot(Fit2, scol = "blue", pch = 16, main = expression(eta), slwd = 1.5, las = 1, length.arrows = 0.07) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.