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)