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