Chapter-2.Rmd

# Döntési fa és véletlen erdő {#Chapter-2}

```{css, echo=FALSE}
p {
  text-align: justify;
}
```

Az egyik talán legfontosabb kérdés a járványügyben, hogy mennyire halálos a vírus. Már a megjelenésekor nagyon hamar ismertté vált a tény: **elsősorban az idősekre és krónikus betegre jelent kockázatot**. Ennek a ténynek empirikus teszteléséhez olyan adattáblára van szükségünk, melyben rögzítésre kerülnek a megfertőzöttek demográfiai adatai és az eset végső kimenete (ezt a fajta adattípust nevezik "line list"-nek). 

## Adatok bemutatása

### Forrás

A fentebb leírt célra lett létrehozva egy nyílt projekt^[Open COVID-19 Data Working Group, *Detailed Epidemiological Data from the COVID-19 Outbreak*, http://virological.org/t/epidemiological-data-from-the-ncov-2019-outbreak-early-descriptions-from-publicly-available-data/337,  letöltve: 2020.12.28.], amely során ezeket az adatokat több országra gyűjtik ki és teszi online elérhetővé.

```{r}
dat <- vroom::vroom('latestdata.csv')
```

Az adattábla összesen `r nrow(dat)` megfigyelést tartalmaz 33 változóval. Azonban a változók jelentős része használhatatlan elemzésre és jelentős mértékben van szükség az adatok tisztítására is. Egyetlen példa említéseként a kimenetben nem egyszerűen *died* vagy *elhunyt* szerepelt, hanem annak számos szinonímája is, amelett, hogy sok esetben nincs elérhető adat.

```{r fig.height=5}
set.seed(1)
dat %>% 
  count(outcome) %>% 
  na.omit() %>% 
  {mutate(., s = runif(n = nrow(.)))} %>%
  ggplot(aes(label = outcome, size = s)) +
  geom_text_wordcloud(color =  'aquamarine4') +
  labs(title = "Az nyers adattáblában a kimenet változó által felvett értékek") +
  theme_minimal()
```

### Adatok elemzéshez való előkészítése

Az adatok megtiszítása előtt el kell döntenünk, hogy mely változókkal is érdemes foglalkozni, mivel sajnos számos változó esetében hiányzik rengeteg sor a táblában.

```{r}
dat %>% 
  apply(2, function(x) sum(is.na(x))) %>% 
  {data.frame(var = names(.), nobs = nrow(dat)-., robs = (nrow(dat)-.)/nrow(dat))} %>% 
  set_names('var', 'Hiánytalan adatok száma', 'Hiánytalanok aránya') %>% 
  gt(
    rowname_col = "var"
  ) %>% 
  fmt_percent(
    columns = vars('Hiánytalanok aránya'),
    dec_mark = ','
  ) %>% 
  tab_header(title = 'Adattáblában lévő hiánytalan megfigyelések aránya')
```

Az adatok hiányosságának és felhasználhatóságának figyelembevételével az alábbi változókat érdemes bevonni a modellbe: \
- **életkor**\
- **nem**\
- **ország**\
- **van-e krónikus betegsége**\
- **fertőződöttség kimutatásának ideje**\
- **földrajzi szélesség**\
- **földrajzi hosszúság**\
- **megbetegedés kimenete**\

A változók tisztítása során egyik legfontosabb, hogy a kimenet oszlopban összevontuk minden elhalálozásra megfelelő szinonímát egységesen az elhunyt kategóriába, és minden felépültnek megfelelőt a felépültbe. Mivel jelen elemzés kutatási kérdése, hogy mi a valószínűsége, hogy valaki túléli, így elvetettünk minden megfigyelést, amely esetében az alany még betegség alatt áll, vagy a kimenet ismeretlen. Az ország oszlopból kontinensként új magyarázó változót határoztunk meg, és az évek esetében ismeretlennek vetük azokat, ahol intervallum került megadásra^[Az intervallumok hossza nem egyezett meg.].

```{r data-prep}
dat %>% 
  select(outcome, age, sex, country, chronic_disease_binary, latitude, longitude,
         date_confirmation) %>% 
  mutate(
    continent = countrycode::countrycode(country, origin = 'country.name',
                                         destination = 'continent'),
    date = as.numeric(difftime(time1 = lubridate::dmy(date_confirmation), 
                               time2 = lubridate::dmy("01/01/2020"), 
                               units = 'days')),
    age = as.numeric(ifelse(str_detect(age, '-'), NA, age))
  ) %>% 
  mutate( # cleaning the depedent var
    outcome = case_when(
      outcome == 'death' ~ 'died',
      outcome == 'died' ~ 'died', 
      outcome == 'Death' ~ 'died',
      outcome == 'dead' ~ 'died',
      outcome == 'Dead' ~ 'died',
      outcome == 'Died' ~ 'died',
      outcome == 'Deceased' ~ 'died',
      outcome == 'discharge' ~ 'survived',
      outcome == 'discharged' ~ 'survived',
      outcome == 'Discharged' ~ 'survived',
      outcome == 'Discharged from hospital' ~ 'survived',
      outcome == 'recovered' ~ 'survived',
      outcome == 'released from quarantine' ~ 'survived',
      outcome == 'recovered' ~ 'survived',
      outcome == 'Recovered' ~ 'survived',
      outcome == 'Recovered' ~ 'survived',
      outcome == 'Recovered' ~ 'survived',
      T ~ 'NA'
    )
  ) %>% 
  filter(outcome != 'NA') %>% 
  mutate(
    died = outcome == 'died'
  ) %>% select(-outcome, -date_confirmation) %>% 
  {dat <<- .}
```

Az így kapott új adattáblában már csupán `r nrow(dat)` megfigyelés és `r ncol(dat)` változó szerepel.

```{r}
dat %>% 
  head %>%
  gt %>% 
  tab_header(title = 'Adattábla modellezésre előkészítve')
```



```{r}
dat %>% 
  select(age, sex, died) %>% 
  GGally::ggpairs(mapping = aes(color = died), 
                  title = 'Életkor, nem és az elhalálozás kapcsolata')
```

```{r}
dat %>% 
  select(continent, date, died) %>% 
  GGally::ggpairs(aes(color = died),
                  title = 'A kontinens a rögzítés dátuma és az elhalálozás kapcsolata')
```

Az ábrákból sajnos gyorsan kiderül, hogy a kimeneti adatra való rászűrést követően az adattáblában már jórészt csak ázsiai országok képviseltetik magukat. Más kontintensekről származó adatok esetében ez az oszlop nem került megfelelően dokumentálásra.

```{r}
dat %>% 
  group_by(continent, country) %>% 
  summarise(n = n()) %>% 
  ungroup() %>% 
  na.omit() %>% 
  mutate(
    r = n/nrow(dat)
  ) %>% 
  arrange(continent) %>% 
  set_names('continent', 'country', 'n', 'n/sum(n)') %>% 
  gt(
    rowname_col = 'country',
    groupname_col = 'continent'
  ) %>% 
  fmt_percent(
    columns = 'n/sum(n)',
  ) %>% 
  summary_rows(groups = T, columns = vars(n), fns = list(TOTAL = 'sum')) %>% 
  summary_rows(groups = T, columns = vars('n/sum(n)'), 
               fns = list(TOTAL = 'sum'), formatter = fmt_percent) %>%
  tab_options(
    summary_row.background.color = "#ACEACE",
    row_group.background.color = "#FFEFDB"
  ) %>% 
  tab_header(title = 'Megfigyelések száma országonként', 
             subtitle = 'Elsődleges adattisztítás után')
```

A táblázatból jól látható, hogy a legtöbb felhasználható megfigyelésünk Indiából származik. Logikusnak tűnik egy dummy változót létrehozni ezért Indiára (india = {TRUE, ha Indiából származik, FALSE egyébként}). Indiára való kontrolálás elméleti megfontolásból is helyesnek hangzik, figyelembe véve, hogy az ország rendkívül szegény, a városok zsúfoltak, orvosokból és felszerelésből, pedig hiány van. Korábban a spanyol nátha idején is kiugróan magas volt Indiában a halálozás^[https://www.economist.com/asia/2020/03/21/if-covid-19-takes-hold-in-india-the-toll-will-be-grim].

```{r add india}
dat %>% 
  mutate(
    india = country == 'India'
  ) %>% 
  select(-country, -continent, -latitude, -longitude) %>% 
  mutate_if(is.character, as.factor) %>% 
  na.omit %>% 
  {dat <<- .}
```

```{r}
dat %>% 
  head() %>% 
  gt %>% 
  tab_header(title = 'Végső adattábla modellezéshez')
```

## Modellbecslés

```{r rpart1, fig.height=7}
rpart1 <- rpart::rpart(dat, formula = as.factor(died) ~ ., cp = .01) 

rattle::fancyRpartPlot(rpart1, palettes = 'OrRd', 
                       sub = 'cp = .01', main = '1. döntési fa') 
```

```{r rpart2, fig.height=7}
rpart2 <- rpart::rpart(dat, formula = as.factor(died) ~ ., cp = .003)

rattle::fancyRpartPlot(rpart2, palettes = 'OrRd', # plot
                       sub = 'cp = .003', main = '2. döntési fa') 
```

A döntési fa alapján még mindig az életkor a legfontosabb paraméter, ugyanakkor láthatjuk, hogy az indiai származás, a dátum és a krónikus betegség is megjelenik, mint meghatározó töréspont. A dátum esetében az értelmezés az alábbi: január 1-e veszi fel az egyes értéket és így az évek számát fejezi ki.


## Keresztvalidáció és véletlen erdő hiperparamétere

```{r}
library(caret)
keresztval <- trainControl(method = "cv", number = 10)
ism_keresztval <- trainControl(method = "repeatedcv", number=10, repeats=20)
cp <- data.frame(cp=seq(from=0, by=0.01, to=0.2))

train(as.factor(died)~., data=na.omit(dat), method="rpart", trControl=keresztval)
train(as.factor(died)~., data=na.omit(dat), method="rpart", trControl=ism_keresztval)

trained_hyperparam <- as.double(train(as.factor(died)~., data=na.omit(dat), method="rpart", trControl=ism_keresztval)$bestTune["cp"])

train(as.factor(died)~., data=na.omit(dat), method="rpart", trControl=keresztval, tuneGrid=cp)
```

A keresztvalidáció alapján a véletlen erdő hiperparaméterét kalibráljuk. Mind a sima, mind az ismétléses módszerrel nagyjából `r trained_hyperparam`
a kalibrációs paraméter értéke.

```{r}
library(randomForest)
```

```{r}
erdo <- randomForest(factor(died) ~ ., data = dat, ntree=500, cp = .069)
```


```{r}
plot(erdo)
varImpPlot(erdo)
```

A modell hibája nagyjából 100 fa után stabilizálódik, így az 500 ismétlést éreztük helyénvalónak.
Ez alapján pedig a legfontosabb paraméter az év, azonban a korábbi döntési fán nem megtalálható
dátum változó is fontos a kimenetel szempontjából.


```{r}
pred_tree <- predict(rpart1, na.omit(dat), type="class")
pred_tree_table <- table(na.omit(dat)$died, pred_tree)

pred_tree_table

pred_forest <- predict(erdo, na.omit(dat), type="class")

pred_forest_table <- table(na.omit(dat)$died, pred_forest)

pred_forest_table

confusionMatrix(pred_tree_table)$overall[c("Accuracy", "Kappa")]
confusionMatrix(pred_forest_table)$overall[c("Accuracy", "Kappa")]

confusionMatrix(pred_tree_table)$byClass[c("Sensitivity", "Specificity")]

confusionMatrix(pred_forest_table)$byClass[c( "Sensitivity", "Specificity")]

```

A döntési fa és véletlen erdő összehasonlításakor mind a kappa, mind pedig a találati arány
értékét összehasonlítva az erdő lenne a jobb modell.