3
%- Also NEED an '\alias' for EACH other topic documented here.
5
Quasi-variances Preprocessing Function
7
%% ~~function to do ... ~~
10
Takes a \code{\link{vglm}} fit or a variance-covariance matrix,
11
and preprocesses it for \code{\link{rcam}} and
12
\code{\link{normal1}} so that quasi-variances can be computed.
15
%% ~~ A concise (1-5 lines) description of what the function does. ~~
18
Qvar(object, factorname = NULL, coef.indices = NULL,
19
labels = NULL, dispersion = NULL, reference.name = "(reference)",
22
%- maybe also 'usage' for other objects documented here.
25
A \code{"\link[=vglmff-class]{vglm}"} object or a variance-covariance
26
matrix, e.g., \code{vcov(vglm.object)}.
27
The former is preferred since it contains all the information
29
If a matrix then \code{factorname} and/or \code{coef.indices}
30
should be specified to identify the factor.
33
%% ~~Describe \code{object} here~~
37
If the \code{\link{vglm}} object contains more than one
38
factor as explanatory variable then this argument should
39
be the name of the factor of interest.
40
If \code{object} is a variance-covariance matrix then
41
this argument should also be specified.
44
%% ~~Describe \code{factor.name} here~~
48
Optional, for labelling the variance-covariance matrix.
50
%% ~~Describe \code{level1.name} here~~
54
Optional, passed into \code{vcov()} with the same argument name.
56
%% ~~Describe \code{level1.name} here~~
58
\item{reference.name}{
60
Label for for the reference level.
62
%% ~~Describe \code{level1.name} here~~
65
Optional numeric vector of length at least 3 specifying
66
the indices of the factor from the variance-covariance
72
an optional vector of estimated coefficients
73
(redundant if \code{object} is a model).
80
Suppose a factor with \eqn{L} levels is an explanatory variable in a
81
regression model. By default, R treats the first level as baseline so
82
that its coefficient is set to zero. It estimates the other \eqn{L-1}
83
coefficients, and with its associated standard errors, this is the
84
conventional output. From the complete variance-covariance matrix one
85
can compute \eqn{L} quasi-variances based on all pairwise difference
86
of the coefficients. They are based on an approximation, and can be
87
treated as uncorrelated. In minimizing the relative (not absolute)
88
errors it is not hard to see that the estimation involves a RCAM
89
(\code{\link{rcam}}) with an exponential link function
90
(\code{\link{explink}}).
93
If \code{object} is a model, then at least one of \code{factorname} or
94
\code{coef.indices} must be non-\code{NULL}. The value of
95
\code{coef.indices}, if non-\code{NULL}, determines which rows and
96
columns of the model's variance-covariance matrix to use. If
97
\code{coef.indices} contains a zero, an extra row and column are
98
included at the indicated position, to represent the zero variances
99
and covariances associated with a reference level. If
100
\code{coef.indices} is \code{NULL}, then \code{factorname} should be
101
the name of a factor effect in the model, and is used in order to
102
extract the necessary variance-covariance estimates.
105
Quasi-variances were first implemented in R with \pkg{qvcalc}.
106
This implementation draws heavily from that.
110
%% ~~ If necessary, more details than the description above ~~
113
A \eqn{L} by \eqn{L} matrix whose \eqn{i}-\eqn{j} element
114
is the logarithm of the variance of the \eqn{i}th coefficient
115
minus the \eqn{j}th coefficient, for all values of \eqn{i}
116
and \eqn{j}. The diagonal elements are abitrary and are set
120
The matrix has an attribute that corresponds to the prior
121
weight matrix; it is accessed by \code{\link{normal1}}
122
and replaces the usual \code{weights} argument.
123
of \code{\link{vglm}}. This weight matrix has ones on
124
the off-diagonals and some small positive number on
133
Overcoming the reference category problem in the
134
presentation of statistical models.
135
\emph{Sociological Methodology} \bold{33}, 1--18.
138
Firth, D. and Menezes, R. X. de (2004)
140
\emph{Biometrika} \bold{91}, 65--80.
148
T. W. Yee, based heavily on \code{qvcalc()} in \pkg{qvcalc}
149
written by David Firth.
155
This is an adaptation of \code{qvcalc()} in \pkg{qvcalc}.
156
It should work for all \code{\link{vglm}}
157
models with one linear predictor, i.e., \eqn{M = 1}.
158
For \eqn{M > 1} the factor should appear only in one of the
162
It is important to set \code{maxit} to be larger than usual for
163
\code{\link{rcam}} since convergence is slow. Upon successful
164
convergence the \eqn{i}th row effect and the \eqn{i}th column effect
165
should be equal. A simple computation involving the fitted and
166
predicted values allows the quasi-variances to be extracted (see
170
A function to plot \emph{comparison intervals} has not been
176
Negative quasi-variances may occur (one of them and
177
only one), though they are rare in practice. If
178
so then numerical problems may occur. See
179
\code{qvcalc()} for more information.
189
\code{\link{normal1}},
190
\code{\link{explink}},
191
\code{qvcalc()} in \pkg{qvcalc},
192
\code{\link[MASS]{ships}}.
195
%% ~~objects to See Also as \code{\link{help}}, ~~~
198
library(MASS) # Get the "ships" data frame
201
detach("package:MASS")
203
Shipmodel <- vglm(incidents ~ type + year + period,
204
quasipoissonff, offset = log(service),
205
# trace = TRUE, model = TRUE,
206
data = ships, subset = (service > 0))
208
# Easiest form of input
209
fit1 <- rcam(Qvar(Shipmodel, "type"), normal1("explink"), maxit = 99)
210
(quasiVar <- exp(diag(fitted(fit1))) / 2) # Version 1
211
(quasiVar <- diag(predict(fit1)[, c(TRUE, FALSE)]) / 2) # Version 2
212
(quasiSE <- sqrt(quasiVar))
214
# Another form of input
215
fit2 <- rcam(Qvar(Shipmodel, coef.ind = c(0,2:5), reference.name = "typeA"),
216
normal1("explink"), maxit = 99)
217
\dontrun{ plotqvar(fit2, lcol = "blue", llwd = 2, las = 1) }
219
# The variance-covariance matrix is another form of input (not recommended)
220
fit3 <- rcam(Qvar(cbind(0, rbind(0, vcov(Shipmodel)[2:5, 2:5])),
221
labels = c("typeA", "typeB", "typeC", "typeD", "typeE"),
222
estimates = c(typeA = 0, coef(Shipmodel)[2:5])),
223
normal1("explink"), maxit = 99)
224
(QuasiVar <- exp(diag(fitted(fit3))) / 2) # Version 1
225
(QuasiVar <- diag(predict(fit3)[, c(TRUE, FALSE)]) / 2) # Version 2
226
(QuasiSE <- sqrt(quasiVar))
227
\dontrun{ plotqvar(fit3) }
229
% Add one or more standard keywords, see file 'KEYWORDS' in the
230
% R documentation directory.
233
% \code{\link[qvcalc:qvcalc]{qvcalc}} in \pkg{qvcalc}
added
removed
man/Qvar.Rd
