1## tests of the simulate.lm method, added Feb 2009
2
3options(digits = 5)
4
5## cases should be named
6hills <- readRDS("hills.rds") # copied from package MASS
7fit1 <- lm(time ~ dist, data = hills)
8set.seed(1)
9simulate(fit1, nsim = 3)
10
11## and weights should be taken into account
12fit2 <- lm(time ~ -1 + dist + climb, hills[-18, ], weight = 1/dist^2)
13coef(summary(fit2))
14set.seed(1); ( ys <- simulate(fit2, nsim = 3) )
15for(i in seq_len(3))
16    print(coef(summary(update(fit2, ys[, i] ~ .))))
17## should be identical to glm(*, gaussian):
18fit2. <- glm(time ~ -1 + dist + climb, family=gaussian, data=hills[-18, ],
19             weight = 1/dist^2)
20set.seed(1); ys. <- simulate(fit2., nsim = 3)
21stopifnot(all.equal(ys, ys.))
22
23## Poisson fit
24load("anorexia.rda") # copied from package MASS
25fit3 <- glm(Postwt ~ Prewt + Treat + offset(Prewt),
26            family = gaussian, data = anorexia)
27coef(summary(fit3))
28set.seed(1)
29simulate(fit3, nsim = 8)
30
31## two-column binomial fit
32ldose <- rep(0:5, 2)
33numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16)
34sex <- factor(rep(c("M", "F"), c(6, 6)))
35SF <- cbind(numdead, numalive = 20 - numdead)
36fit4 <- glm(SF ~ sex + ldose - 1, family = binomial)
37coef(summary(fit4))
38set.seed(1)
39( ys <- simulate(fit4, nsim = 3) )
40for(i in seq_len(3))
41    print(coef(summary(update(fit4, ys[, i] ~ .))))
42
43## same via proportions
44fit5 <- glm(numdead/20 ~ sex + ldose - 1, family = binomial,
45            weights = rep(20, 12))
46set.seed(1)
47( ys <- simulate(fit5, nsim = 3) )
48for(i in seq_len(3))
49    print(coef(summary(update(fit5, ys[, i] ~ .))))
50
51
52## factor binomial fit
53load("birthwt.rda") # copied from package MASS
54bwt <- with(birthwt, {
55    race <- factor(race, labels = c("white", "black", "other"))
56    table(ptl)
57    ptd <- factor(ptl > 0)
58    table(ftv)
59    ftv <- factor(ftv)
60    levels(ftv)[-(1:2)] <- "2+"
61    data.frame(low = factor(low), age, lwt, race,
62               smoke = (smoke > 0), ptd, ht = (ht > 0), ui = (ui > 0), ftv)
63})
64fit6 <- glm(low ~ ., family = binomial, data = bwt)
65coef(summary(fit6))
66set.seed(1)
67ys <- simulate(fit6, nsim = 3)
68ys[1:10, ]
69for(i in seq_len(3))
70    print(coef(summary(update(fit6, ys[, i] ~ .))))
71
72## This requires MASS::gamma.shape
73if(!require("MASS", quietly = TRUE)) {
74    message("skipping tests requiring the MASS package")
75    q()
76}
77
78## gamma fit, from example(glm)
79clotting <- data.frame(u = c(5,10,15,20,30,40,60,80,100),
80                       lot1 = c(118,58,42,35,27,25,21,19,18))
81fit7 <- glm(lot1 ~ log(u), data = clotting, family = Gamma)
82coef(summary(fit7))
83set.seed(1)
84( ys <- simulate(fit7, nsim = 3) )
85for(i in seq_len(3))
86    print(coef(summary(update(fit7, ys[, i] ~ .))))
87
88