The R Project SVN R

Rev

Rev 48011 | Rev 69217 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

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