Thuc hanh ngay 4: Mô hình hồi qui Logistic
Nguyen Ngoc Tuyen
8/7/2019
Task 1: Phân tích mô tả
Nhập dữ liệu
setwd("C:/Users\\HP\\Documents\\R\\Thuc hanh\\Training CR")
t=file.choose()
os=read.csv(t,header=T)
head(os)## id dov gender age dob visit v1 v2 v3 v4 wt bmi ht v5
## 1 3 15/6/89 Male 73 8/6/16 1 0.98 0.88 1.079 1.458 98 32 175 NA
## 2 8 17/4/89 Female 67 11/12/21 1 0.85 0.85 0.966 1.325 72 26 166 18
## 3 9 12/6/90 Male 68 8/1/22 1 0.87 0.84 1.013 1.494 87 26 184 36
## 4 10 4/6/90 Female 62 15/5/28 1 0.62 0.71 0.839 1.214 72 24 173 NA
## 5 23 8/8/89 Male 61 22/9/28 1 0.87 0.60 0.811 1.144 72 24 173 44
## 6 24 3/5/89 Female 76 1/8/13 1 0.76 0.58 0.743 0.980 67 28 156 15
## v6 v7 v8 v9 hipfx timehip
## 1 39.9 1 0 0 0 0.55
## 2 31.0 0 0 0 0 19.68
## 3 28.6 0 0 0 0 5.05
## 4 28.2 1 0 0 0 18.55
## 5 28.9 1 0 0 0 19.37
## 6 33.3 0 0 0 0 12.30attach(os)Phân tích mô tả các biến được chọn theo biến hipfx
require(table1)## Loading required package: table1## Warning: package 'table1' was built under R version 3.5.3##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-table1(~gender+age+v1+v2+v3+v4+wt+ht+bmi|hipfx,data=os)## Warning in table1.formula(~gender + age + v1 + v2 + v3 + v4 + wt + ht
## + : Terms to the right of '|' in formula 'x' define table columns and are
## expected to be factors with meaningful labels.| 0 (n=2599) |
1 (n=189) |
Overall (n=2788) |
|
|---|---|---|---|
| gender | |||
| Female | 1512 (58.2%) | 142 (75.1%) | 1654 (59.3%) |
| Male | 1087 (41.8%) | 47 (24.9%) | 1134 (40.7%) |
| age | |||
| Mean (SD) | 68.9 (6.35) | 75.3 (7.74) | 69.3 (6.65) |
| Median [Min, Max] | 67.0 [60.0, 94.0] | 75.0 [60.0, 96.0] | 68.0 [60.0, 96.0] |
| v1 | |||
| Mean (SD) | 0.820 (0.310) | 0.637 (0.125) | 0.809 (0.305) |
| Median [Min, Max] | 0.810 [0.0900, 12.2] | 0.620 [0.430, 0.970] | 0.800 [0.0900, 12.2] |
| Missing | 51 (2.0%) | 16 (8.5%) | 67 (2.4%) |
| v2 | |||
| Mean (SD) | 0.744 (1.89) | 0.520 (0.129) | 0.730 (1.83) |
| Median [Min, Max] | 0.700 [0.220, 96.0] | 0.500 [0.170, 1.01] | 0.690 [0.170, 96.0] |
| Missing | 50 (1.9%) | 15 (7.9%) | 65 (2.3%) |
| v3 | |||
| Mean (SD) | 0.874 (0.150) | 0.692 (0.123) | 0.862 (0.155) |
| Median [Min, Max] | 0.866 [0.378, 1.51] | 0.683 [0.279, 1.16] | 0.855 [0.279, 1.51] |
| Missing | 50 (1.9%) | 15 (7.9%) | 65 (2.3%) |
| v4 | |||
| Mean (SD) | 1.17 (0.226) | 0.996 (0.218) | 1.16 (0.230) |
| Median [Min, Max] | 1.15 [0.389, 2.38] | 0.976 [0.358, 1.83] | 1.14 [0.358, 2.38] |
| Missing | 56 (2.2%) | 11 (5.8%) | 67 (2.4%) |
| wt | |||
| Mean (SD) | 74.1 (15.3) | 62.4 (13.3) | 73.3 (15.4) |
| Median [Min, Max] | 73.0 [36.0, 149] | 60.5 [35.0, 101] | 72.0 [35.0, 149] |
| Missing | 19 (0.7%) | 13 (6.9%) | 32 (1.1%) |
| ht | |||
| Mean (SD) | 166 (9.23) | 161 (9.72) | 165 (9.34) |
| Median [Min, Max] | 165 [143, 192] | 160 [136, 190] | 165 [136, 192] |
| Missing | 21 (0.8%) | 13 (6.9%) | 34 (1.2%) |
| bmi | |||
| Mean (SD) | 26.9 (4.73) | 23.9 (3.86) | 26.7 (4.74) |
| Median [Min, Max] | 26.0 [15.0, 57.0] | 24.0 [15.0, 36.0] | 26.0 [15.0, 57.0] |
| Missing | 21 (0.8%) | 13 (6.9%) | 34 (1.2%) |
Task 2: Mô hình hồi quy Logistic đơn giản
m=glm(hipfx~gender,family = binomial,data=os)
summary(m)##
## Call:
## glm(formula = hipfx ~ gender, family = binomial, data = os)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4237 -0.4237 -0.4237 -0.2910 2.5232
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.36536 0.08777 -26.949 < 2e-16 ***
## genderMale -0.77567 0.17289 -4.486 7.24e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1382.2 on 2787 degrees of freedom
## Residual deviance: 1360.0 on 2786 degrees of freedom
## AIC: 1364
##
## Number of Fisher Scoring iterations: 5library(epiDisplay)## Warning: package 'epiDisplay' was built under R version 3.5.3## Loading required package: foreign## Loading required package: survival## Loading required package: MASS## Warning: package 'MASS' was built under R version 3.5.3## Loading required package: nnetlogistic.display(m)##
## Logistic regression predicting hipfx
##
## OR(95%CI) P(Wald's test) P(LR-test)
## gender: Male vs Female 0.46 (0.33,0.65) < 0.001 < 0.001
##
## Log-likelihood = -679.98
## No. of observations = 2788
## AIC value = 1363.96Task 3: Chọn biến bằng BMA
Định nghĩa biến x và y
xvars=os[,c("gender", "age", "wt", "ht", "bmi", "v1", "v2", "v3", "v4")]
yvar = os$hipfxTìm biến liên quan
t=file.choose()
os=read.csv(t,header=T)
require(BMA)## Loading required package: BMA## Warning: package 'BMA' was built under R version 3.5.3## Loading required package: leaps## Warning: package 'leaps' was built under R version 3.5.3## Loading required package: robustbase## Warning: package 'robustbase' was built under R version 3.5.3##
## Attaching package: 'robustbase'## The following object is masked from 'package:survival':
##
## heart## Loading required package: inline## Warning: package 'inline' was built under R version 3.5.3## Loading required package: rrcov## Warning: package 'rrcov' was built under R version 3.5.3## Scalable Robust Estimators with High Breakdown Point (version 1.4-7)xvars=os[,c("gender", "age", "wt", "ht", "bmi", "v1", "v2", "v3", "v4")]
yvar = os$hipfx
mh1=bic.glm(xvars,yvar,strict = F,OR=20,glm.family ="binomial")## Warning in bic.glm.data.frame(xvars, yvar, strict = F, OR = 20, glm.family
## = "binomial"): There were 106 records deleted due to NA's## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredsummary(mh1)##
## Call:
## bic.glm.data.frame(x = xvars, y = yvar, glm.family = "binomial", strict = F, OR = 20)
##
##
## 5 models were selected
## Best 5 models (cumulative posterior probability = 1 ):
##
## p!=0 EV SD model 1 model 2
## Intercept 100 -0.2193903 1.369781 -1.060e-01 -2.015e-01
## gender.x 0.0
## .Male 0.0000000 0.000000 . .
## age.x 100.0 0.0555837 0.011979 5.535e-02 5.718e-02
## wt.x 0.0 0.0000000 0.000000 . .
## ht.x 4.8 0.0008384 0.004426 . .
## bmi.x 6.4 -0.0025378 0.011323 . .
## v1.x 11.1 -0.2275109 0.728763 . -2.059e+00
## v2.x 6.5 -0.1607453 0.711598 . .
## v3.x 100.0 -7.9911353 1.178497 -8.357e+00 -6.505e+00
## v4.x 0.0 0.0000000 0.000000 . .
##
## nVar 2 3
## BIC -2.015e+04 -2.015e+04
## post prob 0.712 0.111
## model 3 model 4 model 5
## Intercept -2.962e-01 4.742e-01 -2.750e+00
## gender.x
## .Male . . .
## age.x 5.399e-02 5.519e-02 5.805e-02
## wt.x . . .
## ht.x . . 1.737e-02
## bmi.x . -3.964e-02 .
## v1.x . . .
## v2.x -2.482e+00 . .
## v3.x -6.050e+00 -7.794e+00 -8.860e+00
## v4.x . . .
##
## nVar 3 3 3
## BIC -2.015e+04 -2.015e+04 -2.015e+04
## post prob 0.065 0.064 0.048library(BMA)
imageplot.bma(mh1)
Đánh giá tầm quan trọng
mh2=glm(hipfx~age+v3,family = binomial,data=os)
summary(mh2)##
## Call:
## glm(formula = hipfx ~ age + v3, family = binomial, data = os)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6967 -0.3371 -0.2208 -0.1348 3.4804
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.07342 1.16953 -0.063 0.95
## age 0.05435 0.01183 4.594 4.35e-06 ***
## v3 -8.31357 0.73794 -11.266 < 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: 1293.8 on 2722 degrees of freedom
## Residual deviance: 1016.0 on 2720 degrees of freedom
## (65 observations deleted due to missingness)
## AIC: 1022
##
## Number of Fisher Scoring iterations: 6require(caret)## Loading required package: caret## Warning: package 'caret' was built under R version 3.5.3## Loading required package: lattice##
## Attaching package: 'lattice'## The following object is masked from 'package:epiDisplay':
##
## dotplot## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 3.5.3##
## Attaching package: 'ggplot2'## The following object is masked from 'package:epiDisplay':
##
## alpha##
## Attaching package: 'caret'## The following object is masked from 'package:survival':
##
## clustervarImp(mh2,scale=F)## Overall
## age 4.593955
## v3 11.265944