Bai tap ngay 4
Vu Thien
2022-12-26
1: Đọc dữ liệu “Arrest dataset.csv” vào R dùng hàm read.csv.
arr = read.csv("~/Desktop/R studying/Arrest data.csv")2: Đánh giá đặc điểm các biến số theo tình trạng bị bắt lại “arrest” theo bảng sau:
library(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-table1(~ age + race + educ + prior + work + married + parole + finance | factor(arrest), data = arr)| 0 (N=318) |
1 (N=114) |
Overall (N=432) |
|
|---|---|---|---|
| age | |||
| Mean (SD) | 25.3 (6.31) | 22.8 (5.12) | 24.6 (6.11) |
| Median [Min, Max] | 23.0 [17.0, 44.0] | 21.0 [17.0, 44.0] | 23.0 [17.0, 44.0] |
| race | |||
| black | 277 (87.1%) | 102 (89.5%) | 379 (87.7%) |
| other | 41 (12.9%) | 12 (10.5%) | 53 (12.3%) |
| educ | |||
| Mean (SD) | 3.53 (0.883) | 3.32 (0.656) | 3.48 (0.834) |
| Median [Min, Max] | 3.00 [2.00, 6.00] | 3.00 [2.00, 6.00] | 3.00 [2.00, 6.00] |
| prior | |||
| Mean (SD) | 2.70 (2.55) | 3.77 (3.59) | 2.98 (2.90) |
| Median [Min, Max] | 2.00 [0, 15.0] | 3.00 [0, 18.0] | 2.00 [0, 18.0] |
| work | |||
| no | 123 (38.7%) | 62 (54.4%) | 185 (42.8%) |
| yes | 195 (61.3%) | 52 (45.6%) | 247 (57.2%) |
| married | |||
| married | 45 (14.2%) | 8 (7.0%) | 53 (12.3%) |
| not married | 273 (85.8%) | 106 (93.0%) | 379 (87.7%) |
| parole | |||
| no | 119 (37.4%) | 46 (40.4%) | 165 (38.2%) |
| yes | 199 (62.6%) | 68 (59.6%) | 267 (61.8%) |
| finance | |||
| no | 150 (47.2%) | 66 (57.9%) | 216 (50.0%) |
| yes | 168 (52.8%) | 48 (42.1%) | 216 (50.0%) |
library(compareGroups)
createTable(compareGroups(arrest ~ age + race + educ + prior + work + married + parole + finance, data = arr))##
## --------Summary descriptives table by 'arrest'---------
##
## _________________________________________________
## 0 1 p.overall
## N=318 N=114
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## age 25.3 (6.31) 22.8 (5.12) <0.001
## race: 0.621
## black 277 (87.1%) 102 (89.5%)
## other 41 (12.9%) 12 (10.5%)
## educ 3.53 (0.88) 3.32 (0.66) 0.006
## prior 2.70 (2.55) 3.77 (3.59) 0.004
## work: 0.005
## no 123 (38.7%) 62 (54.4%)
## yes 195 (61.3%) 52 (45.6%)
## married: 0.068
## married 45 (14.2%) 8 (7.02%)
## not married 273 (85.8%) 106 (93.0%)
## parole: 0.660
## no 119 (37.4%) 46 (40.4%)
## yes 199 (62.6%) 68 (59.6%)
## finance: 0.063
## no 150 (47.2%) 66 (57.9%)
## yes 168 (52.8%) 48 (42.1%)
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯3: Đánh giá mối liên quan giữa tuổi “age” và khả năng bị bắt lại “arrest”
3.1 Vẽ biểu đồ hộp thể hiện mối liên quan giữa tuổi và khả năng bị bắt lại
library(ggplot2)
ggplot(arr, aes(group = arrest, x = arrest, y = age, fill = arrest, color = arrest)) + geom_boxplot(colour = "black") + geom_jitter(aes(color = arrest))+ theme(legend.position = "none")+ labs(x = "Arrest", y = "Age (years)")
## 3.2 Sử dụng Student’s t-test để kiểm tra giả thuyết là không có mối
liên quan giữa tuổi và khả năng bị bắt lại. Bạn nhận xét kết quả như thế
nào?
t.test(age ~ arrest, data = arr)##
## Welch Two Sample t-test
##
## data: age by arrest
## t = 4.1789, df = 243.6, p-value = 4.086e-05
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
## 1.317143 3.665975
## sample estimates:
## mean in group 0 mean in group 1
## 25.25472 22.763163.3. Dùng mô hình hồi qui logistic để đánh giá mối liên quan giữa tuổi và khả năng bị bắt lại. Trình bày kết quả theo bảng dưới đây:
m1 = glm(arrest ~ age, family = binomial, data = arr)
summary(m1)##
## Call:
## glm(formula = arrest ~ age, family = binomial, data = arr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9909 -0.8362 -0.7246 1.4131 2.3509
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.95611 0.54357 1.759 0.078586 .
## age -0.08305 0.02296 -3.617 0.000298 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 498.60 on 431 degrees of freedom
## Residual deviance: 482.69 on 430 degrees of freedom
## AIC: 486.69
##
## Number of Fisher Scoring iterations: 4library(epiDisplay)## Loading required package: foreign## Loading required package: survival## Loading required package: MASS## Loading required package: nnet##
## Attaching package: 'epiDisplay'## The following object is masked from 'package:ggplot2':
##
## alphalogistic.display(m1)##
## Logistic regression predicting arrest
##
## OR(95%CI) P(Wald's test) P(LR-test)
## age (cont. var.) 0.92 (0.88,0.96) < 0.001 < 0.001
##
## Log-likelihood = -241.3436
## No. of observations = 432
## AIC value = 486.68724: Đánh giá mối liên quan giữa hỗ trợ tài chánh “finance” và khả năng bị bắt lại “arrest”
library(ggplot2)
library(ggthemes)
ggplot(data = arr, aes(x = age, y = arrest)) + geom_point(alpha = 0.15) + geom_smooth(method = "glm", method.args = list(family = "binomial")) ## `geom_smooth()` using formula 'y ~ x'
## 4.1 Dùng mô hình hồi qui logistic đánh giá mối liên quan giữa hỗ trợ
tài chánh và khả năng bị bắt lại. Tóm tắt kết quả theo bảng dưới
đây:
library(rms)## Loading required package: Hmisc## Loading required package: lattice##
## Attaching package: 'lattice'## The following object is masked from 'package:epiDisplay':
##
## dotplot## Loading required package: Formula##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:table1':
##
## label, label<-, units## The following objects are masked from 'package:base':
##
## format.pval, units## Loading required package: SparseM##
## Attaching package: 'SparseM'## The following object is masked from 'package:base':
##
## backsolve## Registered S3 method overwritten by 'rms':
## method from
## print.lrtest epiDisplay##
## Attaching package: 'rms'## The following object is masked from 'package:epiDisplay':
##
## lrtestm2 = glm(arrest ~ finance, family = binomial, data = arr)
summary(m2)##
## Call:
## glm(formula = arrest ~ finance, family = binomial, data = arr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.854 -0.854 -0.709 1.540 1.734
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8210 0.1477 -5.558 2.73e-08 ***
## financeyes -0.4318 0.2205 -1.959 0.0502 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 498.60 on 431 degrees of freedom
## Residual deviance: 494.73 on 430 degrees of freedom
## AIC: 498.73
##
## Number of Fisher Scoring iterations: 4logistic.display(m2)##
## Logistic regression predicting arrest
##
## OR(95%CI) P(Wald's test) P(LR-test)
## finance (cont. var.) 0.65 (0.42,1) 0.05 0.049
##
## Log-likelihood = -247.3642
## No. of observations = 432
## AIC value = 498.72834.2 Vẽ biểu đồ hộp và dùng Student’s t-test để so sánh khác biệt giữa tuổi và khả năng hỗ trợ tài chánh. Bạn nhận xét gì về mối liên hệ giữa tuổi và khả năng hỗ trợ tài chánh.
# Biểu đồ hộp
ggplot(arr, aes(group = finance, x = finance, y = age, fill = finance, color = finance)) + geom_boxplot(colour = "black") + geom_jitter(aes(color = arrest))+theme(legend.position = "none")+ labs(x = "Finance", y = "Age (years)")
t.test(age ~ finance, data = arr)##
## Welch Two Sample t-test
##
## data: age by finance
## t = -1.2759, df = 423.8, p-value = 0.2027
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
## -1.9054289 0.4054289
## sample estimates:
## mean in group no mean in group yes
## 24.22222 24.972224.3 Đánh giá ảnh hưởng của hỗ trợ tài chánh lên khả năng bị bắt lại sau khi hiệu chỉnh cho khác biệt về độ tuổi theo bảng sau:
m3 = glm(arrest ~ finance + age, family = binomial, data = arr)
summary(m3)##
## Call:
## glm(formula = arrest ~ finance + age, family = binomial, data = arr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0633 -0.8193 -0.7107 1.3320 2.2949
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.11913 0.55462 2.018 0.043608 *
## financeyes -0.39905 0.22408 -1.781 0.074935 .
## age -0.08197 0.02309 -3.550 0.000385 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 498.60 on 431 degrees of freedom
## Residual deviance: 479.49 on 429 degrees of freedom
## AIC: 485.49
##
## Number of Fisher Scoring iterations: 4library(Publish)## Loading required package: prodlimpublish(m3)## Variable Units OddsRatio CI.95 p-value
## finance no Ref
## yes 0.67 [0.43;1.04] 0.0749350
## age 0.92 [0.88;0.96] 0.0003852Bạn có nhận xét gì về kết quả này? Kết quả này có gì khác biệt so với kết quả 4.1?
5: Xây dựng mô hình tối ưu dự báo khả năng bị bắt lại
library(BMA)## Loading required package: leaps## Loading required package: robustbase##
## Attaching package: 'robustbase'## The following object is masked from 'package:survival':
##
## heart## Loading required package: inline## Loading required package: rrcov## Scalable Robust Estimators with High Breakdown Point (version 1.7-2)yvar = arr[, "arrest"]
xvars = arr[, c("finance", "age", "race", "work", "married", "parole", "prior", "educ")]
bma = bic.glm(xvars, yvar, glm.family = "binomial")
summary(bma)##
## Call:
## bic.glm.data.frame(x = xvars, y = yvar, glm.family = "binomial")
##
##
## 12 models were selected
## Best 5 models (cumulative posterior probability = 0.7985 ):
##
## p!=0 EV SD model 1 model 2 model 3
## Intercept 100 0.817350 0.81303 5.022e-01 9.561e-01 1.466e+00
## finance 14.6 -0.058569 0.16599 . . .
## age 100.0 -0.078012 0.02328 -7.796e-02 -8.305e-02 -7.692e-02
## race 3.1 -0.010614 0.08693 . . .
## work 4.9 -0.014130 0.08348 . . .
## married 5.0 0.027326 0.15084 . . .
## parole 2.1 -0.002227 0.03709 . . .
## prior 75.0 0.077046 0.05457 1.046e-01 . 9.433e-02
## educ 20.8 -0.062332 0.14083 . . -2.802e-01
##
## nVar 2 1 3
## BIC -2.129e+03 -2.127e+03 -2.126e+03
## post prob 0.395 0.122 0.107
## model 4 model 5
## Intercept 6.565e-01 2.054e+00
## finance -4.057e-01 .
## age -7.648e-02 -8.108e-02
## race . .
## work . .
## married . .
## parole . .
## prior 1.056e-01 .
## educ . -3.351e-01
##
## nVar 3 2
## BIC -2.126e+03 -2.126e+03
## post prob 0.095 0.079imageplot.bma(bma)
6: Đánh giá mô hình tiên lượng
6.1 Sử dụng package “caret” để đánh giá mô hình tiên lượng
library(caret)##
## Attaching package: 'caret'## The following object is masked from 'package:survival':
##
## clustersample = createDataPartition(arr$arrest, p = 0.6, list = FALSE)
train = arr[sample,]
test = arr[-sample,]
m4 = train(factor(arrest) ~ age + prior, data = train, method = "glm", family = "binomial")
summary(m4)##
## Call:
## NULL
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2572 -0.8495 -0.7394 1.3305 2.2043
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.48130 0.65935 0.730 0.46541
## age -0.06767 0.02599 -2.604 0.00923 **
## prior 0.07954 0.04256 1.869 0.06163 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 312.40 on 259 degrees of freedom
## Residual deviance: 299.78 on 257 degrees of freedom
## AIC: 305.78
##
## Number of Fisher Scoring iterations: 4arrest.pred = predict(m4, newdata = test, type = "raw")
confusionMatrix(data = arrest.pred, reference = factor(test$arrest))## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 133 38
## 1 0 1
##
## Accuracy : 0.7791
## 95% CI : (0.7096, 0.8387)
## No Information Rate : 0.7733
## P-Value [Acc > NIR] : 0.4703
##
## Kappa : 0.0391
##
## Mcnemar's Test P-Value : 1.947e-09
##
## Sensitivity : 1.00000
## Specificity : 0.02564
## Pos Pred Value : 0.77778
## Neg Pred Value : 1.00000
## Prevalence : 0.77326
## Detection Rate : 0.77326
## Detection Prevalence : 0.99419
## Balanced Accuracy : 0.51282
##
## 'Positive' Class : 0
## 6.2 Sử dụng package “rms” để đánh giá mô hình tiên lượng.
library(rms)
m5 = lrm(arrest ~ age + prior, data = arr, x = TRUE, y = TRUE)
validate = validate(m5, B = 100)
validate## index.orig training test optimism index.corrected n
## Dxy 0.3132 0.3179 0.3078 0.0100 0.3032 100
## R2 0.0800 0.0894 0.0764 0.0130 0.0670 100
## Intercept 0.0000 0.0000 -0.0423 0.0423 -0.0423 100
## Slope 1.0000 1.0000 0.9619 0.0381 0.9619 100
## Emax 0.0000 0.0000 0.0161 0.0161 0.0161 100
## D 0.0540 0.0613 0.0514 0.0099 0.0441 100
## U -0.0046 -0.0046 0.0009 -0.0055 0.0009 100
## Q 0.0586 0.0659 0.0505 0.0154 0.0432 100
## B 0.1825 0.1820 0.1839 -0.0019 0.1844 100
## g 0.6416 0.6747 0.6195 0.0552 0.5864 100
## gp 0.1135 0.1175 0.1098 0.0076 0.1058 100Bạn có nhận xét gì về hai kết quả này.
7: Tính AUC cho mô hình hồi qui logistic
prob = predict(m4, newdata = test, type = "prob")
pred = data.frame(prob, test$arrest)
head(pred)## X0 X1 test.arrest
## 2 0.5250946 0.4749054 1
## 5 0.6378052 0.3621948 0
## 6 0.7278530 0.2721470 0
## 7 0.7703901 0.2296099 1
## 9 0.6295389 0.3704611 0
## 11 0.7387649 0.2612351 0library(pROC)## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following object is masked from 'package:epiDisplay':
##
## ci## The following objects are masked from 'package:stats':
##
## cov, smooth, varroc(pred$test.arrest, pred$X1)## Setting levels: control = 0, case = 1## Setting direction: controls < cases##
## Call:
## roc.default(response = pred$test.arrest, predictor = pred$X1)
##
## Data: pred$X1 in 133 controls (pred$test.arrest 0) < 39 cases (pred$test.arrest 1).
## Area under the curve: 0.6827plot(smooth(roc(pred$test.arrest, pred$X1)))## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Dùng package “pROC” tính AUC và vẽ đường cong ROC cho mô hình tiên
lượng khả năng bị bắt lại.
8: Xây dựng nomogram để dự báo khả năng bị bắt lại. Diễn giải cách sử dụng nomogram này.
ddist = datadist(arr)
options(datadist = "ddist")
m6 = lrm(arrest ~ age + prior, data = arr)
p = nomogram(m6, fun = function(x)1/(1+exp(-x)), fun.at = c(0.001, 0.01, 0.05, seq(0.1, 0.9, by = 0.1), 0.95, 0.99, 0.999), funlabel = "Probability of arrest")
plot(p, cex.axis = 0.6, lmgp = 0.2)
# 9: Bạn hãy ghi lại tất cả những hàm/lệnh trên trong RMarkdown và share
trên mạng rpubs.com/tài khoản của bạn.