library(outbreaks)
## Warning: package 'outbreaks' was built under R version 3.4.4
data("fluH7N9_china_2013")

class(fluH7N9_china_2013)
## [1] "data.frame"
View(fluH7N9_china_2013)
str(fluH7N9_china_2013)
## 'data.frame':    136 obs. of  8 variables:
##  $ case_id                : Factor w/ 136 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ date_of_onset          : Date, format: "2013-02-19" "2013-02-27" ...
##  $ date_of_hospitalisation: Date, format: NA "2013-03-03" ...
##  $ date_of_outcome        : Date, format: "2013-03-04" "2013-03-10" ...
##  $ outcome                : Factor w/ 2 levels "Death","Recover": 1 1 1 NA 2 1 1 1 NA 1 ...
##  $ gender                 : Factor w/ 2 levels "f","m": 2 2 1 1 1 1 2 2 2 2 ...
##  $ age                    : Factor w/ 61 levels "?","15","2","21",..: 58 7 11 18 20 9 54 14 39 20 ...
##  $ province               : Factor w/ 13 levels "Anhui","Beijing",..: 11 11 1 8 8 8 8 13 13 11 ...
head(fluH7N9_china_2013)
##   case_id date_of_onset date_of_hospitalisation date_of_outcome outcome
## 1       1    2013-02-19                    <NA>      2013-03-04   Death
## 2       2    2013-02-27              2013-03-03      2013-03-10   Death
## 3       3    2013-03-09              2013-03-19      2013-04-09   Death
## 4       4    2013-03-19              2013-03-27            <NA>    <NA>
## 5       5    2013-03-19              2013-03-30      2013-05-15 Recover
## 6       6    2013-03-21              2013-03-28      2013-04-26   Death
##   gender age province
## 1      m  87 Shanghai
## 2      m  27 Shanghai
## 3      f  35    Anhui
## 4      f  45  Jiangsu
## 5      f  48  Jiangsu
## 6      f  32  Jiangsu
fluH7N9_china_2013$age[fluH7N9_china_2013$age == '?'] <- NA
fluH7N9_china_2013$age <- as.numeric(fluH7N9_china_2013$age)

fluH7N9_china_2013$age
##   [1] 58  7 11 18 20  9 54 14 39 20 36 24 39 15 34 51 46 38 31 27 39 56  5
##  [24] 36 35 49 23 51 48 53 43 46 37 46 54 40  8 28 38 46 26 25 57 42 28 49
##  [47] 44 37 14 10 37 36 35 47 51 45 26 50 22  6 33 40 33 28  4 44 28 29 35
##  [70] 30 44 19 41 NA NA 27  3 59 13 46 57 16 14  6 52 26 41 15 26 17 20 38
##  [93] 28 11 13 17 47 48 40 30 51 53 40 26  9 12 61 55 35 25 28 41 33 22 14
## [116] 37 48 21 12 33 36 14 26 52  8 52 15 30 41 41 60 51 32  2 34 23
fluH7N9_china_2013$case_id <- paste('case',fluH7N9_china_2013$case_id, sep = '_')


library(tidyr)
## Warning: package 'tidyr' was built under R version 3.4.4
fluH7N9_china_2013_gather  <- fluH7N9_china_2013 %>%
    gather(Group, Date, date_of_onset:date_of_outcome)

fluH7N9_china_2013[fluH7N9_china_2013$case_id == 'case_1',]
##   case_id date_of_onset date_of_hospitalisation date_of_outcome outcome
## 1  case_1    2013-02-19                    <NA>      2013-03-04   Death
##   gender age province
## 1      m  58 Shanghai
fluH7N9_china_2013_gather[fluH7N9_china_2013$case_id == 'case_1',]
##     case_id outcome gender age province                   Group       Date
## 1    case_1   Death      m  58 Shanghai           date_of_onset 2013-02-19
## 137  case_1   Death      m  58 Shanghai date_of_hospitalisation       <NA>
## 273  case_1   Death      m  58 Shanghai         date_of_outcome 2013-03-04
nrow(fluH7N9_china_2013)
## [1] 136
nrow(fluH7N9_china_2013_gather)
## [1] 408
fluH7N9_china_2013_gather$Group <- factor(fluH7N9_china_2013_gather$Group, levels =
c("date_of_onset", "date_of_hospitalisation", "date_of_outcome"))

str(fluH7N9_china_2013_gather)
## 'data.frame':    408 obs. of  7 variables:
##  $ case_id : chr  "case_1" "case_2" "case_3" "case_4" ...
##  $ outcome : Factor w/ 2 levels "Death","Recover": 1 1 1 NA 2 1 1 1 NA 1 ...
##  $ gender  : Factor w/ 2 levels "f","m": 2 2 1 1 1 1 2 2 2 2 ...
##  $ age     : num  58 7 11 18 20 9 54 14 39 20 ...
##  $ province: Factor w/ 13 levels "Anhui","Beijing",..: 11 11 1 8 8 8 8 13 13 11 ...
##  $ Group   : Factor w/ 3 levels "date_of_onset",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Date    : Date, format: "2013-02-19" "2013-02-27" ...
library(plyr)
## Warning: package 'plyr' was built under R version 3.4.4
fluH7N9_china_2013_gather$Group <-mapvalues(fluH7N9_china_2013_gather$Group, from =
c("date_of_onset", "date_of_hospitalisation", "date_of_outcome"),
to = c("Date of onset", "Date of hospitalisation", "Date of outcome"))


levels(fluH7N9_china_2013_gather$province)
##  [1] "Anhui"     "Beijing"   "Fujian"    "Guangdong" "Hebei"    
##  [6] "Henan"     "Hunan"     "Jiangsu"   "Jiangxi"   "Shandong" 
## [11] "Shanghai"  "Taiwan"    "Zhejiang"
fluH7N9_china_2013_gather$province <- mapvalues(fluH7N9_china_2013_gather$province,from =
c("Anhui", "Beijing", "Fujian", "Guangdong", "Hebei", "Henan", "Hunan", "Jiangxi", "Shandong", "Taiwan"),
to = rep("Other", 10))

levels(fluH7N9_china_2013_gather$province)
## [1] "Other"    "Jiangsu"  "Shanghai" "Zhejiang"
levels(fluH7N9_china_2013_gather$gender) <- c(levels(fluH7N9_china_2013_gather$gender),
"unknown")

fluH7N9_china_2013_gather$gender[is.na(fluH7N9_china_2013_gather$gender)] <- "unknown"

fluH7N9_china_2013_gather$province <- factor(fluH7N9_china_2013_gather$province, levels =
c("Jiangsu", "Shanghai", "Zhejiang", "Other"))

str(fluH7N9_china_2013_gather)
## 'data.frame':    408 obs. of  7 variables:
##  $ case_id : chr  "case_1" "case_2" "case_3" "case_4" ...
##  $ outcome : Factor w/ 2 levels "Death","Recover": 1 1 1 NA 2 1 1 1 NA 1 ...
##  $ gender  : Factor w/ 3 levels "f","m","unknown": 2 2 1 1 1 1 2 2 2 2 ...
##  $ age     : num  58 7 11 18 20 9 54 14 39 20 ...
##  $ province: Factor w/ 4 levels "Jiangsu","Shanghai",..: 2 2 4 1 1 1 1 3 3 2 ...
##  $ Group   : Factor w/ 3 levels "Date of onset",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Date    : Date, format: "2013-02-19" "2013-02-27" ...

資料視覺化

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.4.4
ggplot(data = fluH7N9_china_2013_gather, aes(x = Date, y = age, fill = outcome)) +
stat_density2d(aes(alpha = ..level..), geom = "polygon") +
geom_jitter(aes(color = outcome, shape = gender), size = 1.5) +
geom_rug(aes(color = outcome)) +
labs(
fill = "Outcome",
color = "Outcome",
alpha = "Level",
shape = "Gender",
x = "Date in 2013",
y = "Age",
title = "2013 Influenza A H7N9 cases in China"
) 
## Warning: Removed 149 rows containing non-finite values (stat_density2d).
## Warning: Removed 149 rows containing missing values (geom_point).

ggplot(data = fluH7N9_china_2013_gather, aes(x = Date, y = age, fill = outcome)) +
stat_density2d(aes(alpha = ..level..), geom = "polygon") +
geom_jitter(aes(color = outcome, shape = gender), size = 1.5) +
geom_rug(aes(color = outcome)) +
labs(
fill = "Outcome",
color = "Outcome",
alpha = "Level",
shape = "Gender",
x = "Date in 2013",
y = "Age",
title = "2013 Influenza A H7N9 cases in China"
) +
facet_grid(Group ~ province) +
scale_shape_manual(values = c(15, 16, 17)) +
scale_color_brewer(palette="Set1", na.value = "grey50") +
scale_fill_brewer(palette="Set1")
## Warning: Removed 149 rows containing non-finite values (stat_density2d).

## Warning: Removed 149 rows containing missing values (geom_point).

table(fluH7N9_china_2013$outcome, fluH7N9_china_2013$gender)
##          
##            f  m
##   Death    9 22
##   Recover 12 34
ggplot(data = fluH7N9_china_2013_gather, aes(x = Date, y = age, color = outcome)) +
geom_point(aes(color = outcome, shape = gender), size = 1.5, alpha = 0.6) +
geom_path(aes(group = case_id)) +
facet_wrap( ~ province, ncol = 2) +
scale_shape_manual(values = c(15, 16, 17)) +
scale_color_brewer(palette="Set1", na.value = "grey50") +
scale_fill_brewer(palette="Set1") +
labs(
color = "Outcome",
shape = "Gender",
x = "Date in 2013",
y = "Age",
title = "Time from onset of flu to outcome"
)
## Warning: Removed 149 rows containing missing values (geom_point).
## Warning: Removed 98 rows containing missing values (geom_path).

fluH7N9_china_2013_gather_2 <- fluH7N9_china_2013_gather[, -4] %>% gather(group_2, value, gender:province)
## Warning: attributes are not identical across measure variables;
## they will be dropped
fluH7N9_china_2013_gather_2$value <-
mapvalues(fluH7N9_china_2013_gather_2$value, from = c("m", "f", "unknown","Other"), to = c("Male", "Female", "Unknown gender", "Other province"))

fluH7N9_china_2013_gather_2$value <-
factor(fluH7N9_china_2013_gather_2$value, levels = c("Female", "Male","Unknown gender", "Jiangsu", "Shanghai", "Zhejiang", "Other province"))

str(fluH7N9_china_2013_gather_2)
## 'data.frame':    816 obs. of  6 variables:
##  $ case_id: chr  "case_1" "case_2" "case_3" "case_4" ...
##  $ outcome: Factor w/ 2 levels "Death","Recover": 1 1 1 NA 2 1 1 1 NA 1 ...
##  $ Group  : Factor w/ 3 levels "Date of onset",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Date   : Date, format: "2013-02-19" "2013-02-27" ...
##  $ group_2: chr  "gender" "gender" "gender" "gender" ...
##  $ value  : Factor w/ 7 levels "Female","Male",..: 2 2 1 1 1 1 2 2 2 2 ...
p1 <- ggplot(data = fluH7N9_china_2013_gather_2, aes(x = value, fill = outcome, color = outcome)) +
geom_bar(position = "dodge", alpha = 0.7, size = 1) +
scale_fill_brewer(palette="Set1", na.value = "grey50") +
scale_color_brewer(palette="Set1", na.value = "grey50") +
labs(
color = "",
fill = "",
x = "",
y = "Count",
title = "2013 Influenza A H7N9 cases in China")

p1

p2 <- ggplot(data = fluH7N9_china_2013_gather, aes(x = age, fill = outcome, color = outcome)) +
geom_density(alpha = 0.3, size = 1) +
geom_rug() +
scale_color_brewer(palette="Set1", na.value = "grey50") +
scale_fill_brewer(palette="Set1", na.value = "grey50") +
labs(
color = "",
fill = "",
x = "Age",
y = "Density",
title = "Age distribution of flu cases"
)

p2
## Warning: Removed 6 rows containing non-finite values (stat_density).

long to wide or wide to long

library(reshape2)
## Warning: package 'reshape2' was built under R version 3.4.4
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
data("airquality")
names(airquality) <- tolower(names(airquality))
head(airquality)
##   ozone solar.r wind temp month day
## 1    41     190  7.4   67     5   1
## 2    36     118  8.0   72     5   2
## 3    12     149 12.6   74     5   3
## 4    18     313 11.5   62     5   4
## 5    NA      NA 14.3   56     5   5
## 6    28      NA 14.9   66     5   6
aqm <- melt(airquality, id=c("month", "day"), na.rm=TRUE)
head(aqm)
##   month day variable value
## 1     5   1    ozone    41
## 2     5   2    ozone    36
## 3     5   3    ozone    12
## 4     5   4    ozone    18
## 6     5   6    ozone    28
## 7     5   7    ozone    23
newdata <- dcast(aqm, month ~ variable, mean, margins = c("month", "variable"))

抽取特徵

data(fluH7N9_china_2013)
dataset <- fluH7N9_china_2013

dataset $age[which(dataset$age == "?")] <- NA
dataset$age <- as.numeric(as.character(dataset$age))

dataset$case_id <- paste("case", dataset$case_id, sep = "_")

head(dataset)
##   case_id date_of_onset date_of_hospitalisation date_of_outcome outcome
## 1  case_1    2013-02-19                    <NA>      2013-03-04   Death
## 2  case_2    2013-02-27              2013-03-03      2013-03-10   Death
## 3  case_3    2013-03-09              2013-03-19      2013-04-09   Death
## 4  case_4    2013-03-19              2013-03-27            <NA>    <NA>
## 5  case_5    2013-03-19              2013-03-30      2013-05-15 Recover
## 6  case_6    2013-03-21              2013-03-28      2013-04-26   Death
##   gender age province
## 1      m  87 Shanghai
## 2      m  27 Shanghai
## 3      f  35    Anhui
## 4      f  45  Jiangsu
## 5      f  48  Jiangsu
## 6      f  32  Jiangsu
dataset$hospital <- as.factor(ifelse(is.na(dataset$date_of_hospitalisation), 0, 1))

dataset$gender_f <- as.factor(ifelse(dataset$gender == "f", 1, 0))


dataset$province_Jiangsu <- as.factor(ifelse(dataset$province == "Jiangsu", 1, 0))
dataset$province_Shanghai <- as.factor(ifelse(dataset$province == "Shanghai", 1, 0))
dataset$province_Zhejiang <- as.factor(ifelse(dataset$province == "Zhejiang", 1, 0))

dataset$province_other <- as.factor(ifelse(dataset$province == "Zhejiang" | dataset$province
== "Jiangsu" | dataset$province == "Shanghai", 0, 1))

head(dataset)
##   case_id date_of_onset date_of_hospitalisation date_of_outcome outcome
## 1  case_1    2013-02-19                    <NA>      2013-03-04   Death
## 2  case_2    2013-02-27              2013-03-03      2013-03-10   Death
## 3  case_3    2013-03-09              2013-03-19      2013-04-09   Death
## 4  case_4    2013-03-19              2013-03-27            <NA>    <NA>
## 5  case_5    2013-03-19              2013-03-30      2013-05-15 Recover
## 6  case_6    2013-03-21              2013-03-28      2013-04-26   Death
##   gender age province hospital gender_f province_Jiangsu province_Shanghai
## 1      m  87 Shanghai        0        0                0                 1
## 2      m  27 Shanghai        1        0                0                 1
## 3      f  35    Anhui        1        1                0                 0
## 4      f  45  Jiangsu        1        1                1                 0
## 5      f  48  Jiangsu        1        1                1                 0
## 6      f  32  Jiangsu        1        1                1                 0
##   province_Zhejiang province_other
## 1                 0              0
## 2                 0              0
## 3                 0              1
## 4                 0              0
## 5                 0              0
## 6                 0              0
#as.Date(as.character(dataset$date_of_outcome), format = "%Y-%m-%d") - as.Date(as.character(dataset$date_of_onset), format = "%Y-%m-%d")

dataset$days_onset_to_outcome <- as.numeric(as.character(gsub(" days", "",
as.Date(as.character(dataset$date_of_outcome), format = "%Y-%m-%d") -
as.Date(as.character(dataset$date_of_onset), format = "%Y-%m-%d"))))


dataset$days_onset_to_hospital <- as.numeric(as.character(gsub(" days", "",
as.Date(as.character(dataset$date_of_hospitalisation), format = "%Y-%m-%d") -
as.Date(as.character(dataset$date_of_onset), format = "%Y-%m-%d"))))

summary(dataset$date_of_onset)
##         Min.      1st Qu.       Median         Mean      3rd Qu. 
## "2013-02-19" "2013-03-31" "2013-04-08" "2013-04-07" "2013-04-14" 
##         Max.         NA's 
## "2013-07-27"         "10"
dataset$early_onset <- as.factor(ifelse(dataset$date_of_onset < summary(dataset$date_of_onset)[[3]], 1, 0))

dataset$early_outcome <- as.factor(ifelse(dataset$date_of_outcome < summary(dataset$date_of_outcome)[[3]], 1, 0))

#dataset


dataset <- dataset[,c('case_id','outcome'
, 'age', 'hospital','gender_f'
,'province_Jiangsu'
,'province_Shanghai'
,'province_Zhejiang'
,'province_other'
,'days_onset_to_outcome'
,'days_onset_to_hospital'
,'early_onset'
,'early_outcome')]

head(dataset)
##   case_id outcome age hospital gender_f province_Jiangsu province_Shanghai
## 1  case_1   Death  87        0        0                0                 1
## 2  case_2   Death  27        1        0                0                 1
## 3  case_3   Death  35        1        1                0                 0
## 4  case_4    <NA>  45        1        1                1                 0
## 5  case_5 Recover  48        1        1                1                 0
## 6  case_6   Death  32        1        1                1                 0
##   province_Zhejiang province_other days_onset_to_outcome
## 1                 0              0                    13
## 2                 0              0                    11
## 3                 0              1                    31
## 4                 0              0                    NA
## 5                 0              0                    57
## 6                 0              0                    36
##   days_onset_to_hospital early_onset early_outcome
## 1                     NA           1             1
## 2                      4           1             1
## 3                     10           1             1
## 4                      8           1          <NA>
## 5                     11           1             0
## 6                      7           1             1
sum(is.na(dataset))
## [1] 277
library(Amelia)
## Warning: package 'Amelia' was built under R version 3.4.4
## Loading required package: Rcpp
## Warning: package 'Rcpp' was built under R version 3.4.4
## ## 
## ## Amelia II: Multiple Imputation
## ## (Version 1.7.5, built: 2018-05-07)
## ## Copyright (C) 2005-2018 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ##
missmap(dataset)

library(mice)
## Warning: package 'mice' was built under R version 3.4.4
## Loading required package: lattice
## 
## Attaching package: 'mice'
## The following object is masked from 'package:tidyr':
## 
##     complete
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
dataset_impute <- mice(data = dataset[, -2], print = FALSE)
## Warning: Number of logged events: 151
#dataset_impute

missmap(complete(dataset_impute, "long"))

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
#dataset[,c(1,2)]
#right_join(dataset[,c(1,2)], )
datasets_complete <- right_join(dataset[, c(1, 2)],
complete(dataset_impute, "long"),
by = "case_id") %>% select(-.id)

missmap(datasets_complete)

建立模型

train_index <- which(is.na(datasets_complete$outcome))

nrow(datasets_complete)
## [1] 680
datasets_complete$.imp           <- NULL
datasets_complete$province_other <- NULL
length(train_index)
## [1] 285
train_data <- datasets_complete[-train_index, ]
test_data  <- datasets_complete[train_index, -2]

head(train_data)
##   case_id outcome age hospital gender_f province_Jiangsu province_Shanghai
## 1  case_1   Death  87        0        0                0                 1
## 2  case_2   Death  27        1        0                0                 1
## 3  case_3   Death  35        1        1                0                 0
## 5  case_5 Recover  48        1        1                1                 0
## 6  case_6   Death  32        1        1                1                 0
## 7  case_7   Death  83        1        0                1                 0
##   province_Zhejiang days_onset_to_outcome days_onset_to_hospital
## 1                 0                    13                      4
## 2                 0                    11                      4
## 3                 0                    31                     10
## 5                 0                    57                     11
## 6                 0                    36                      7
## 7                 0                    20                      9
##   early_onset early_outcome
## 1           1             1
## 2           1             1
## 3           1             1
## 5           1             0
## 6           1             1
## 7           1             1
set.seed(42)
idx <- sample.int(2 , nrow(train_data), p = c(0.7, 0.3), replace=TRUE)
val_train_data <- train_data[idx == 1, ]
val_test_data  <- train_data[idx == 2, ]

dim(val_train_data)
## [1] 273  12
dim(val_test_data)
## [1] 122  12
#val_train_data[, -1]


library(rpart)
fit <- rpart(outcome ~., data = val_train_data[,-1])
plot(fit, margin = 0.1)
text(fit)

library(rpart.plot)
## Warning: package 'rpart.plot' was built under R version 3.4.4
rpart.plot(fit)

predicted <- predict(fit, val_test_data, type = 'class')

table(val_test_data$outcome, predicted)
##          predicted
##           Death Recover
##   Death      45       6
##   Recover    10      61
sum(val_test_data$outcome == predicted) / length(predicted)
## [1] 0.8688525

評估特徵重要性

library(caret)
## Warning: package 'caret' was built under R version 3.4.4
control <- trainControl(method = "repeatedcv", number = 10, repeats = 3)

set.seed(42)
model1 <- train(outcome ~ ., data = val_train_data[,-1], method = "rpart",preProcess = NULL, trControl = control)

importance1 <- varImp(model1, scale=TRUE)
plot(importance1)

Logistic Regression

fit <- glm(outcome ~ ., data = val_train_data[,-1], family = 'binomial')

summary(fit)
## 
## Call:
## glm(formula = outcome ~ ., family = "binomial", data = val_train_data[, 
##     -1])
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4192  -0.8580   0.4280   0.8088   1.7569  
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             4.126283   0.710054   5.811 6.20e-09 ***
## age                    -0.038345   0.008062  -4.756 1.97e-06 ***
## hospital1              -0.178303   0.329025  -0.542    0.588    
## gender_f1              -0.415759   0.343030  -1.212    0.226    
## province_Jiangsu1       0.069791   0.539090   0.129    0.897    
## province_Shanghai1     -0.034750   0.448794  -0.077    0.938    
## province_Zhejiang1     -0.053870   0.502262  -0.107    0.915    
## days_onset_to_outcome  -0.016929   0.014190  -1.193    0.233    
## days_onset_to_hospital  0.030334   0.040251   0.754    0.451    
## early_onset1           -0.035302   0.425105  -0.083    0.934    
## early_outcome1         -2.252748   0.446894  -5.041 4.63e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 367.30  on 272  degrees of freedom
## Residual deviance: 280.85  on 262  degrees of freedom
## AIC: 302.85
## 
## Number of Fisher Scoring iterations: 5
fit.step <- step(fit)
## Start:  AIC=302.85
## outcome ~ age + hospital + gender_f + province_Jiangsu + province_Shanghai + 
##     province_Zhejiang + days_onset_to_outcome + days_onset_to_hospital + 
##     early_onset + early_outcome
## 
##                          Df Deviance    AIC
## - province_Shanghai       1   280.86 300.86
## - early_onset             1   280.86 300.86
## - province_Zhejiang       1   280.86 300.86
## - province_Jiangsu        1   280.87 300.87
## - hospital                1   281.14 301.14
## - days_onset_to_hospital  1   281.41 301.41
## - days_onset_to_outcome   1   282.23 302.23
## - gender_f                1   282.32 302.32
## <none>                        280.85 302.85
## - age                     1   308.54 328.54
## - early_outcome           1   309.54 329.54
## 
## Step:  AIC=300.86
## outcome ~ age + hospital + gender_f + province_Jiangsu + province_Zhejiang + 
##     days_onset_to_outcome + days_onset_to_hospital + early_onset + 
##     early_outcome
## 
##                          Df Deviance    AIC
## - province_Zhejiang       1   280.86 298.86
## - early_onset             1   280.87 298.87
## - province_Jiangsu        1   280.89 298.89
## - hospital                1   281.15 299.15
## - days_onset_to_hospital  1   281.42 299.42
## - days_onset_to_outcome   1   282.23 300.23
## - gender_f                1   282.32 300.32
## <none>                        280.86 300.86
## - early_outcome           1   309.54 327.54
## - age                     1   309.55 327.55
## 
## Step:  AIC=298.86
## outcome ~ age + hospital + gender_f + province_Jiangsu + days_onset_to_outcome + 
##     days_onset_to_hospital + early_onset + early_outcome
## 
##                          Df Deviance    AIC
## - early_onset             1   280.88 296.88
## - province_Jiangsu        1   280.91 296.91
## - hospital                1   281.15 297.15
## - days_onset_to_hospital  1   281.49 297.49
## - days_onset_to_outcome   1   282.23 298.23
## - gender_f                1   282.32 298.32
## <none>                        280.86 298.86
## - age                     1   309.55 325.55
## - early_outcome           1   309.57 325.57
## 
## Step:  AIC=296.88
## outcome ~ age + hospital + gender_f + province_Jiangsu + days_onset_to_outcome + 
##     days_onset_to_hospital + early_outcome
## 
##                          Df Deviance    AIC
## - province_Jiangsu        1   280.93 294.93
## - hospital                1   281.18 295.18
## - days_onset_to_hospital  1   281.56 295.56
## - gender_f                1   282.35 296.35
## - days_onset_to_outcome   1   282.76 296.76
## <none>                        280.88 296.88
## - age                     1   309.92 323.92
## - early_outcome           1   327.25 341.25
## 
## Step:  AIC=294.93
## outcome ~ age + hospital + gender_f + days_onset_to_outcome + 
##     days_onset_to_hospital + early_outcome
## 
##                          Df Deviance    AIC
## - hospital                1   281.29 293.29
## - days_onset_to_hospital  1   281.62 293.62
## - gender_f                1   282.37 294.37
## - days_onset_to_outcome   1   282.78 294.78
## <none>                        280.93 294.93
## - age                     1   311.92 323.92
## - early_outcome           1   327.45 339.45
## 
## Step:  AIC=293.29
## outcome ~ age + gender_f + days_onset_to_outcome + days_onset_to_hospital + 
##     early_outcome
## 
##                          Df Deviance    AIC
## - days_onset_to_hospital  1   282.04 292.04
## - gender_f                1   282.79 292.79
## <none>                        281.29 293.29
## - days_onset_to_outcome   1   283.83 293.83
## - age                     1   312.10 322.10
## - early_outcome           1   330.23 340.23
## 
## Step:  AIC=292.04
## outcome ~ age + gender_f + days_onset_to_outcome + early_outcome
## 
##                         Df Deviance    AIC
## - gender_f               1   283.21 291.21
## <none>                       282.04 292.04
## - days_onset_to_outcome  1   284.07 292.07
## - age                    1   312.25 320.25
## - early_outcome          1   331.66 339.66
## 
## Step:  AIC=291.21
## outcome ~ age + days_onset_to_outcome + early_outcome
## 
##                         Df Deviance    AIC
## <none>                       283.21 291.21
## - days_onset_to_outcome  1   285.87 291.87
## - age                    1   312.48 318.48
## - early_outcome          1   333.25 339.25
#val_test_data$outcome
predict(fit.step, val_test_data) >0
##     1     2     5     8    12    14    17    20    21    26    29    30 
## FALSE  TRUE  TRUE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE 
##    36    43    44    45    51    57    64    72    73    75    81    83 
## FALSE FALSE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE FALSE FALSE  TRUE 
##    98   106   111   119   131   143   155   156   157   161   170   171 
## FALSE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE FALSE FALSE  TRUE FALSE  TRUE 
##   172   173   182   185   195   208   210   211   213   214   218   219 
## FALSE  TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE 
##   225   230   247   252   253   257   261   264   269   275   289   296 
##  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE 
##   298   299   301   302   307   308   309   315   321   331   332   333 
## FALSE FALSE  TRUE  TRUE  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE FALSE 
##   337   347   361   370   388   391   405   406   409   425   432   438 
##  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE  TRUE 
##   443   451   453   454   458   463   481   489   490   502   527   529 
##  TRUE FALSE  TRUE FALSE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE 
##   537   543   547   549   556   558   561   562   563   573   579   588 
##  TRUE  TRUE FALSE  TRUE FALSE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE 
##   590   593   601   604   618   622   631   635   638   651   667   669 
## FALSE  TRUE FALSE  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE 
##   673   679 
##  TRUE  TRUE
predicted <- ifelse(predict(fit.step, val_test_data) >0, 1, 0 )
predicted <- factor(predicted, levels = c(0,1), labels = c('Recover', 'Death'))
table(val_test_data$outcome, predicted)
##          predicted
##           Recover Death
##   Death        31    20
##   Recover      15    56

RandomForest

library(rpart)
fit1 <- rpart(outcome ~., data = val_train_data[,-1])
predicted1 <- predict(fit1, val_test_data, type = 'class')
table(val_test_data$outcome, predicted1)
##          predicted1
##           Death Recover
##   Death      45       6
##   Recover    10      61
library(randomForest)
## Warning: package 'randomForest' was built under R version 3.4.4
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:dplyr':
## 
##     combine
## The following object is masked from 'package:ggplot2':
## 
##     margin
fit2 <- randomForest(outcome ~., data = val_train_data[,-1], ntree = 100)
predicted2 <- predict(fit2, val_test_data, type = 'class')
table(val_test_data$outcome, predicted2)
##          predicted2
##           Death Recover
##   Death      49       2
##   Recover     4      67
library(e1071)
## Warning: package 'e1071' was built under R version 3.4.4
fit3 <- svm(outcome ~., data = val_train_data[,-1])
predicted3 <- predict(fit3, val_test_data, type = 'class')
table(val_test_data$outcome, predicted3)
##          predicted3
##           Death Recover
##   Death      36      15
##   Recover     7      64
#predict(fit1, val_test_data)

control <- trainControl(method = "repeatedcv", number = 10, repeats = 3)
set.seed(42)
model2 <- train(outcome ~ ., data = val_train_data[,-1], method = "rf", preProcess = NULL, trControl = control)
importance2 <- varImp(model2, scale=TRUE)
plot(importance2)

使用ROC Curve 比較模型

library(rpart)
fit1 <- rpart(outcome ~., data = val_train_data[,-1])

library(ROCR)
## Warning: package 'ROCR' was built under R version 3.4.4
## Loading required package: gplots
## Warning: package 'gplots' was built under R version 3.4.4
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
predictions1 <- predict(fit1, val_test_data, type="prob")
pred.to.roc1 <- predictions1[, 2]
pred.rocr1 <- prediction(pred.to.roc1, as.factor(val_test_data$outcome))
perf.rocr1 <- performance(pred.rocr1, measure = "auc", x.measure = "cutoff")
perf.tpr.rocr1 <- performance(pred.rocr1, "tpr","fpr")



library(randomForest)
fit2 <- randomForest(outcome ~., data = val_train_data[,-1], ntree = 100 )

predictions2 <- predict(fit2, val_test_data, type="prob")
pred.to.roc2 <- predictions2[, 2]
pred.rocr2 <- prediction(pred.to.roc2, as.factor(val_test_data$outcome))
perf.rocr2 <- performance(pred.rocr2, measure = "auc", x.measure = "cutoff")
perf.tpr.rocr2 <- performance(pred.rocr2, "tpr","fpr")


library(e1071)
fit3 <- svm(outcome ~., data = val_train_data[,-1],probability = TRUE)
?svm
## starting httpd help server ...
##  done
predictions3 <- predict(fit3, val_test_data, probability=TRUE)
pred.to.roc3 <- attr(predictions3, "probabilities")[, 2]
pred.rocr3 <- prediction(pred.to.roc3, as.factor(val_test_data$outcome))
perf.rocr3 <- performance(pred.rocr3, measure = "auc", x.measure = "cutoff")
perf.tpr.rocr3 <- performance(pred.rocr3, "tpr","fpr")

plot(perf.tpr.rocr1,main='ROC Curve', col=1)
legend(0.7, 0.2, c('rpart','randomforest', 'svm'), 1:3)
plot(perf.tpr.rocr2, col=2, add=TRUE)
plot(perf.tpr.rocr3, col=3, add=TRUE)