1 코로나19 데이터 분석의 시작

코로나 바이러스 전파현황을 시각화해 현재 코로나 바이러스 현황을 정리해 나 자신과 주변사람에게 알리는 것이 충분한 의미가 있다고 생각했다.
코드 자체는 3월 말에 만들어두었지만, 최대한 깔끔한 폼에 갖추어 정리해 데이터를 이해하고 현상을 해석해보았다.
데이콘 등에서 코로나 시각화 대회를 진행해 의미있는 시각화 결과물들을 많이있다. 뿐만 아니라 현재도 많은 서비스들이 코로나와 관련한 현황을 전달하고 있다.
이 분석은 주로 시간 흐름에 따른 코로나 발생현황에 초점을 맞추고 있다.



## 1.1 라이브러리 불러오기

library(dplyr)
library(ggplot2)
library(lubridate)
library(plotly)
library(reshape2)
library(gmodels)

데이터 조작을 하는 dplyr, 시각화를 위한 ggplot2, 날짜를 다루는 lubridate,
자료를 인터랙티브하게 만들어주는 plotly, 데이터 구조를 바꿀 때 쓰는 reshape2 패키지를 불러온다.

1.2 데이터 불러오기

patientinfo<-read.csv("Patientinfo.csv")
time<-read.csv("Time.csv")
time_age<-read.csv("TimeAge.csv")
policy<-read.csv("Policy.csv")

공개된 여러 데이터가 더 있지만, 위의 4개 데이터만을 활용해 분석해보고자 한다
데이터의 출처 아래와 같다.
https://www.kaggle.com/kimjihoo/coronavirusdataset

1.3 데이터 타입 변경 및 새로운 변수 생성 

time$date<-as.Date(time$date)

time<- time %>%
        mutate(confirmed_rate= confirmed/test,
               released_rate = released/confirmed,
               inpatient_number = confirmed-released-deceased)

datebreak<-seq(as.Date("2020-01-20"), as.Date("2020-05-31"), by="1 week")

str(patientinfo) time데이터는 날짜, 시간, 테스트, 음성, 양성, 격리해제, 사망자의 수 정보를 가지고 있다.
해당 데이터를 이용해 확진자의 비율, 격리해제자의 비율, 확진자의 수에서 격리해제자와 사망자의 수를 빼어 실제환자 수 변수를 만들어내었다.
또 datebreak에 데이터가 시작하는 날짜와 끝나는 날까지를 1주일 간격으로 나타냈다. datebreak를 통해 x축에 날짜가 1주일 간격으로 보여질 수 있도록 했다.

2 시각화


## 2.1 테스트 대비 양성비율

confirmed_rate<-time %>%
  ggplot(aes(x=date, y=confirmed_rate))+geom_line(color="red")+geom_point()+scale_x_date(breaks=datebreak)+ylim(c(0.0,0.10))+
  theme(axis.text.x = element_text(angle=30, hjust=1))+ggtitle("테스트 대비 양성비율")

ggplotly(confirmed_rate)

위 그래프를 통해 우리는 테스트 대비 확진비율이 3월 초이후부터 감소하고 있음을 알 수 있다. 5월 31일의 경우 테스트대비 양성비율은 1.25%다. 

tail(time)
##           date time   test negative confirmed released deceased confirmed_rate
## 128 2020-05-26    0 839475   806206     11225    10275      269     0.01337145
## 129 2020-05-27    0 852876   820550     11265    10295      269     0.01320825
## 130 2020-05-28    0 868666   834952     11344    10340      269     0.01305910
## 131 2020-05-29    0 885120   849161     11402    10363      269     0.01288187
## 132 2020-05-30    0 902901   865162     11441    10398      269     0.01267138
## 133 2020-05-31    0 910822   876060     11468    10405      270     0.01259082
##     released_rate inpatient_number
## 128     0.9153675              681
## 129     0.9138926              701
## 130     0.9114951              735
## 131     0.9088756              770
## 132     0.9088366              774
## 133     0.9073073              793

5월 31일까지 누적 코로나 검사 수는 91만 건, 이 중 음성은 87만건, 양성은 1만 1천 건임을 확인할 수 있다.

2.2 일별 확진자 수

daily_added_confirmed<-diff(time$confirmed)
daily_added_confirmed<-append(daily_added_confirmed,1,after = 0)
daily_added_confirmed<-as.numeric(daily_added_confirmed)
time<-time %>%
  mutate(daily_added_confirmed = daily_added_confirmed)

daily<-time %>%
  ggplot(aes(x=date, y=daily_added_confirmed))+geom_point(size="0.5")+geom_line(color="red")+ggtitle("일별 신규 확진자 수")+
  scale_x_date(breaks=datebreak)+theme(axis.text.x = element_text(angle=30, hjust=1))+
  geom_vline(xintercept = as.numeric(time$date[c(63,107)]), color="blue", linetype=2)+
  geom_vline(xintercept = as.numeric(time$date[77]), color="purple", linetype=3)


ggplotly(daily)

일별 확진자 수와 사회적 거리두기 기간을 비교해보았다.
사회적 거리두기는 3월 22일 시작해서 5월 5일 종료되어 생활 속 거리두기로 권고수칙이 바뀌었다. 이는 policy 데이터를 통해 확인할 수 있다.
파란색 점선으로 시작일과 종료일을 표시했고, 보라색 선을 사회적 거리두기가 시작일로부터, 코로나 바이러스의 잠복기로 알려진 2주 후가 지난 날을 표시한 것이다.
사회적 거리두기 이후 2주동안의 일별 확진자 수는 잠복기간의 영향이 있어서인지 사회적 거리두기 2주 전과 비슷한 100여명 전후로 비슷한 수준이다.
그러나 잠복기간이 끝난 이후부터 사회적 거리두기가 종료될 때까지 일별 확진자 수가 50명 아래로 감소하며 사회적 거리두기의 정부 권고수칙이 효과가 있는 것으로 보인다.
추가 확진자 수가 제일 적은 날은 증가폭이 4명에 그쳤다. 코로나 확진자 수가 한자리수로 줄어들 때는 적어도 대한민국에서는 코로나가 끝날 것 같은 생각이었다.


## 2.3 일별 누적 확진자 수

daily_accumu_confirmed<-time %>%
  ggplot(aes(x=date, y= confirmed, group=1))+geom_point(size="0.5")+geom_line(color="red")+ggtitle("일별 누적 확진자 수")+
  scale_x_date(breaks=datebreak)+theme(axis.text.x = element_text(angle=30, hjust=1))
ggplotly(daily_accumu_confirmed)

누적 확진자 수를 그래프로 그려보면 위와 같다. 로지스틱 회귀선 같은 느낌이다. 이제 어느정도 특정값을 향해 수렴하면 좋겠다.

2.4 누적확진자 대비 누적완치자 비율

cured_rate<-time %>%
  ggplot(aes(x=date, y=released_rate, group=1))+geom_point()+geom_line(color="red")+ggtitle("누적확진자 대비 누적완치자 비율") +
  scale_x_date(breaks=datebreak) +theme(axis.text.x = element_text(angle=30, hjust=1)) +
  scale_y_continuous(limits= c(0.40,1))

ggplotly(cured_rate)

누적 확진자 대비 누적완치자의 비율을 나타냈다. 40%이상일 때부터 나타내었는데 3월 25일 그 비율이 40%를 넘어 처음으로
그 점이 나타난다. 5월 31일자 기준 그 비율은 90.7%이다. 약 9.3%가량은 여전히 치료중이다.

2.5 누적확진자, 누적완치자, 누적 사망자 추이

time_reshape<- melt(time, id.vars = "date", measure.vars = c("confirmed","released","deceased"))

ggplotly(time_reshape %>%
  ggplot(aes(x=date, y= value, col=variable))+geom_point()+geom_line()+ggtitle("누적 확진자, 누적 완치자, 누적 사망자 추이")+
  scale_x_date(breaks=datebreak) +theme(axis.text.x = element_text(angle=30, hjust=1))+ylab(""))

누적 확진자, 완치자, 사망자의 수를 그래프로 그렸다. 5월 중순부터 확진자 수의 증가폭이 다시 오름새를 보이고 있다.
반면 완치된(격리해제된) 수의 증가폭은 상대적으로 완만해보이는 것을 확인할 수 있다.

2.6 코로나19 실질환자 수 추이

#치료중인 인원 치료인원=누적 확진자 수- 누적 완치자-누적 사망자
real_patient_number<-time %>%
  ggplot(aes(x=date, y=inpatient_number, group=1))+geom_point()+geom_line(color="red")+ggtitle("코로나19 실질환자 수 추이")+
  scale_x_date(breaks=datebreak)+theme(axis.text.x = element_text(angle=30, hjust=1))+ylab("실질 환자 수")

ggplotly(real_patient_number)

국내 코로나 실질환자 수의 추이를 확인해보자. 실질 환자 수=확진자 수-격리해제자 수-사망자 수로 구했다.
실질 환자 수는 5월 26일 681명까지 감소했다가 다시 증가하고 있음을 알 수 있다.

3 연령대별 코로나 현황

3.1 연령대별 누적 확진자 추이

time_age$date<-as.Date(time_age$date)
age_accumu_confirmed<-time_age %>%
  ggplot(aes(x=date, y=confirmed, color=age, group=age))+geom_line()+geom_point()+ggtitle("연령대별 누적 확진자 추이")+
scale_x_date(breaks=datebreak) +theme(axis.text.x = element_text(angle=30, hjust=1))+ylab("누적 확진자 수")

ggplotly(age_accumu_confirmed)

연령대 별 누적 확진자 수의 추이를 보면 20대, 50대, 40대 순으로 많다.
신천지교회에서의 집단감염으로 20대가 많이 확진되었지만 그 이후에도 확진자 증가 수는 20대가 가장 많다.

3.2 연령대별 누적 사망자 추이

age_accumu_deceased<-time_age %>%
  ggplot(aes(x=date, y= deceased, color=age, group=age))+geom_line()+geom_point()+ggtitle("연령대별 누적 사망자 추이")+
  scale_x_date(breaks=datebreak) +theme(axis.text.x = element_text(angle=30, hjust=1))+ylab("누적 사망자 수")


ggplotly(age_accumu_deceased)

가장 많은 확진자 규모를 가진 20대의 사망자는 아직없다. 반면 두번째로 확진자 수가 적은 80대의 경우
사망자 수가 5월 31일 기준 131명으로 가장많다. 80대의 확진자 수인 498명 대비 사망자 수는 131명으로
국내 80대의 코로나 치명율은 26%에 이른다. 면역력이 상대적으로 약한 고령자에게 치명적이라 해석할 수 있다.
2번째로 사망자 수가 많은 70대 군의 경우 5월 31일 기준 확진자 수는 725명, 사망자 수는 80명으로
그 치명률은 11%다.

3.3 연령대별 코로나19 치명률

age_fatality_rate<-time_age %>%
  ggplot(aes(x=date, y= deceased/confirmed, color=age))+geom_line()+geom_point()+ggtitle("연령대별 코로나19 치명률")+
  scale_x_date(breaks=datebreak) +theme(axis.text.x = element_text(angle=30, hjust=1))
ggplotly(age_fatality_rate)

이 비율을 시각화하면 위와 같은 그래프로 표현할 수 있다. 60대이상의 고령군의 사망률의 변화 추이를 확인 할 수 있다. 

str(patientinfo)
## 'data.frame':    4004 obs. of  18 variables:
##  $ patient_id        : num  1e+09 1e+09 1e+09 1e+09 1e+09 ...
##  $ global_num        : int  2 5 6 7 9 10 11 13 19 21 ...
##  $ sex               : Factor w/ 3 levels "","female","male": 3 3 3 3 2 2 3 3 3 2 ...
##  $ birth_year        : int  1964 1987 1964 1991 1992 1966 1995 1992 1983 1960 ...
##  $ age               : Factor w/ 13 levels "","0s","100s",..: 9 7 9 5 5 9 5 5 7 10 ...
##  $ country           : Factor w/ 15 levels "Bangladesh","Canada",..: 8 8 8 8 8 8 8 8 8 8 ...
##  $ province          : Factor w/ 17 levels "Busan","Chungcheongbuk-do",..: 16 16 16 16 16 16 16 16 16 16 ...
##  $ city              : Factor w/ 157 levels "","Andong-si",..: 41 93 91 94 122 91 91 33 129 122 ...
##  $ disease           : logi  NA NA NA NA NA NA ...
##  $ infection_case    : Factor w/ 30 levels "","Bonghwa Pureun Nursing Home",..: 23 23 5 23 5 5 5 23 23 5 ...
##  $ infection_order   : int  1 1 2 1 2 3 3 1 2 3 ...
##  $ infected_by       : Factor w/ 415 levels "","#11185","#11243",..: 1 1 296 1 8 9 9 1 1 9 ...
##  $ contact_number    : Factor w/ 84 levels "","-","0","1",..: 76 44 24 83 29 58 3 3 71 67 ...
##  $ symptom_onset_date: Factor w/ 100 levels ""," ","2020-01-19",..: 4 1 1 5 1 1 1 1 1 1 ...
##  $ confirmed_date    : Factor w/ 120 levels "","2020-01-20",..: 3 6 6 6 7 7 7 9 11 11 ...
##  $ released_date     : Factor w/ 106 levels "","2020-02-05",..: 2 19 10 7 13 10 5 13 12 17 ...
##  $ deceased_date     : Factor w/ 39 levels "","2020-02-19",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ state             : Factor w/ 3 levels "deceased","isolated",..: 3 3 3 3 3 3 3 3 3 3 ...

4 완치자의 치료기간 알아보기

4.1 연령대별 치료기간

released_dataset<- patientinfo %>%
  filter(state=="released") 

released_dataset$released_date<-as.Date(released_dataset$released_date)
released_dataset$confirmed_date<-as.Date(released_dataset$confirmed_date)
released_dataset<-released_dataset %>%
  mutate(period=released_date-confirmed_date)
released_dataset$period<-as.numeric(released_dataset$period)

released_dataset<-released_dataset %>%
  filter(!is.na(period))

age_cure_period <- released_dataset %>%
  ggplot(aes(x=age, y=period))+geom_boxplot()+geom_point(alpha=0.15)

ggplotly(age_cure_period)

patientinfo데이터에서 완치자들의 데이터만을 선택해 그들의 회복기간을 구했다.회복기간=완치일-양성판정일로 정의했다
위의 표를 통해 각 연령대별 회복기간의 최소값, 1사분위 수, 중간값, 3사분위 수, 최댓값을 확인할 수 있다.
70대, 80대에서의 치료기간이 다른 연령대의 치료기간보다 더 길다는 것을 알 수 있다.

 released_dataset %>%
  ggplot(aes(x=age, y=period))+geom_violin(trim = FALSE)

바이올린 그래프로 치료기간을 보면 분포의 퍼진정도를 더 쉽게 이해할 수 있다. 60대,70대,80대의 경우 치료기간의 분포가
다른 연령층에 비해 더 크고, 치료기간 자체가 더 긴 경향을 알 수 있다. 다시말해 표준편차가 크다.

5 치료기간이 나이에 따라 유의미한 차이가 있을까

5.1 차이 알아보기

aggregate(data=released_dataset, period~age, mean)
##     age   period
## 1       13.44444
## 2    0s 20.68750
## 3  100s 66.00000
## 4   10s 20.67213
## 5   20s 22.88976
## 6   30s 22.62944
## 7   40s 23.33621
## 8   50s 23.69580
## 9   60s 26.41104
## 10  70s 30.89333
## 11  80s 35.29310
## 12  90s 28.78571
aggregate(data=released_dataset, period~age, sd)
##     age    period
## 1        6.984109
## 2    0s  9.477825
## 3  100s        NA
## 4   10s  9.076933
## 5   20s 10.821199
## 6   30s 11.279372
## 7   40s 11.099662
## 8   50s 11.487374
## 9   60s 11.971020
## 10  70s 16.097317
## 11  80s 17.844365
## 12  90s 11.060523

20대의 치료기간과 60대의 치료기간의 차이를 알아보고자 한다.
20대 평균 치료기간은 22.88일 표준편차는 10.82일이다.60대의 경우 평균치료기간은 26.41일 표준편차는 11.97일이다.
언뜻보기에 평균치료기간과 표준편차가 있어보이는데 과연 유의미한 차이인지 알아보자.

age_20_age_60<- released_dataset %>%
  filter(age=="20s"|age=="60s")
var.test(age_20_age_60[age_20_age_60$age=="20s",19], age_20_age_60[age_20_age_60$age=="60s",19])
## 
##  F test to compare two variances
## 
## data:  age_20_age_60[age_20_age_60$age == "20s", 19] and age_20_age_60[age_20_age_60$age == "60s", 19]
## F = 0.81712, num df = 380, denom df = 162, p-value = 0.1192
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.625188 1.053298
## sample estimates:
## ratio of variances 
##          0.8171249

20대와 60대의 등분산성 차이를 검정했다. p-value가 0.11로 유의수준인 0.05보다 크므로 등분산성을 만족한다.
따라서 두 집단은 분산의 차이가 없다.

t.test(period~age, data=age_20_age_60, var.equal=TRUE, conf.level=0.95)
## 
##  Two Sample t-test
## 
## data:  period by age
## t = -3.3661, df = 542, p-value = 0.0008166
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.576214 -1.466345
## sample estimates:
## mean in group 20s mean in group 60s 
##          22.88976          26.41104

20대와 60대의 치료기간을 비교한 결과 p-value는 0.0008로 유의수준인 0.05보다 낮아
20대의 평균치료기간과 60대의 평균치료기간간의 통계적으로 유의미한 차이가 있다고 할 수 있다.

## 5.2 회귀식 만들기

summary(lm(period~birth_year, data=released_dataset))
## 
## Call:
## lm(formula = period ~ birth_year, data = released_dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.508  -7.994  -2.076   5.648  57.851 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 318.5978    34.1709   9.324   <2e-16 ***
## birth_year   -0.1492     0.0173  -8.625   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.79 on 1217 degrees of freedom
##   (274 observations deleted due to missingness)
## Multiple R-squared:  0.0576, Adjusted R-squared:  0.05683 
## F-statistic: 74.39 on 1 and 1217 DF,  p-value: < 2.2e-16

출생연도에 따라 치료기간의 차이가 있는지를 알아보기 위한 회귀식을 만들어 확인해보자. 결정계수는 약 5.6%, p-value는 일반적인 유의수준인 0.05보다 낮다.

  plot(period~birth_year, data=released_dataset)

cor(released_dataset$period, released_dataset$birth_year, use = 'complete.obs',method = 'pearson')
## [1] -0.2400079

상관계수는 -0.24 태어난 연도와 입원기간은 약한 음의 상관관계가 있음을 알 수 있다.
다시 말해 나이와 입원기간은 약한 양의 상관관계에 있다.

cor.test(released_dataset$period, released_dataset$birth_year)
## 
##  Pearson's product-moment correlation
## 
## data:  released_dataset$period and released_dataset$birth_year
## t = -8.6249, df = 1217, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2922168 -0.1863727
## sample estimates:
##        cor 
## -0.2400079


5.3 한계

우리는 20대와 60대의 치료기간이 유의미하게 차이가 있나 확인을 했다.
그리고 나이와 입원기간이 약한 양의 상관관계가 있다고 봤다.
이 데이터의 한계는 관측치가 편향되어 있을 수 있다는 점을 꼽을 수 있다.
환자 정보를 다룬 patientinfo의 관측치는 4004개이다. 현재 대한민국의 누적 확진자는 1만 1천여명을 넘어간다.
전체 환자 중에서 1/3이 약간 넘는 수준이다. 안타깝게도 patientinfo는 전체 확진자의 정보를 담지 못하고 있다.
그렇기 때문에 어떠한 기준을 가지고 patientinfo에 확진자정보가 담겨있는지 알 수 없다. 따라서 편향된 기준이 있을 수 있다.

6 성별과 환자의 상태간의 연관성 알아보기

성별(남/녀)과 환자의 상태(사망/격리/완치)간의 연관성을 알아보기 위해 카이제곱검정을 하려고 한다.
카이제곱검정에서
귀무가설은 성별은 환자의 상태와 연관성이 없다.
대립가설은 성별은 환자의 상태와 연관성이 있다.
로 해석할 수 있다.
gmodels라이브러리에서 crosstable함수를 이용해 카이제곱검정을 할 수 있다.
성별변수에 결측치가 있어, 결측치를 제외한 남성과 여성만 선택을 해준다.

test<-patientinfo %>%
  filter(sex=="female"|sex=="male") 

CrossTable(test$sex, test$state, expected = TRUE, chisq = TRUE)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |              Expected N |
## | Chi-square contribution |
## |           N / Row Total |
## |           N / Col Total |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  3674 
## 
##  
##              | test$state 
##     test$sex |  deceased |  isolated |  released | Row Total | 
## -------------|-----------|-----------|-----------|-----------|
##       female |        27 |       714 |      1280 |      2021 | 
##              |    40.706 |   736.559 |  1243.735 |           | 
##              |     4.615 |     0.691 |     1.057 |           | 
##              |     0.013 |     0.353 |     0.633 |     0.550 | 
##              |     0.365 |     0.533 |     0.566 |           | 
##              |     0.007 |     0.194 |     0.348 |           | 
## -------------|-----------|-----------|-----------|-----------|
##         male |        47 |       625 |       981 |      1653 | 
##              |    33.294 |   602.441 |  1017.265 |           | 
##              |     5.642 |     0.845 |     1.293 |           | 
##              |     0.028 |     0.378 |     0.593 |     0.450 | 
##              |     0.635 |     0.467 |     0.434 |           | 
##              |     0.013 |     0.170 |     0.267 |           | 
## -------------|-----------|-----------|-----------|-----------|
## Column Total |        74 |      1339 |      2261 |      3674 | 
##              |     0.020 |     0.364 |     0.615 |           | 
## -------------|-----------|-----------|-----------|-----------|
## 
##  
## Statistics for All Table Factors
## 
## 
## Pearson's Chi-squared test 
## ------------------------------------------------------------
## Chi^2 =  14.14328     d.f. =  2     p =  0.0008488399 
## 
## 
##

결과를 보면 p-value는 0.0008로 유의수준인 0.05보다 충분히 작아 귀무가설을 기각할 수 있다.
따라서 성별과 환자의 상태는 연관이 있다고 해석할 수 있다.
이 차이는 환자의 상태 3가지(사망, 격리, 완치) 모두에서 차이가 있다는 것은 아니고 1가지 상태에서만 차이가 있을 수도 있다.
어느 부분에서 차이가 있는지 알기 위해서는 추가검정을 시행해야한다.
(본페로니 사후검정을 통해 알 수 있다는 사실을 구글링과 질문을 통해 알았다. 해당 방식은 범주형 변수를 검정할 때 사용해 딱 적절하다고 이해했는데, 에러코드가 나서 원인 파악중이다.)

7 결론

시간흐름에 따라 데이터를 분석한 결과, 강력한 사회적 거리두기 정책은 효과적인 것으로 보인다.
하지만 생활 속 거리두기로 그 정책이 바뀐 이후, 매일 두자릿수 이상의 확진자가 발생한다.
확진자 대비 격리해제자의 비율은 90% 정도 수준에 수렴하고 있다.
연령군(20대와 60대)에 따라 치료기간의 차이를 알아보니 유의미한 차이가 있음을 알 수 있었다.
하지만 연령은 치료기간에 대해 큰 설명변수가 아닌 것도 알 수 있었다. 연령대는 치료기간에 대해 약 5%의 설명력을 가지고 있었다.