library(VIM)
## Loading required package: colorspace
## Loading required package: grid
## Loading required package: data.table
## VIM is ready to use. 
##  Since version 4.0.0 the GUI is in its own package VIMGUI.
## 
##           Please use the package to use the new (and old) GUI.
## Suggestions and bug-reports can be submitted at: https://github.com/alexkowa/VIM/issues
## 
## Attaching package: 'VIM'
## The following object is masked from 'package:datasets':
## 
##     sleep
library(mice)
## Loading required package: lattice
## 
## Attaching package: 'mice'
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
md.pattern(sleep)

##    BodyWgt BrainWgt Pred Exp Danger Sleep Span Gest Dream NonD   
## 42       1        1    1   1      1     1    1    1     1    1  0
## 9        1        1    1   1      1     1    1    1     0    0  2
## 3        1        1    1   1      1     1    1    0     1    1  1
## 2        1        1    1   1      1     1    0    1     1    1  1
## 1        1        1    1   1      1     1    0    1     0    0  3
## 1        1        1    1   1      1     1    0    0     1    1  2
## 2        1        1    1   1      1     0    1    1     1    0  2
## 2        1        1    1   1      1     0    1    1     0    0  3
##          0        0    0   0      0     4    4    4    12   14 38
aggr(sleep,prop=T, numbers=T)

matrixplot(sleep, interactive = F, sort= "Sleep")

# đoán giá trị trống, dùng hàm mice để đoán giá trị trống, seed: theo đúng thứ tụ ngẫu nhiên
msleep=mice(sleep, seed=1234, printFlag = F)
## Warning: Number of logged events: 5
# kiểm định giá trị imputed
msleep$imp$Dream
##      1   2   3   4   5
## 1  0.5 0.3 0.5 1.3 1.4
## 3  1.4 4.1 2.3 1.4 1.5
## 4  3.4 3.4 2.2 3.6 3.6
## 14 0.5 0.3 0.5 0.3 0.3
## 24 1.2 1.0 3.4 1.3 6.1
## 26 2.1 5.6 2.4 2.8 2.2
## 30 2.3 1.2 2.3 2.4 2.4
## 31 3.4 0.6 0.5 0.5 2.3
## 47 1.4 2.0 1.4 1.5 3.6
## 53 0.5 0.5 0.0 0.6 0.5
## 55 0.5 0.3 1.4 2.7 3.4
## 62 2.3 3.1 6.6 5.6 3.4
# tạo một data mới
isleep=complete(msleep, action=3)
head(isleep)
##    BodyWgt BrainWgt NonD Dream Sleep Span Gest Pred Exp Danger
## 1 6654.000   5712.0  3.3   0.5   3.3 38.6  645    3   5      3
## 2    1.000      6.6  6.3   2.0   8.3  4.5   42    3   1      3
## 3    3.385     44.5 10.4   2.3  12.5 14.0   60    1   1      1
## 4    0.920      5.7 14.3   2.2  16.5  2.0   25    5   2      3
## 5 2547.000   4603.0  2.1   1.8   3.9 69.0  624    3   5      4
## 6   10.550    179.5  9.1   0.7   9.8 27.0  180    4   4      4
# bước 5 phân tích bằng mô hình

dư liệu hipfacture

t="C:\\Users\\pc\\Downloads\\CAC NGHIEN CUU\\BAI GIANG CAC MON\\GS TUAN\\BG 12.6.19\\Hip fracture data.csv"
dat=read.csv(t,na.strings = "")
#chú thích: na.srings là báo cho mình có dữ liệu NA
dim(dat)
## [1] 2788   20
#chú thích: có 2788 dữ liệu và 20 dữ liệu NA. 
#chọn những đối tượng mà v2<2 đưa vào dataframe"hh" hay tạo ra bộ dữ liệu nhỏ hơn là dat1 là subset của dat. 
dat1=subset(dat, v2<2.0)
# tạo tiếp một subset hip chỉ gồm một số biến từ hipfx đến bmi
hip=dat1[, c("hipfx", "gender", "age", "v1", "v2", "v3", "v4", "v5", "v6", "v7", "v8", "v9", "bmi")]

imputation

library(mice) 
library(VIM)
#xem phần có missing data. 
md.pattern(hip)

##      hipfx gender age v2 v3 v7 v8 v9 v1 bmi v4  v6  v5    
## 2315     1      1   1  1  1  1  1  1  1   1  1   1   1   0
## 229      1      1   1  1  1  1  1  1  1   1  1   1   0   1
## 12       1      1   1  1  1  1  1  1  1   1  1   0   1   1
## 125      1      1   1  1  1  1  1  1  1   1  1   0   0   2
## 29       1      1   1  1  1  1  1  1  1   1  0   1   1   1
## 1        1      1   1  1  1  1  1  1  1   1  0   1   0   2
## 2        1      1   1  1  1  1  1  1  1   1  0   0   0   3
## 3        1      1   1  1  1  1  1  1  1   0  1   1   1   1
## 4        1      1   1  1  1  1  1  1  1   0  1   0   0   3
## 1        1      1   1  1  1  1  1  1  0   1  1   1   1   1
## 1        1      1   1  1  1  1  1  1  0   1  1   0   0   3
##          0      0   0  0  0  0  0  0  2   7 32 144 362 547
ms.hip=mice(hip, seed=1234, printFlag = F)
ihip=complete(ms.hip)
#ihip: là có tất cả NA

#chia data thành 2 nhóm training và testing dùng package caret, vơi chỉ biến hip là có NA nên training và testing trên biến hip.

library(caret)
## Warning: package 'caret' was built under R version 3.5.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.5.3
## 
## Attaching package: 'ggplot2'
## The following object is masked _by_ '.GlobalEnv':
## 
##     msleep
index=createDataPartition(ihip$hipfx, p=0.7, list=F)
training=ihip[index, ]
testing=ihip[-index,]
dim(ihip)
## [1] 2722   13
dim(testing)
## [1] 816  13

xây dựng mô hình trên training data

# chỉ đưa hai biến age và v3 là hai biến có liên quan đến hipfx, fit là dữ liệu training
fit=train(hipfx~age+v3, data=training, method="glm", family=binomial)
## Warning in train.default(x, y, weights = w, ...): You are trying to do
## regression and your outcome only has two possible values Are you trying to
## do classification? If so, use a 2 level factor as your outcome column.
summary(fit)
## 
## Call:
## NULL
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6467  -0.3399  -0.2302  -0.1410   3.4351  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.53502    1.38010  -0.388    0.698    
## age          0.05757    0.01410   4.082 4.46e-05 ***
## v3          -7.96972    0.85465  -9.325  < 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: 917.41  on 1905  degrees of freedom
## Residual deviance: 725.57  on 1903  degrees of freedom
## AIC: 731.57
## 
## Number of Fisher Scoring iterations: 6

Test model trên testing data

#Tìm giá trị tiên lượng
testing$pred=predict(fit, newdata = testing, type="raw")
head(testing)
##    hipfx gender age   v1   v2    v3    v4 v5   v6 v7 v8 v9 bmi        pred
## 2      0 Female  67 0.85 0.85 0.966 1.325 18 31.0  0  0  0  26 0.012409233
## 5      0   Male  61 0.87 0.60 0.811 1.144 44 28.9  1  0  0  24 0.029686736
## 6      0 Female  76 0.76 0.58 0.743 0.980 15 33.3  0  0  0  28 0.110909559
## 7      0   Male  64 0.94 0.79 1.006 1.376 26 33.4  1  0  0  32 0.007627676
## 9      0 Female  76 0.50 0.46 0.654 0.874 20 29.7  0  0  0  21 0.202268050
## 13     0   Male  62 1.16 0.94 1.155 1.441 31 30.5  0  0  0  33 0.002084924
# so sánh giá trị thực tế alf hipfx và giá trị dự đoán
library(pROC)
## Warning: package 'pROC' was built under R version 3.5.3
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following object is masked from 'package:colorspace':
## 
##     coords
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
mm=roc(testing$hipfx, testing$pred)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
auc(mm)
## Area under the curve: 0.8638
ci(mm)
## 95% CI: 0.8168-0.9109 (DeLong)
plot(mm)