library(foreign)
## Warning: package 'foreign' was built under R version 3.2.5
tr=read.spss("C:/Users/BINH THANG/Dropbox/R - Learning/trung/SO-LIEU-NGHIEN-CUU-XU-LY-01.6.sav", to.data.frame=T)
## Warning in read.spss("C:/Users/BINH THANG/Dropbox/R - Learning/trung/
## SO-LIEU-NGHIEN-CUU-XU-LY-01.6.sav", : C:/Users/BINH THANG/Dropbox/R -
## Learning/trung/SO-LIEU-NGHIEN-CUU-XU-LY-01.6.sav: Unrecognized record type
## 7, subtype 18 encountered in system file
a=subset(tr, select=c("HOIQUY", "C2", "C4","C7","C8","C10", "C11" , "C12", "C13" , "C14" , "C15", "C17", "C18", "C19X1", "C19X2", "C19X3", "C19X4","C19X5" , "C19X6" , "B1" , "E1"))
dinh nghiax bien so
library(plyr)
## Warning: package 'plyr' was built under R version 3.2.5
a$HOIQUY <- mapvalues(a$HOIQUY, from = c("DUONG TINH", "AM TINH"), to = c("1", "0"))
a$HOIQUY=as.factor(a$HOIQUY)
a$C2=as.factor(a$C2)
a$C4=as.factor(a$C4)
a$C7=as.factor(a$C7)
a$C8=as.factor(a$C8)
a$C10=as.factor(a$C10)
a$C11=as.factor(a$C11)
a$C12=as.factor(a$C12)
a$C13=as.factor(a$C13)
a$C14=as.factor(a$C14)
a$C15=as.factor(a$C15)
a$C17=as.factor(a$C17)
a$C18=as.factor(a$C18)
a$C19X1=as.factor(a$C19X1)
a$C19X2=as.factor(a$C19X2)
a$C19X3=as.factor(a$C19X3)
a$C19X4=as.factor(a$C19X4)
a$C19X5=as.factor(a$C19X5)
a$C19X6=as.factor(a$C19X6)
a$B1=as.factor(a$B1)
a$E1=as.factor(a$E1)
a$C2=relevel(a$C2, ref='2')
a$C4=relevel(a$C4, ref='2')
a$C7=relevel(a$C7, ref='2')
a$C8=relevel(a$C8, ref='2')
a$C10=relevel(a$C10, ref='2')
a$C11=relevel(a$C11, ref='2')
a$C12=relevel(a$C12, ref='2')
a$C13=relevel(a$C13, ref='2')
a$C14=relevel(a$C14, ref='2')
a$C15=relevel(a$C15, ref='2')
a$C17=relevel(a$C17, ref='2')
a$C18=relevel(a$C18, ref='2')
a$C19X1=relevel(a$C19X1, ref='2')
a$C19X2=relevel(a$C19X2, ref='2')
a$C19X3=relevel(a$C19X3, ref='2')
a$C19X4=relevel(a$C19X4, ref='2')
a$C19X5=relevel(a$C19X5, ref='2')
a$C19X6=relevel(a$C19X6, ref='2')
a$B1=relevel(a$B1, ref='KHONG NHON')
a$E1=as.factor(a$E1)
phuong phap 1: su dung mo hinh BMA dua vao BIC
yvar = a[,1]
xvars = a[,-1]
#install.packages("BMA")
library(BMA)
## Warning: package 'BMA' was built under R version 3.2.5
## Loading required package: survival
## Warning: package 'survival' was built under R version 3.2.5
## Loading required package: leaps
## Warning: package 'leaps' was built under R version 3.2.5
## Loading required package: robustbase
## Warning: package 'robustbase' was built under R version 3.2.5
##
## 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.2.5
## Loading required package: rrcov
## Warning: package 'rrcov' was built under R version 3.2.5
## Scalable Robust Estimators with High Breakdown Point (version 1.4-3)
bma.search=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 2 records deleted due to NA's
summary(bma.search)
##
## Call:
## bic.glm.data.frame(x = xvars, y = yvar, glm.family = "binomial", strict = F, OR = 20)
##
##
## 21 models were selected
## Best 5 models (cumulative posterior probability = 0.6601 ):
##
## p!=0 EV SD model 1 model 2 model 3
## Intercept 100 -1.447556 0.21882 -1.3573 -1.3768 -1.6493
## C2.x 4.2
## .1 0.010487 0.06264 . . .
## C4.x 7.4
## .1 0.024620 0.10146 . . .
## C7.x 13.6
## .1 0.053149 0.15300 . . .
## C8.x 0.0
## .1 0.000000 0.00000 . . .
## C10.x 0.0
## .1 0.000000 0.00000 . . .
## C11.x 0.0
## .1 0.000000 0.00000 . . .
## C12.x 100.0
## .1 0.647483 0.18035 0.6654 0.6673 0.6509
## C13.x 5.1
## .1 0.011971 0.06541 . . .
## C14.x 0.0
## .1 0.000000 0.00000 . . .
## C15.x 3.9
## .1 0.009282 0.05850 . . .
## C17.x 0.0
## .1 0.000000 0.00000 . . .
## C18.x 0.0
## .1 0.000000 0.00000 . . .
## C19X1.x 98.6
## .1 -0.764524 0.24405 -0.7553 -0.7980 -0.7464
## C19X2.x 0.0
## .1 0.000000 0.00000 . . .
## C19X3.x 43.0
## .1 0.713542 0.94536 . 1.6435 .
## C19X4.x 1.3
## .1 0.003295 0.04555 . . .
## C19X5.x 1.6
## .1 -0.010786 0.10962 . . .
## C19X6.x 2.6
## .1 -0.345474 100.95638 . . .
## B1.x 100.0
## .NHON 0.855789 0.17614 0.8640 0.8600 0.8640
## E1.x 19.7
## .1 0.080864 0.18498 . . 0.4185
##
## nVar 3 4 4
## BIC -3405.3397 -3404.9470 -3403.6235
## post prob 0.232 0.191 0.098
## model 4 model 5
## Intercept -1.4055 -1.6575
## C2.x
## .1 . .
## C4.x
## .1 . .
## C7.x
## .1 0.4039 .
## C8.x
## .1 . .
## C10.x
## .1 . .
## C11.x
## .1 . .
## C12.x
## .1 0.6015 0.6522
## C13.x
## .1 . .
## C14.x
## .1 . .
## C15.x
## .1 . .
## C17.x
## .1 . .
## C18.x
## .1 . .
## C19X1.x
## .1 -0.7899 -0.7890
## C19X2.x
## .1 . .
## C19X3.x
## .1 . 1.6245
## C19X4.x
## .1 . .
## C19X5.x
## .1 . .
## C19X6.x
## .1 . .
## B1.x
## .NHON 0.8401 0.8587
## E1.x
## .1 . 0.4055
##
## nVar 4 5
## BIC -3402.9569 -3402.9112
## post prob 0.070 0.069
#logistic model 1
log <- glm(HOIQUY ~ C12+ C19X1+ C19X3+ B1+ E1 ,family=binomial(logit),data=a)
library("epiDisplay")
## Warning: package 'epiDisplay' was built under R version 3.2.5
## Loading required package: MASS
## Loading required package: nnet
logistic.display(log)
##
## Logistic regression predicting HOIQUY : 1 vs 0
##
## crude OR(95%CI) adj. OR(95%CI) P(Wald's test)
## C12: 1 vs 2 1.68 (1.2,2.34) 1.89 (1.34,2.69) < 0.001
##
## C19X1: 1 vs 2 0.56 (0.37,0.86) 0.47 (0.3,0.73) < 0.001
##
## C19X3: 1 vs 2 4.71 (1.21,18.4) 5.02 (1.24,20.41) 0.024
##
## B1: NHON vs KHONG NHON 2.17 (1.56,3.02) 2.37 (1.68,3.35) < 0.001
##
## E1: 1 vs 0 1.53 (1.06,2.21) 1.47 (1.01,2.16) 0.045
##
## P(LR-test)
## C12: 1 vs 2 < 0.001
##
## C19X1: 1 vs 2 < 0.001
##
## C19X3: 1 vs 2 0.017
##
## B1: NHON vs KHONG NHON < 0.001
##
## E1: 1 vs 0 0.042
##
## Log-likelihood = -389.8063
## No. of observations = 653
## AIC value = 791.6127
require(moonBook)
## Loading required package: moonBook
## Warning: package 'moonBook' was built under R version 3.2.5
library(survival)
ORplot(log,type=0,show.CI=TRUE,xlab="Odd ratio, 95% Confidence Interval",main="Odds Ratio")
phuong phap 2: dua vao bien so quan trong trong mo hinh
#mo hinh toan bo bien so
log <- glm(HOIQUY ~ .,family=binomial(logit),data=a)
#danh gia bien nao quan trong
library(caret)
## Warning: package 'caret' was built under R version 3.2.5
## Loading required package: lattice
##
## Attaching package: 'lattice'
## The following object is masked from 'package:moonBook':
##
## densityplot
## The following object is masked from 'package:epiDisplay':
##
## dotplot
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.2.5
##
## Attaching package: 'ggplot2'
## The following object is masked from 'package:epiDisplay':
##
## alpha
##
## Attaching package: 'caret'
## The following object is masked from 'package:survival':
##
## cluster
varImp(log, scale = FALSE)
## Overall
## C21 0.75272180
## C41 1.06637385
## C71 1.94170650
## C81 1.01826345
## C101 1.00155071
## C111 0.55463437
## C121 2.73462723
## C131 1.72317650
## C141 0.89078212
## C151 1.23966684
## C171 0.34198109
## C181 0.68381609
## C19X11 2.96911709
## C19X21 0.20271750
## C19X31 2.35528642
## C19X41 0.06406907
## C19X51 0.50852179
## C19X61 0.02140906
## B1NHON 4.33824727
## E11 1.77677737
#chon lai cac bien cho mo hinh cuoi cung
log1 <- glm(HOIQUY ~ C7+ C12+ C13+ C15+ C19X1+ C19X3+ B1+ E1 ,family=binomial(logit),data=a)
library("epiDisplay")
logistic.display(log1)
##
## Logistic regression predicting HOIQUY : 1 vs 0
##
## crude OR(95%CI) adj. OR(95%CI) P(Wald's test)
## C7: 1 vs 2 1.67 (1.16,2.41) 1.43 (0.97,2.13) 0.073
##
## C12: 1 vs 2 1.68 (1.2,2.34) 1.65 (1.14,2.39) 0.008
##
## C13: 1 vs 2 1.15 (0.83,1.59) 1.31 (0.93,1.86) 0.124
##
## C15: 1 vs 2 1.58 (1.14,2.21) 1.26 (0.88,1.81) 0.211
##
## C19X1: 1 vs 2 0.56 (0.37,0.86) 0.45 (0.29,0.71) < 0.001
##
## C19X3: 1 vs 2 4.71 (1.21,18.4) 4.32 (1.04,18.04) 0.045
##
## B1: NHON vs KHONG NHON 2.17 (1.56,3.02) 2.38 (1.67,3.37) < 0.001
##
## E1: 1 vs 0 1.53 (1.06,2.21) 1.45 (0.99,2.14) 0.057
##
## P(LR-test)
## C7: 1 vs 2 0.074
##
## C12: 1 vs 2 0.007
##
## C13: 1 vs 2 0.123
##
## C15: 1 vs 2 0.211
##
## C19X1: 1 vs 2 < 0.001
##
## C19X3: 1 vs 2 0.034
##
## B1: NHON vs KHONG NHON < 0.001
##
## E1: 1 vs 0 0.055
##
## Log-likelihood = -386.3921
## No. of observations = 653
## AIC value = 790.7842
ORplot(log1,type=0,show.CI=TRUE,xlab="Odd ratio, 95% Confidence Interval",main="Odds Ratio")
