| ## 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")) 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] ~ .)))) |
| |