Tai nạn phi thuyền Challenger

temp = c(66, 70, 69, 80, 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, NA, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0,
0, 0, 0, 0, 1, 0, 1)
dat = data.frame(temp, damage)
dat
##    temp damage
## 1    66      0
## 2    70      1
## 3    69      0
## 4    80     NA
## 5    68      0
## 6    67      0
## 7    72      0
## 8    73      0
## 9    70      0
## 10   57      1
## 11   63      1
## 12   70      1
## 13   78      0
## 14   67      0
## 15   53      1
## 16   67      0
## 17   75      0
## 18   70      0
## 19   81      0
## 20   76      0
## 21   79      0
## 22   75      1
## 23   76      0
## 24   58      1

hiển thị mối liên quan

library(ggplot2)
p = ggplot(data=dat, aes(x = temp, y = damage))
p + geom_point(alpha = 0.15) + geom_smooth(method = "glm",
method.args = list(family = "binomial"))
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Warning: Removed 1 rows containing missing values (geom_point).

# mô hình hồi qui logistic

m = glm(damage ~ temp, family=binomial, data=dat)
summary(m)
## 
## Call:
## glm(formula = damage ~ temp, family = binomial, data = dat)
## 
## 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
##   (1 observation deleted due to missingness)
## AIC: 24.315
## 
## Number of Fisher Scoring iterations: 5

Hiển thị kết quả mô hình hồi qui logistic

Có thể dùng package “GGally” để vẽ biểu đồ odds ratio cho yếu tố nguy cơ

library(GGally); library(ggplot2)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
ggcoef(m, exponentiate=T, exclude_intercept=T, vline_color = "red", errorbar_color = "blue", errorbar_height = 0.10)

ggcoef(m, exponentiate = TRUE, exclude_intercept=T, mapping = aes(x = estimate, y = term, size = p.value)) + scale_size_continuous(trans = "reverse")

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.