hw10data = read.table("./data/ex12-45.TXT", header = T)
aggLevel = hw10data[1:12, 1]
race = as.factor(hw10data[1:12, 2])
race = relevel(race, ref="White")
deathYes = hw10data[1:12, 4]
deathNo = hw10data[13:24, 4]
hw10dataNew = data.frame(aggLevel, race, deathYes, deathNo); hw10dataNew
## aggLevel race deathYes deathNo
## 1 1 White 2 60
## 2 1 Black 1 181
## 3 2 White 2 15
## 4 2 Black 1 21
## 5 3 White 6 7
## 6 3 Black 2 9
## 7 4 White 9 3
## 8 4 Black 2 4
## 9 5 White 9 0
## 10 5 Black 4 3
## 11 6 White 17 0
## 12 6 Black 4 0
hw10data = read.table("./data/ex12-45.TXT", header = T)
aggLevel = as.factor(hw10data[1:12, 1])
race = as.factor(hw10data[1:12, 2])
race = relevel(race, ref="White")
deathYes = hw10data[1:12, 4]
deathNo = hw10data[13:24, 4]
mod = glm(cbind(deathYes, deathNo)~1+aggLevel+race, family = binomial(link = "logit"))
summary(mod)
##
## Call:
## glm(formula = cbind(deathYes, deathNo) ~ 1 + aggLevel + race,
## family = binomial(link = "logit"))
##
## Deviance Residuals:
## 1 2 3 4 5 6 7
## 0.02705 -0.03705 -0.27695 0.46062 -0.22255 0.33222 0.02846
## 8 9 10 11 12
## -0.03695 1.21437 -0.55797 0.00006 0.00007
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.4207 0.6144 -5.567 2.59e-08 ***
## aggLevel2 1.6090 0.8506 1.892 0.05855 .
## aggLevel3 3.3902 0.7474 4.536 5.74e-06 ***
## aggLevel4 4.5004 0.7858 5.727 1.02e-08 ***
## aggLevel5 5.8814 0.9128 6.443 1.17e-10 ***
## aggLevel6 26.2636 8772.8073 0.003 0.99761
## raceBlack -1.7409 0.5426 -3.208 0.00134 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 212.2838 on 11 degrees of freedom
## Residual deviance: 2.2391 on 5 degrees of freedom
## AIC: 38.105
##
## Number of Fisher Scoring iterations: 19
Set aggravation level 1 as the reference level, then the odds ratio for the other five levels are:
hw10data = read.table("./data/ex12-45.TXT", header = T)
aggLevel = hw10data[1:12, 1]
race = as.factor(hw10data[1:12, 2])
race = relevel(race, ref="White")
deathYes = hw10data[1:12, 4]
deathNo = hw10data[13:24, 4]
mod = glm(cbind(deathYes, deathNo)~1+aggLevel+race, family = binomial(link = "logit"))
summary(mod)
##
## Call:
## glm(formula = cbind(deathYes, deathNo) ~ 1 + aggLevel + race,
## family = binomial(link = "logit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.93570 -0.22548 0.05142 0.65620 1.01444
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.8653 0.6004 -8.103 5.37e-16 ***
## aggLevel 1.5397 0.1867 8.246 < 2e-16 ***
## raceBlack -1.8106 0.5361 -3.377 0.000732 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 212.2838 on 11 degrees of freedom
## Residual deviance: 3.8816 on 9 degrees of freedom
## AIC: 31.747
##
## Number of Fisher Scoring iterations: 4
beta0 = coef(mod)['(Intercept)']
betaAggLevel = coef(mod)['aggLevel']
betaRace = coef(mod)['raceBlack']
Model: \(logit(\hat{\pi})\) = -4.8653284 + 1.5396614*(aggLevel) -1.8106466*(race=black)
pVal = coef(summary(mod))['aggLevel', 'Pr(>|z|)']; pVal
## [1] 1.64349e-16
Yes. p-value = 1.64349e-16, indicating a strong evidence of association between the severity of the crime and the probability of receiving the death penalty.
dataBlack = hw10dataNew[race=="Black", ]
modBlack = glm(cbind(deathYes, deathNo)~aggLevel, family = binomial(link = "logit"), data=dataBlack)
summary(modBlack)
##
## Call:
## glm(formula = cbind(deathYes, deathNo) ~ aggLevel, family = binomial(link = "logit"),
## data = dataBlack)
##
## Deviance Residuals:
## 2 4 6 8 10 12
## -0.3966 0.3456 0.6146 -0.1022 -0.6614 0.9132
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.2319 0.9015 -6.913 4.75e-12 ***
## aggLevel 1.4067 0.2466 5.705 1.17e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 57.5843 on 5 degrees of freedom
## Residual deviance: 1.9364 on 4 degrees of freedom
## AIC: 16.973
##
## Number of Fisher Scoring iterations: 4
dataWhite = hw10dataNew[race=="White", ]
modWhite = glm(cbind(deathYes, deathNo)~aggLevel, family = binomial(link = "logit"), data=dataWhite)
summary(modWhite)
##
## Call:
## glm(formula = cbind(deathYes, deathNo) ~ aggLevel, family = binomial(link = "logit"),
## data = dataWhite)
##
## Deviance Residuals:
## 1 3 5 7 9 11
## 0.2315 -0.1673 0.1000 -0.5411 0.8681 0.5192
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.2531 0.8480 -6.194 5.85e-10 ***
## aggLevel 1.6811 0.2842 5.916 3.30e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 106.2819 on 5 degrees of freedom
## Residual deviance: 1.4074 on 4 degrees of freedom
## AIC: 16.236
##
## Number of Fisher Scoring iterations: 4
library(HH)
#black_y = predict(modBlack, data.frame(aggLevel=3), type="response"); black_y
black_CI = interval(modBlack, newdata=data.frame(aggLevel=3), type="response"); black_CI
## fit ci.low ci.hi pi.low pi.hi
## 1 0.1179736 0.04227748 0.2883896 0.006684025 0.7266728
# white_y = predict(modWhite, data.frame(aggLevel=3), type="response"); white_y
white_CI = interval(modWhite, newdata=data.frame(aggLevel=3), type="response"); white_CI
## fit ci.low ci.hi pi.low pi.hi
## 1 0.4477287 0.240278 0.6751238 0.04142656 0.938302