Bài thực hành 1
Nghiên cứu về smoking và ung thư phổi
Chúng ta cùng quay lại nghiên cứu giữa hút thuốc lá và ung thư phổi. Gọi P là xác suất bị ung thư phổi, smoke là biến số hút thuốc lá với giá trị 0 là không hút thuốc và giá trị 1 là hút thuốc.
cancer <- c(1,1,0,0)
smoke <- c(1,0,1,0)
n <- c(647,2,622, 27)
mydata <- data.frame(cancer, smoke, n)
mydata## cancer smoke n
## 1 1 1 647
## 2 1 0 2
## 3 0 1 622
## 4 0 0 27Bây giờ mô hình của chúng ta là: \[log\Big(\frac{P}{1-P}\Big)=\alpha + \beta *smoke\] Trong R, chúng ta sử dụng hàm glm() để xây dựng mô hình hồi quy Logistic
m <- glm(cancer ~ smoke, family = binomial, weights = n, data = mydata)
summary(m)##
## Call:
## glm(formula = cancer ~ smoke, family = binomial, data = mydata,
## weights = n)
##
## Deviance Residuals:
## 1 2 3 4
## 29.524 3.271 -29.783 -1.964
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.6027 0.7320 -3.556 0.000377 ***
## smoke 2.6421 0.7341 3.599 0.000319 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1799.4 on 3 degrees of freedom
## Residual deviance: 1773.3 on 2 degrees of freedom
## AIC: 1777.3
##
## Number of Fisher Scoring iterations: 4Ở nội dung trên, chúng ta sử dụng thêm đối số weight = n để glm biết rằng mỗi dòng trong dữ liệu không phải là một đối tượng mà là n đối tượng. Như vậy, phương trình hồi quy của chúng ta thu được: \[log\Big(\frac{P}{1-P}\Big)=-2.6 + 2.64 *smoke\] Hệ số hồi quy (log odds ration) là 2.64 với sai số chuẩn 0.73, trị số z = 3.6 cao hơn 2.0; và do đó có ý nghĩa thống kê (p=0.0003). Nói cách khác mối liên quan giữa ung thư phổi và hút thuốc lá có ý nghĩa thống kê. Chúng ta có: \[log(OR)=2.64\] \[OR=e^{2.64}=14.0\] Để tính khoảng tin cậy, chúng ta phải sử dụng sai số chuẩn 0.73. Phần dưới của khoảng tin cậy 95% là: \[e^{2.64-1.96*0.73}=3.35\] và phần trên: \[e^{2.64+1.96*0.73}=58.6\] Tuy nhiên, trong thực tế chúng ta có thể sử dụng package epiDisplay và hàm logistic.display để tính toán chính xác hơn và nhanh hơn:
library(epiDisplay)## Loading required package: foreign## Loading required package: survival##
## Attaching package: 'survival'## The following object is masked _by_ '.GlobalEnv':
##
## cancer## Loading required package: MASS## Loading required package: nnetlogistic.display(m)##
## Logistic regression predicting cancer
##
## OR(95%CI) P(Wald's test) P(LR-test)
## smoke: 1 vs 0 14.04 (3.33,59.2) < 0.001 < 0.001
##
## Log-likelihood = -886.6352
## No. of observations = 4
## AIC value = 1777.2704Kết quả trên cho thấy odds mắc ung thư phổi ở nhóm hút thuốc lá cao hơn nhóm không hút thuốc lá đến 14 lần, khoảng tin cậy 95% dao động từ 3.3 đến 59.2.
Bài thực hành số 2
Nghiên cứu về nồng độ chất béo trong nước xốt
Trong nghiên cứu này, biến dự đoán (nồng độ chất béo) là một biến liên tục và biến phụ thuộc là biến có dạng nhị phân (thích=1/không thích=0). Gọi biến phụ thuộc là “like” và biến dự đoán là “concen”, khởi tạo dữ liệu trong R như sau:
concen <- c(1.35, 1.6, 1.75, 1.85, 1.95, 2.05, 2.15, 2.25, 2.35)
like <- c(13,19, 67, 45, 71, 50, 35, 7, 1)
dislike <- c(0,0,2,5,8,20,31,49,12)
# Xac suat thich la
total <- like+dislike
prob <- round(like/total,2)
mydata <- data.frame(concen, like, dislike, prob)
head(mydata) ## concen like dislike prob
## 1 1.35 13 0 1.00
## 2 1.60 19 0 1.00
## 3 1.75 67 2 0.97
## 4 1.85 45 5 0.90
## 5 1.95 71 8 0.90
## 6 2.05 50 20 0.71# Xay dung bieu do tuong quan
plot(prob ~ concen, pch =16, xlab="Concentration", data = mydata, col ="red")
Biểu đồ tương quan này cho thấy, nồng độ chất béo càng tăng thì xác suất thích càng thấp. Tuy nhiên mối liên quan này không phải là đường thẳng tuyến tính mà là đường cong. Mô hình hồi quy có thể được miêu tả như sau: \[log\Big(\frac{P}{1-P}\Big)=\alpha + \beta *concen\] Ước tính tham số của mô hình này, sử dụng hàm glm()
m <- glm(prob ~ concen, family = binomial, weights = total, data = mydata)## Warning in eval(family$initialize): non-integer #successes in a binomial glm!summary(m)##
## Call:
## glm(formula = prob ~ concen, family = binomial, data = mydata,
## weights = total)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.85103 -0.73974 0.07971 0.36586 1.40775
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 22.733 2.266 10.032 <2e-16 ***
## concen -10.679 1.083 -9.862 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 199.9537 on 8 degrees of freedom
## Residual deviance: 9.0126 on 7 degrees of freedom
## AIC: 37.101
##
## Number of Fisher Scoring iterations: 5Kết quả cho thấy, mối liên quan trên có ý nghĩa thống kê. Phương trình thu được là: \[log\Big(\frac{P}{1-P}\Big)=22.733-10.679*concen\] Như vậy có thể hiểu là mỗi một đơn vị concen tăng 1 (tức 1 gram/L) thì OR thích nước xốt là exp(-10.679)=0.0000234, rất thấp đến mức không thực tế. Trong thực tế thì nồng độ chất béo tăng rất hạn chế trong khoảng 1.35 đến 2.35 gram/L. Do đó chúng ta cần tính OR trên 0.01 gram/L. Để làm được thì chúng ta sử dụng hàm I(concen/0.01) để cho biết hàm glm() sẽ có hệ số hồi quy được tính trên 0.01 gram/L
m1 <- glm(prob ~ I(concen/0.01), family = binomial, weights = total,
data = mydata)
summary(m1)##
## Call:
## glm(formula = prob ~ I(concen/0.01), family = binomial, data = mydata,
## weights = total)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.85103 -0.73974 0.07971 0.36586 1.40775
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 22.73345 2.26601 10.032 <2e-16 ***
## I(concen/0.01) -0.10679 0.01083 -9.862 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 199.9537 on 8 degrees of freedom
## Residual deviance: 9.0126 on 7 degrees of freedom
## AIC: 37.101
##
## Number of Fisher Scoring iterations: 5Phương trình hồi quy Logistic bây giờ là: \[log\Big(\frac{P}{1-P}\Big)=22.73345-0.10679 *concen\] Kết quả khi này cho thấy OR = exp(-0.10679) = 0.9, có nghĩa là mỗi 0.01 gram/L tăng nồng độ chất béo trong xốt, thì odds nước xốt giảm 10% (thích lúc này chỉ còn 90%).
plot(prob ~ concen, pch=16, data = mydata, col = "red", xlab = "Concentration")
lines(concen, m1$fitted, type = "l", col ="blue")
Biểu đồ cho thấy mô hình dự đoán khá tốt, khá khớp với bộ dữ liệu. Cùng tính phần dư và so sánh giá trị thực và giá trị dự đoán của mô hình
real_value <- mydata$prob
fitted_value <- m1$fitted.values
residual_value <- real_value - fitted_value
residual_value_percent = paste((residual_value/real_value)*100, "%")
reusult <- data.frame(real_value, fitted_value, residual_value, residual_value_percent)Bài thực hành số 3
Phi thuyền Challenger
Với nội dung cuối này, chúng ta sẽ sử dụng mô hình hồi quy Logistic để dự đoán về tai nạn của phi thuyền Challenger. Biến dự đoán của chúng ta ở đây là nhiệt độ, là một biến liên tục. Biến cần dự đoán, hay biến phụ thuộc là P- xác suất vòng O-Ring bị hỏng. Xây dựng mô hình hồi quy như sau: \[log\Big(\frac{P}{1-P}\Big)=\alpha+\beta*temp\]
temp <- c(66,70,69,68,67,72,73,70,57,63,70,78,67,53,67,75,70,81,76,79,75,76,58)
damage <- c(0,1,0,0,0,0,0,0,1,1,1,0,0,1,0,0,0,0,0,0,1,0,1)
data_challenger <- data.frame(temp, damage)
challenger_logistic_model <- glm(damage ~ temp, family = binomial, data = data_challenger)
summary(challenger_logistic_model)##
## Call:
## glm(formula = damage ~ temp, family = binomial, data = data_challenger)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0611 -0.7613 -0.3783 0.4524 2.2175
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 15.0429 7.3786 2.039 0.0415 *
## temp -0.2322 0.1082 -2.145 0.0320 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 28.267 on 22 degrees of freedom
## Residual deviance: 20.315 on 21 degrees of freedom
## AIC: 24.315
##
## Number of Fisher Scoring iterations: 5Mô hình thu được có phương trình như sau: \[log\Big(\frac{P}{1-P}\Big) = 15.0429 - 0.2322 *temp\] Như vậy mỗi một độ F tăng thì nguy cơ (odds) hư hỏng giảm 20%, vì có dấu trừ trong phương trình.
# Tính xác suất hư hỏng cho mỗi phi thuyền dựa vào mô hình này
data_challenger$predicted <- predict(challenger_logistic_model, type = "response") # Thêm tham số vào là response để kết quả đưa ra là xác suấtSử dụng thư viện trực quan ggplot2 để vẽ đường hồi quy Logistic
library(ggplot2)##
## Attaching package: 'ggplot2'## The following object is masked from 'package:epiDisplay':
##
## alphaggplot(data = data_challenger, aes(x = temp, y = damage)) +
geom_point() +
geom_line(aes(x=temp, y = predicted), col = "red", size = 0.4)