The R program for chapter 17.
We will only be using the depress data set so we will attach it.
attach(depress)
Creating the dichotomous variable inc.high for high income (income > 19).
depress$inc.high <- 1*(income > 19)
Creating the factor variables for sex and inc.high.
depress$sex.fac <- factor(sex, labels=c("male", "female"))
depress$inc.high.fac <- factor(inc.high, labels=c("low", "high"))
Table 17.1, p. 447.
Classification of individuals by income and sex.
table(depress$sex.fac, depress$inc.high.fac)
low high
male 54 57
female 125 58
Creating the factor variable for treat.
depress$treat.fac <- factor(treat, labels=c("No", "Yes"))
Table 17.2, p. 448.
table(depress$sex.fac, depress$inc.high.fac, depress$treat.fac)
, , No
low high
male 20 21
female 73 34
, , Yes
low high
male 34 36
female 52 24
Creating the factor variable for cases.
depress$cases.fac <- factor(cases, labels=c("low", "high"))
Table 17.3, p. 449.
table(depress$sex.fac, depress$inc.high.fac, depress$treat.fac,
depress$cases.fac)
, , No, low
low high
male 17 18
female 57 26
, , Yes, low
low high
male 31 35
female 38 22
, , No, high
low high
male 3 3
female 16 8
, , Yes, high
low high
male 3 1
female 14 2
Output for chi-squared test, bottom p. 451.
chisq.test(depress$inc.high, depress$sex, correct=F) Pearson's chi-square test without Yates' continuity correction data: depress$inc.high and depress$sex X-square = 11.2104, df = 1, p-value = 0.0008
Table 17.7, p. 455.
table1 <- table(depress$sex, depress$inc.high)
loglin(table1, list(1:2), param=T, fit=T)
$fit:
0 1
1 54 57
2 125 58
$param:
$constant:
[1] 4.230198
$"1":
1 2
-0.2141804 0.2141804
$"2":
0 1
0.1784509 -0.1784509
$"1 X 2":
0 1
1 -0.2054845 0.2054845
2 0.2054845 -0.2054845
Output for three-way model, p. 461.
table3 <- table(depress$sex, depress$inc.high, depress$treat)
loglin(table3, list(1:3), param=T, fit=T)
$fit:
, , 1
0 1
1 20 21
2 73 34
, , 2
0 1
1 34 36
2 52 24
$param:
$constant:
[1] 3.512031
$"1":
1 2
-0.2244979 0.2244979
$"2":
0 1
0.1789175 -0.1789175
$"1 X 2":
0 1
1 -0.2054047 0.2054047
2 0.2054047 -0.2054047
$"3":
1 2
-0.04776275 0.04776275
$"1 X 3":
1 2
1 -0.2196434 0.2196434
2 0.2196434 -0.2196434
$"2 X 3":
1 2
0 -0.00009030104 0.00009030104
1 0.00009030104 -0.00009030104
$"1 X 2 X 3":
, , 1
0 1
1 0.002182364 -0.002182364
2 -0.002182364 0.002182364
, , 2
0 1
1 -0.002182364 0.002182364
2 0.002182364 -0.002182364
Calculations, top p. 470.
Note: We need to be modeling variables that are coded (0, 1) rather than (1, 2). Thus, we want the change both treat and predictor sex to be (0, 1) rather than (1, 2).
depress$sex1 <- depress$sex - 1 depress$treat1 <- depress$treat -1 glm(treat1 ~ sex1 + inc.high, depress, family=binomial) Coefficients: (Intercept) sex1 inc.high 0.5359694 -0.8774062 -0.002075705
Unless you will continue to work with the depress data then it might be a good idea to detach the data set.
detach(depress)