R로 코로나 예측모델 만들기

changsubchang·2020년 5월 1일
2

이번 학기에는 통계수업을 듣는다.
바이오통계를 전공 / 가르치시던 교수님이 가르치셔서 이번에는 흥미로운 주제로 코로나 모델을 짜는 것을 숙제로 내주셨다. 파이썬도 R도 못하는 쪼렙이지만, 이래저래 검색해가며 코로나 모델을 짜봤는데 간략하게 소개하려고 한다.

코로나 감염 모델

학교에서 푸아송 분포를 배웠다. 푸아송 분포는 아래처럼 생긴 함수로, 분포의 평균과 분산이 동일한 함수다.

엄청 거대하게 위키피디아를 스크린샷해오긴 했는데, 핵심은 평균이랑 분산 = 람다라고 함 = 이랑 같기 때문에, 평균이 커질수록, 보다 넓게 퍼지는 특성을 띈다는 점이다. (분산이 커지니까)

과제의 조건은 다음과 같았다

  1. 감염자가 접촉하는 사람의 수는 푸아송 분포 (람다가 10) 을 따른다.
  2. 접촉하면 감염률이 p=0.2 (접촉한다고 다 감염되는건 아니니까)
  3. 처음에 감염자는 5명
  4. 5일을 하나의 시간 단위 로 !
  5. 감염의 사건은 한 시간단위에서만 유효함. 한 시간 단위가 지나면 원래 감염자는 다다음 시간단위에 영향력을 미치지 않음 (죽거나, 감염력이 없어지거나, 낫거나 라고 생각하면 간편)
  6. 각 사례를 100개씩 돌려봄 (즉, 하나의 사례가 100번 돌아가게 해서 평균적으로 어떤지 보려고 함)
  7. 케이스를 3개로 나눔
  • A. 기본 케이스
  • B. 50일째부터 갑자기 마스크 착용 시작. 그래서 감염률이 p=0.14 로 감소
  • C. 60일째부터 갑자기 사회적 거리두기 시작. 그래서 접촉인원인 푸아송 분포의 람다가 3으로 감소

실제로 코딩할때는 R알못이라 고생꽤나 했지만, 완성은 지었고 결과는 아래와 같음

A. 기본케이스

본 케이스에서는, 75일째 (15 unit)이후 급격하게 감염자수가 증가함을 확인할 수 있다. 즉 아무런 Mitigation이 계획되고 수행되지 않는다면, 100일째 370~620만명 정도의 감염자가 발생한다고 말할 수 있다. 다만, 최악의 경우에는 1000만명 이상의 신규 감염자가 발생할 수 있다는 점에서 5명으로 시작한 감염자수가 100일만에 200만배로 급증함을 확인할 수 있다. 본 데이터상에서는 아직 Peak 점에 도달했다고 말하기 어려우며, 동시에 Flatten 되었다고 보기도 어렵다.

options(scipen=10000)
N<-20 #number of timeframes (100 days divided by 5 days each)
D<-100 #number of datasets
OUT<-matrix(0,N,D)
for(a in 1:D){
  K<-5 #those infected in the initial timeline
  for(i in 1:N){
    C<-rpois(1, 10*K) #aggregate of individual poisson distributions
    C_infected<-rbinom(1, C, 0.2) #aggregate of individual bernouille distributions
    K<-C_infected #assigning number of infected for next timeframe
    OUT[i,a]<-c(C_infected)
  }
}
boxplot.matrix(OUT, use.cols=FALSE, main="COVID19 infected cases", xlab="Timeframe, 1unit = 5days", ylab="Infected")
![](https://media.vlpt.us/images/changsubchang/post/c25aaf35-7a8f-4b22-a3f9-ae0b6624ca6c/image.png)

B. 50일째부터 갑자기 마스크 착용 시작. 그래서 감염률이 p=0.14 로 감소

본 케이스에서는 (a)에 비해서 확산세가 더디게 나타남을 확인할 수 있다. 특히 50일째 이후인 51일째부터는 (Time unit 상 11부터) 마스크 착용으로 인해 감염세가 둔화되었다. 즉, 접촉인원은 동일하게 유지되지만 감염자의 비말을 통한 감염 확률이 기존 0.2에서 0.14로 줄어들게 된 것이다. 이를 통해 (a)에서 사태가 급격하게 전환되는 75일째 이전에 조치를 취했으므로 증가폭을 상대적으로 빠른 시일에 줄일 수 있었으며, 결과적으로는 100일째 8~18만명 정도의 감염자로 신규 감염자수를 큰 폭으로 줄일 수 있었다. 단, 본 시뮬레이션에서도 Peak 혹은 Plateau에 도달함을 확인할 수는 없었다.

N<-20
D<-100
OUT<-matrix(0,N,D) # N x 2 container matrix with zero
for(a in 1:D){
  K<-5
  for(i in 1:10){
    C<-rpois(1, 10*K)
    C_infected<-rbinom(1, C, 0.2)
    K<-C_infected
    OUT[i,a]<-c(C_infected)
  }
  for(i in 11:N){
    C<-rpois(1, 10*K)
    C_infected<-rbinom(1, C, 0.14)
    K<-C_infected
    OUT[i,a]<-c(C_infected)
  }
}
boxplot.matrix(OUT, use.cols=FALSE, main="COVID19 infected cases", xlab="Timeframe, 1unit = 5days", ylab="Infected")
![](https://media.vlpt.us/images/changsubchang/post/aa46fe57-4323-4c9c-9f41-29733d81cda7/image.png)

C. 60일째부터 갑자기 사회적 거리두기 시작. 그래서 접촉인원인 푸아송 분포의 람다가 3으로 감소

마지막 시뮬레이션에서는 놀라운 결과를 확인할 수 있다. (b)의 사례에 더하여, 60일째 이후 (Time unit 상 13부터) 절대적인 접촉자수가 평균 10명에서 3으로 줄어든 것이다. 이를 통해 마스크 착용으로 인한 감염률 감소와, 접촉자수가 줄어듦으로 인한 급격한 감염자 감소의 효과를 얻을 수 있었다. 본 사례에서는 특히 Peak 도달 이후 Flatten 에서 더 나아가 심지어 신규 감염자수가 점차 줄어듦을 확인할 수 있다. 특히 사회적 거리두기를 시작한 13 이후에 신규 감염자수가 현저하게 줄어들기 시작함을 확인할 수 있다. 이후 감소세가 지속되어 100일째 되는 시점에는 신규 감염자수가 거의 없게 됨을 확인할 수 있다. 절대적인 접촉자수가 10명에서 3명으로 줄어들게 됨으로써, 마스크 착용대비 현저하게 높은 방역 효과를 달성함을 확인할 수 있었다.

N<-20
D<-100
OUT<-matrix(0,N,D) # N x 2 container matrix with zero
for(a in 1:D){
  K<-5
  for(i in 1:10){
    C<-rpois(1, 10*K)
    C_infected<-rbinom(1, C, 0.2)
    K<-C_infected
    OUT[i,a]<-c(C_infected)
  }
  for(i in 11:12){
    C<-rpois(1, 10*K)
    C_infected<-rbinom(1, C, 0.14)
    K<-C_infected
    OUT[i,a]<-c(C_infected)
  }
  for(i in 13:N){
    C<-rpois(1, 3*K)
    C_infected<-rbinom(1, C, 0.14)
    K<-C_infected
    OUT[i,a]<-c(C_infected)    
  }
}
boxplot.matrix(OUT, use.cols=FALSE, main="COVID19 infected cases", xlab="Timeframe, 1unit = 5days", ylab="Infected")
![](https://media.vlpt.us/images/changsubchang/post/54b85286-e6e1-4cee-b722-bd48719258ff/image.png)

과제의 결과 및 한계

과제 수행의 편의를 위해서 몇몇 요소들로 제한되었지만, 아래와 같은 요소들이 추가적으로 고려된다면 보다 현실적인 시뮬레이션 결과를 얻을 수 있을 것이다.

1. 잠복기에 대한 고려

코로나바이러스는 평균적으로 5일의 잠복기가 있다고 알려져 있으며, 특히 무증상 감염자의 경우 본인도 모르는 와중에 바이러스를 전파할 수 있다. 이에 잠복기를 고려하여, 개별 Timeframe에서 해당 감염의 영향력이 끝나는게 아니라, 특정 감염자의 영향력을 평균 10일(2 Time unit) 을 고려한다면 보다 현실적인 결과가 나올 수 있다.

2. 개별적 인구통계학적 특성에 따른 감염률의 차등적 접근

바이러스 보균자의 나이, 기저질환 유무등에 따라 감염률을 차등적으로 적용할 수 있다. 예측 모델이기 때문에 일반화하여 적용하는 것이 어렵지만, 지역 / 연령 등에 따라 몇몇 인구집단을 Grouping 하여 차등적인 감염률을 적용할 수 있을 것이다. (예. 노인이 많이 거주하는 시골 / 지방 지역, 혹은 병원 등 지역에서 발병하는 경우 차등적인 감염률 적용

3. 접촉인원의 Lambda 현실화

2번과 유사한 맥락에서, 인구 통계학적 특성을 감안한 차등적인 Lambda 적용 또한 고려할 수 있다. 도시의 접촉인원과, 지방의 접촉인원이 상이할 것이며, 한국에서 2월에 집단감염이 일어난 교회 등 작은 공간에서 대규모의 인원이 밀집되는 곳에서의 접촉인원 수는 평균적인 수준을 상회할 것이다. 본 바이러스의 전염성이 강하고, 밀접접촉시 높은 확률로 감염되는 만큼 접촉자수의 예측력을 높이면 보다 정확한 시뮬레이션 결과를 도출할 수 있을 것이다.

4. Dynamic 한 Parameter update

시뮬레이션 모델에 투입된 2개의 Parameter, 감염률과 평균접촉인원인 Lambda를 매일마다 업데이트 하여 모델을 수정한다면, 해당 시점에서 사용 가능한 데이터로 얻을 수 있는 최적의 시뮬레이션 결과를 얻을 수 있을 것이다. 특히, 평균접촉인원은 접촉자의 수에 따라 Exponential 하게 증가할 수 있기 때문에 Dynamic한 업데이트 수행이 필수적으로 요구된다.

할 때는 스트레스 많이받았는데, 생각보다 재밌는 결과가 나와서 흥미로웠던 과제였다. 도움이 되셨길~ 코로나 극복!

profile
데린이임니다

1개의 댓글

comment-user-thumbnail
2020년 5월 2일

재밌네여 잘 읽었습니다.

답글 달기