confintvglm.RdComputes confidence intervals (CIs)
for one or more parameters in a fitted model.
Currently the object must be a
"vglm" object.
confintvglm(object, parm, level = 0.95, method = c("wald", "profile"),
trace = NULL, ...)A fitted model object.
Same as confint.
Character.
The default is the first method.
Abbreviations are allowed.
Currently "profile" is basically working;
and it is likely to be more accurate especially for
small samples, as it is based on a profile log likelihood,
however it is computationally intensive.
Logical. If TRUE then one can monitor the
computation as it progresses (because it is expensive).
The default is the orginal model's trace value
(see vglm.control).
Setting FALSE suppresses all intermediate output.
The default for
this methods function is based on confint.default
and assumes
asymptotic normality. In particular,
the coef and
vcov methods functions are used for
vglm-class objects.
When method = "profile" the function
profilevglm
is called to do the profiling. The code is very heavily
based on profile.glm
which was originally written by
D. M. Bates and W. N. Venables (For S in 1996)
and subsequently corrected by B. D. Ripley.
Sometimes the profiling method can give problems, for
example, cumulative requires the \(M\)
linear predictors not to intersect in the data cloud.
Such numerical problems are less common when
method = "wald", however, it is well-known
that inference based on profile likelihoods is generally
more accurate than Wald, especially when the sample size
is small.
The deviance (deviance(object)) is used if possible,
else the difference
2 * (logLik(object) - ell) is computed,
where ell are the values of the loglikelihood on a grid.
For
Wald CIs and
rrvglm-class
objects, currently an error message is produced because
I haven't gotten around to write the methods function;
it's not too hard, but am too busy!
An interim measure is to
coerce the object into a "vglm" object,
but then the confidence intervals will tend to be too narrow because
the estimated constraint matrices are treated as known.
For
Wald CIs and
vgam-class
objects, currently an error message is produced because
the theory is undeveloped.
Same as confint.
The order of the values of argument method may change
in the future without notice.
The functions
plot.profile.glm
and
pairs.profile.glm
from MASS
appear to work with output from this function.
vcovvlm,
summaryvglm,
confint,
profile.glm,
lrt.stat.vlm,
wald.stat,
plot.profile.glm,
pairs.profile.glm.
# Example 1: this is based on a glm example
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3, 1, 9); treatment <- gl(3, 3)
glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())
vglm.D93 <- vglm(counts ~ outcome + treatment, family = poissonff)
confint(glm.D93) # needs MASS to be present on the system
#> Waiting for profiling to be done...
#> 2.5 % 97.5 %
#> (Intercept) 2.6958215 3.36655581
#> outcome2 -0.8577018 -0.06255840
#> outcome3 -0.6753696 0.08244089
#> treatment2 -0.3932548 0.39325483
#> treatment3 -0.3932548 0.39325483
confint.default(glm.D93) # based on asymptotic normality
#> 2.5 % 97.5 %
#> (Intercept) 2.7095672 3.37947764
#> outcome2 -0.8505027 -0.05800787
#> outcome3 -0.6707552 0.08478093
#> treatment2 -0.3919928 0.39199279
#> treatment3 -0.3919928 0.39199279
confint(vglm.D93)
#> 2.5 % 97.5 %
#> (Intercept) 2.7095672 3.37947764
#> outcome2 -0.8505027 -0.05800787
#> outcome3 -0.6707552 0.08478093
#> treatment2 -0.3919928 0.39199279
#> treatment3 -0.3919928 0.39199279
confint(vglm.D93) - confint(glm.D93) # Should be all 0s
#> Waiting for profiling to be done...
#> 2.5 % 97.5 %
#> (Intercept) 0.013745734 0.012921826
#> outcome2 0.007199164 0.004550529
#> outcome3 0.004614421 0.002340033
#> treatment2 0.001262037 -0.001262037
#> treatment3 0.001262036 -0.001262036
confint(vglm.D93) - confint.default(glm.D93) # based on asympt. normality
#> 2.5 % 97.5 %
#> (Intercept) 9.862333e-12 -9.865886e-12
#> outcome2 3.233246e-10 -3.233231e-10
#> outcome3 1.765693e-10 -1.765685e-10
#> treatment2 2.162424e-10 -2.162407e-10
#> treatment3 1.785788e-10 -1.785767e-10
# Example 2: simulated negative binomial data with multiple responses
ndata <- data.frame(x2 = runif(nn <- 100))
ndata <- transform(ndata, y1 = rnbinom(nn, mu = exp(3+x2), size = exp(1)),
y2 = rnbinom(nn, mu = exp(2-x2), size = exp(0)))
fit1 <- vglm(cbind(y1, y2) ~ x2, negbinomial, data = ndata, trace = TRUE)
#> Iteration 1: loglikelihood = -686.41283
#> Iteration 2: loglikelihood = -681.75536
#> Iteration 3: loglikelihood = -681.61541
#> Iteration 4: loglikelihood = -681.61531
#> Iteration 5: loglikelihood = -681.61531
coef(fit1)
#> (Intercept):1 (Intercept):2 (Intercept):3 (Intercept):4 x2:1
#> 2.8798610 1.2293548 2.1593298 0.3344336 1.2198740
#> x2:2
#> -1.3628664
coef(fit1, matrix = TRUE)
#> loglink(mu1) loglink(size1) loglink(mu2) loglink(size2)
#> (Intercept) 2.879861 1.229355 2.159330 0.3344336
#> x2 1.219874 0.000000 -1.362866 0.0000000
confint(fit1)
#> 2.5 % 97.5 %
#> (Intercept):1 2.66075946 3.0989626
#> (Intercept):2 0.92725911 1.5314505
#> (Intercept):3 1.80408739 2.5145722
#> (Intercept):4 -0.04950863 0.7183759
#> x2:1 0.83295723 1.6067907
#> x2:2 -2.02944689 -0.6962860
confint(fit1, "x2:1") # This might be improved to "x2" some day...
#> 2.5 % 97.5 %
#> x2:1 0.8329572 1.606791
if (FALSE) { # \dontrun{
confint(fit1, method = "profile") # Computationally expensive
confint(fit1, "x2:1", method = "profile", trace = FALSE)
} # }
fit2 <- rrvglm(y1 ~ x2, negbinomial(zero = NULL), data = ndata)
confint(as(fit2, "vglm")) # Too narrow (SEs are biased downwards)
#> 2.5 % 97.5 %
#> (Intercept):1 2.6478171 3.102459
#> (Intercept):2 0.7183359 1.334768
#> x2 0.8471583 1.611676