88 lines
2.6 KiB
R
88 lines
2.6 KiB
R
## tests of the simulate.lm method, added Feb 2009
|
|
|
|
options(digits = 5)
|
|
|
|
## cases should be named
|
|
hills <- readRDS("hills.rds") # copied from package MASS
|
|
fit1 <- lm(time ~ dist, data = hills)
|
|
set.seed(1)
|
|
simulate(fit1, nsim = 3)
|
|
|
|
## and weights should be taken into account
|
|
fit2 <- lm(time ~ -1 + dist + climb, hills[-18, ], weight = 1/dist^2)
|
|
coef(summary(fit2))
|
|
set.seed(1); ( ys <- simulate(fit2, nsim = 3) )
|
|
for(i in seq_len(3))
|
|
print(coef(summary(update(fit2, ys[, i] ~ .))))
|
|
## should be identical to glm(*, gaussian):
|
|
fit2. <- glm(time ~ -1 + dist + climb, family=gaussian, data=hills[-18, ],
|
|
weight = 1/dist^2)
|
|
set.seed(1); ys. <- simulate(fit2., nsim = 3)
|
|
stopifnot(all.equal(ys, ys.))
|
|
|
|
## Poisson fit
|
|
load("anorexia.rda") # copied from package MASS
|
|
fit3 <- glm(Postwt ~ Prewt + Treat + offset(Prewt),
|
|
family = gaussian, data = anorexia)
|
|
coef(summary(fit3))
|
|
set.seed(1)
|
|
simulate(fit3, nsim = 8)
|
|
|
|
## two-column binomial fit
|
|
ldose <- rep(0:5, 2)
|
|
numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16)
|
|
sex <- factor(rep(c("M", "F"), c(6, 6)))
|
|
SF <- cbind(numdead, numalive = 20 - numdead)
|
|
fit4 <- glm(SF ~ sex + ldose - 1, family = binomial)
|
|
coef(summary(fit4))
|
|
set.seed(1)
|
|
( ys <- simulate(fit4, nsim = 3) )
|
|
for(i in seq_len(3))
|
|
print(coef(summary(update(fit4, ys[, i] ~ .))))
|
|
|
|
## same via proportions
|
|
fit5 <- glm(numdead/20 ~ sex + ldose - 1, family = binomial,
|
|
weights = rep(20, 12))
|
|
set.seed(1)
|
|
( ys <- simulate(fit5, nsim = 3) )
|
|
for(i in seq_len(3))
|
|
print(coef(summary(update(fit5, ys[, i] ~ .))))
|
|
|
|
|
|
## factor binomial fit
|
|
load("birthwt.rda") # copied from package MASS
|
|
bwt <- with(birthwt, {
|
|
race <- factor(race, labels = c("white", "black", "other"))
|
|
table(ptl)
|
|
ptd <- factor(ptl > 0)
|
|
table(ftv)
|
|
ftv <- factor(ftv)
|
|
levels(ftv)[-(1:2)] <- "2+"
|
|
data.frame(low = factor(low), age, lwt, race,
|
|
smoke = (smoke > 0), ptd, ht = (ht > 0), ui = (ui > 0), ftv)
|
|
})
|
|
fit6 <- glm(low ~ ., family = binomial, data = bwt)
|
|
coef(summary(fit6))
|
|
set.seed(1)
|
|
ys <- simulate(fit6, nsim = 3)
|
|
ys[1:10, ]
|
|
for(i in seq_len(3))
|
|
print(coef(summary(update(fit6, ys[, i] ~ .))))
|
|
|
|
## This requires MASS::gamma.shape
|
|
if(!require("MASS", quietly = TRUE)) {
|
|
message("skipping tests requiring the MASS package")
|
|
q()
|
|
}
|
|
|
|
## gamma fit, from example(glm)
|
|
clotting <- data.frame(u = c(5,10,15,20,30,40,60,80,100),
|
|
lot1 = c(118,58,42,35,27,25,21,19,18))
|
|
fit7 <- glm(lot1 ~ log(u), data = clotting, family = Gamma)
|
|
coef(summary(fit7))
|
|
set.seed(1)
|
|
( ys <- simulate(fit7, nsim = 3) )
|
|
for(i in seq_len(3))
|
|
print(coef(summary(update(fit7, ys[, i] ~ .))))
|
|
|