Prior and Working Weights of a VGLM fit
Returns either the prior weights or working weights of a VGLM object.
weightsvglm(object, type = c("prior", "working"),
matrix.arg = TRUE, ignore.slot = FALSE,
deriv.arg = FALSE, ...)object |
a model object from the VGAM R package
that inherits from a vector generalized linear model (VGLM),
e.g., a model of class |
type |
Character, which type of weight is to be returned? The default is the first one. |
matrix.arg |
Logical, whether the answer is returned as a matrix. If not, it will be a vector. |
ignore.slot |
Logical. If |
deriv.arg |
Logical. If |
... |
Currently ignored. |
Working weights are used by the IRLS algorithm. They correspond to the second derivatives of the log-likelihood function with respect to the linear predictors. The working weights correspond to positive-definite weight matrices and are returned in matrix-band form, e.g., the first M columns correspond to the diagonals, etc.
If one wants to perturb the linear predictors then the
fitted.values slots should be assigned to the object
before calling this function. The reason is that,
for some family functions,
the variable mu is used directly as one of the parameter
estimates, without recomputing it from eta.
If type = "working" and deriv = TRUE then a
list is returned with the two components described below.
Otherwise the prior or working weights are returned depending
on the value of type.
deriv |
Typically the first derivative of the
log-likelihood with respect to the linear predictors.
For example, this is the variable |
weights |
The working weights. |
This function is intended to be similar to
weights.glm (see glm).
Thomas W. Yee
glm,
vglmff-class,
vglm.
pneumo <- transform(pneumo, let = log(exposure.time))
(fit <- vglm(cbind(normal, mild, severe) ~ let,
cumulative(parallel = TRUE, reverse = TRUE), data = pneumo))
depvar(fit) # These are sample proportions
weights(fit, type = "prior", matrix = FALSE) # Number of observations
# Look at the working residuals
nn <- nrow(model.matrix(fit, type = "lm"))
M <- ncol(predict(fit))
wwt <- weights(fit, type = "working", deriv = TRUE) # In matrix-band format
wz <- m2a(wwt$weights, M = M) # In array format
wzinv <- array(apply(wz, 3, solve), c(M, M, nn))
wresid <- matrix(NA, nn, M) # Working residuals
for (ii in 1:nn)
wresid[ii, ] <- wzinv[, , ii, drop = TRUE] %*% wwt$deriv[ii, ]
max(abs(c(resid(fit, type = "work")) - c(wresid))) # Should be 0
(zedd <- predict(fit) + wresid) # Adjusted dependent vectorPlease choose more modern alternatives, such as Google Chrome or Mozilla Firefox.