## Overdispersed Poisson loglinear model for ship damage data
## from McCullagh and Nelder (1989), Sec 6.3.2
if (require(MASS)) {
data(ships)
ships$year <- as.factor(ships$year)
ships$period <- as.factor(ships$period)
shipmodel <- glm(formula = incidents ~ type + year + period,
family = quasipoisson,
data = ships,
subset = (service > 0),
offset = log(service))
shiptype.qv <- qvcalc(shipmodel, "type")
## We can plot "comparison intervals" as follows:
## plot(shiptype.qv, xlab = "ship type")
## An equivalent result by using the coef.indices argument instead:
## shiptype.qv2 <- qvcalc(shipmodel, coef.indices = c(0, 2:5))
summary(shiptype.qv, digits = 4)
}
## Example of a "coxph" model
if(require(survival)) {
data("veteran", package = "survival")
cancer_model <- coxph(Surv(time,status) ~ celltype, data = veteran)
celltype_qv <- qvcalc(cancer_model, "celltype")
summary(celltype_qv)
}
## Example of a "survreg" model
if(require(survival)) {
data("veteran", package = "survival")
cancer_model2 <- survreg(Surv(time,status) ~ celltype, data = veteran,
dist = "weibull")
celltype_qv2 <- qvcalc(cancer_model2, "celltype")
summary(celltype_qv2)
}
## Based on an example from ?itempar
if(require(psychotools)) {
data("VerbalAggression", package = "psychotools")
raschmod <- raschmodel(VerbalAggression$resp2)
ip1 <- itempar(raschmod)
qv1 <- qvcalc(ip1)
summary(qv1) }
## Example of a negative quasi variance
## Requires the "car" package
if (FALSE) {
library(car)
data(Prestige)
attach(Prestige)
mymodel <- lm(prestige ~ type + education)
library(qvcalc)
type.qvs <- qvcalc(mymodel, "type")
## Warning message:
## In sqrt(qv) : NaNs produced
summary(type.qvs)
## Model call: lm(formula = prestige ~ type + education)
## Factor name: type
## estimate SE quasiSE quasiVar
## bc 0.000000 0.000000 2.874361 8.261952
## prof 6.142444 4.258961 3.142737 9.876793
## wc -5.458495 2.690667 NaN -1.022262
## Worst relative errors in SEs of simple contrasts (%): 0 0
## Worst relative errors over *all* contrasts (%): 0 0
plot(type.qvs)
## Error in plot.qv(type.qvs) : No comparison intervals available,
## since one of the quasi variances is negative. See ?qvcalc for more.
}
Run the code above in your browser using DataLab