-
Notifications
You must be signed in to change notification settings - Fork 0
/
RinA CH13 Code.txt
executable file
·96 lines (78 loc) · 3.02 KB
/
RinA CH13 Code.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
#--------------------------------------------------#
# R in Action: Chapter 13 #
# requires that the AER, robust, and qcc packages #
# have been installed #
# install.packages(c('AER', 'robust', 'qcc')) #
#--------------------------------------------------#
# --Logistic Regression--
# get summary statistics
data(Affairs, package = "AER")
summary(Affairs)
table(Affairs$affairs)
# create binary outcome variable
Affairs$ynaffair[Affairs$affairs > 0] <- 1
Affairs$ynaffair[Affairs$affairs == 0] <- 0
Affairs$ynaffair <- factor(Affairs$ynaffair, levels = c(0,
1), labels = c("No", "Yes"))
table(Affairs$ynaffair)
# fit full model
fit.full <- glm(ynaffair ~ gender + age + yearsmarried +
children + religiousness + education + occupation + rating,
data = Affairs, family = binomial())
summary(fit.full)
# fit reduced model
fit.reduced <- glm(ynaffair ~ age + yearsmarried +
religiousness + rating, data = Affairs, family = binomial())
summary(fit.reduced)
# compare models
anova(fit.reduced, fit.full, test = "Chisq")
# interpret coefficients
coef(fit.reduced)
exp(coef(fit.reduced))
# calculate probability of extramariatal affair by marital ratings
testdata <- data.frame(rating = c(1, 2, 3, 4, 5),
age = mean(Affairs$age), yearsmarried = mean(Affairs$yearsmarried),
religiousness = mean(Affairs$religiousness))
testdata$prob <- predict(fit.reduced, newdata = testdata,
type = "response")
testdata
# calculate probabilites of extramariatal affair by age
testdata <- data.frame(rating = mean(Affairs$rating),
age = seq(17, 57, 10), yearsmarried = mean(Affairs$yearsmarried),
religiousness = mean(Affairs$religiousness))
testdata$prob <- predict(fit.reduced, newdata = testdata,
type = "response")
testdata
# evaluate overdispersion
fit <- glm(ynaffair ~ age + yearsmarried + religiousness +
rating, family = binomial(), data = Affairs)
fit.od <- glm(ynaffair ~ age + yearsmarried + religiousness +
rating, family = quasibinomial(), data = Affairs)
pchisq(summary(fit.od)$dispersion * fit$df.residual,
fit$df.residual, lower = F)
# --Poisson Regression--
# look at dataset
data(breslow.dat, package = "robust")
names(breslow.dat)
summary(breslow.dat[c(6, 7, 8, 10)])
# plot distribution of post-treatment seizure counts
opar <- par(no.readonly = TRUE)
par(mfrow = c(1, 2))
attach(breslow.dat)
hist(sumY, breaks = 20, xlab = "Seizure Count", main = "Distribution of Seizures")
boxplot(sumY ~ Trt, xlab = "Treatment", main = "Group Comparisons")
par(opar)
# fit regression
fit <- glm(sumY ~ Base + Age + Trt, data = breslow.dat,
family = poisson())
summary(fit)
# interpret model parameters
coef(fit)
exp(coef(fit))
# evaluate overdispersion
library(qcc)
qcc.overdispersion.test(breslow.dat$sumY, type = "poisson")
# fit model with quasipoisson
fit.od <- glm(sumY ~ Base + Age + Trt, data = breslow.dat,
family = quasipoisson())
summary(fit.od)