SIR загвар

Халдварт өвчний дэгдэлтийн SIR загвар нь зөвхөн өвчинд мэдрэмтгий буюу susceptible, халдвар авсан буюу infected, эдгэрсэн буюу recovered гэсэн гурав төлөв бүхий энгийн загвар юм.

\begin{align*}
\frac{dS}{dt} &= -\beta S I\\
\frac{dI}{dt} &= \beta S I - \gamma I\\
\frac{dR}{dt} &= \gamma I
\end{align*}

Энд \(S\), \(I\), \(R\) үсгээр харгалзан susceptible буюу өвчинд мэдрэмтгий, infected буюу халдвар авсан, recovered буюу эдгэрсэн гэсэн бүлгүүдийн тоон хэмжээг тэмдэглэсэн. Харин \(t\) нь цаг хугацааг илэрхийлнэ. \(\beta\) болон \(\gamma\) нь харгалзан халварлах болон эдгэрэх эрчим бөгөөд тус загварын параметрүүд юм.

Мөн төрөлт болон жам ёсны нас баралт, дархлаажуулалт буюу вакцинжуулалт зэргийг тусгасан загварууд байдаг. Харин COVID-19 вирусын хувьд одоогоор вакцин гараагүй байгааг анхаарч улмаар дэгдэлтийнх нь үед хүн амын тоо их өөрчлөгдөхгүй гэж үзвэл SIR загварыг хэрэглэх боломжтой.

Тус загварын хувьд хүн амын тоо яагаад тогтмол байдгийг тайлбарлая. Хүн бүр \(S\), \(I\), \(R\) төлвийн аль нэгд байх нь гарцаагүй тул нийт хүн амын тоо \(N=S+I+R\) болно. Нөгөө талаас SIR загварыг илэрхийлсэн гурван тэгшитгэлүүдийг нэмбэл уламжлалуудын нийлбэр тэгтэй тэнцүү болно. \[\frac{dS}{dt}+\frac{dI}{dt}+\frac{dR}{dt}=-\beta S I+\beta S I - \gamma I+\gamma I=0\]Гэтэл \(\frac{dN}{dt}=\frac{dS}{dt}+\frac{dI}{dt}+\frac{dR}{dt}\) байх тул дээрх чанар ёсоор \(\frac{dN}{dt}=0\) болно. Харин тогтмол хэмжигдэхүүний уламжлал л тэгтэй тэнцүү байдаг тул \(N\) буюу хүн амын тоо тогтмол болно.

Загвар дифференциал тэгшитгэлээр тодорхойлогдсон динамик систем хэлбэртэй тул системийн тогтвортой байдлыг шинжлэх зайлшгүй шаардлагатай. Ингээд \((S,I)\) фазын хавтгай байгуулав.

Фазын хавтгайг харвал систем асимптот тогтвортой ажээ.

SIS загвар

SIR загварын хувьд эдгэрсэн хүн ахин өвчлөхгүй буюу дархлаа тогтдог гэж үздэгийг анхаарах хэрэгтэй. Өөрөөр хэлбэл халдварт өвчний үед зөвхөн \[S\rightarrow I\rightarrow R\]процесс явагдана гэж тооцдог. Харин эдгэрсэн хүмүүс ахин халдвар авах буюу \[S\rightarrow I\rightarrow R\rightarrow S\]байдал их бол үүнийг нь тооцсон дараах SIS загвар авч үзнэ.

\begin{align*}
\frac{dS}{dt} &= - \beta S I + \gamma I\\
\frac{dI}{dt} &= \beta S I - \gamma I
\end{align*}

Энэхүү SIS загвар нь SIR загвартай адил \(\beta\), \(\gamma\) гэсэн параметрүүдтэй. SIS загварын фазын хавтгайг дор харуулав.

Фазын хавтгайг харвал систем асимптот тогтвортой байна. Гэхдээ өвчний дэгдэлт бүрэн намжихгүй буюу байнгын тогтвортой өвчлөл байсаар байж болох дүр зураг ажиглагдаж байна.

Загварын параметрийн үнэлэлт

Халдварт өвчний дэгдэлтийн симуляц хийх эсвэл прогноз гаргахын тулд \(\beta\) болон \(\gamma\) параметрүүдийг өгөх шаардлагатай. Иймд бидний нэн тэргүүний зорилт бол эдгээр параметрүүдийг үнэлэх явдал юм.

Загварын параметрүүдийг дундаж квадрат алдааг минимумчлах байдлаар оптимзацын аргаар үнэлэх боломжтой. Загварын параметрүүдийг үнэлэхийн тулд халдварын тоо хэмжээг хугацааны нэг ижил алхамд бүртгэсэн гэе. Тухайлбал COVID-19 цар тахлын хувьд халдвар авсан хүний тоог өдөр өдрөөр бүртгэснээ улс орон бүр мэдээлж байгаа билээ. Иймд тус вирусээр халдварласан хүний тоог бүртгэсэн \(I_1,I_2,\ldots,I_n\) хэлбэртэй өгөгдөл байгаа гэе. Тэгвэл загварын тусламжтай өдөр өдрийн халдвартай хүний тоог үнэлж болно. Тийнхүү үнэлсэн утгуудыг \(\hat I_1,\hat I_2,\ldots,\hat I_n\) гэж тэмдэглэе. Тэгвэл үнэлэлтийн алдааг дараах байдлаар хэмжиж болно. \[RSS=\sum_{i=1}^{n}(I_i-\hat I_i)^2\]Улмаар алдааг хамгийн бага байлгах \(\beta\) болон \(\gamma\) параметрүүдийн утгыг олох ёстой. Өөрөөр хэлбэл \(RSS(\beta,\gamma)\rightarrow\text{min}\) бодлого бодно. Эндээс олдох шийд глобал минимумын цэг байхын тулд \(RSS(\beta,\gamma)\) функц гүдгэр буюу график нь тогоо шиг хэлбэртэй байх ёстой. SIR загварын талаар уншиж судалж байхад уг функцийн график тогоо шиг хэлбэртэй байхыг харуулсан туршилтын үр дүн бүхий веб хуудас RPubs.com дээр байсантай тааралдаж байлаа.

R програм ашиглаж хийсэн загварын үнэлгээ болон COVID-19 цар тахлын прогноз

SIR загварыг ажиллуулахын тулд хүн амын тоо болон өвчтэй хүний тоо хэрэгтэй. Хүн амын тоог www.worldometers.info сайтаас авсан. Харин COVID-19 вирусын тархалтын мэдээг бүх улс орны хувьд өдөр өдрөөр тооцож гаргасан өгөгдөл олсон. Иймд SIR загварыг хүссэн улсынхаа хувьд бас дэлхийн хэмжээнд ажиллуулах боломжтойг энд тэмдэглэе. Ингээд дэлхийн улс орнуудаас жишээ болгон Япон, Өмнөд Солонгос улсуудыг авч үзье.

country <- "Japan"; N <- 126400000
country <- "Korea, South"; N <- 51200000

COVID-19 вирусын нөхцөл байдлын мэдээг www.kaggle.com сайтаас татан авсан. Ингээд R програм дээр бичиж ажиллуулсан кодоо хэсэг хэсгээр нь товч тайлбартайгаар оруулъя.

Шаардлагатай багц суулгах

SIR загварыг тодорхойлох дифференциал тэгшитгэлүүдийн системийг бодоход R програмын deSolve багц ашигласан.

install.packages("deSolve")

Өгөгдөл оруулах

COVID-19-ийн мэдээллийг дээр дурдсан веб сайтаас татан авч хадгалсанаа дараах байдлаар ачаалав.

infected <- read.csv(
  file = "~/Downloads/time_series_covid19_confirmed_global.csv",
  check.names = FALSE)
deaths <- read.csv(
  file = "~/Downloads/time_series_covid19_deaths_global.csv",
  check.names = FALSE)
recovered <- read.csv(
  file = "~/Downloads/time_series_covid19_recovered_global.csv",
  check.names = FALSE)

Өгөгдөл ялгаж авах

Сонгон авсан улсынхаа мэдээллийг дараах байдлаар ялгаж авч улмаар вектор хэлбэртэй болгож хувиргав.

infected <- as.vector(t(infected[which(infected$`Country/Region` == country),-{1:4}]))
deaths <- as.vector(t(deaths[which(deaths$`Country/Region` == country),-{1:4}]))
recovered <- as.vector(t(recovered[which(recovered$`Country/Region` == country),-{1:4}]))

Ялгаж авсан өгөгдлөө ашиглаж вирусын халдвар, эдгэрэлт болон нас баралтын түүх харуулсан диаграмм байгуулав.

matplot(
  x = 1:length(infected),
  y = cbind(infected, deaths, recovered),
  main = paste0("COVID-19\n", country), xlab = "Day", ylab = "Infected, Deaths, Recovered",
  type = "l", lwd = 1, lty = 1, bty = "l", col = c("blue", "red", "green"))
legend(
  x = 0, y = max(infected, deaths, recovered), legend = c("Infected", "Deaths", "Recovered"),
  pch = 1, col = c("blue", "red", "green"), bty = "n")

Байгуулсан диаграммын нэгийг жишээ болгон харуулъя.

Өгөгдөл бэлдэх

Дээрх өгөгдөлд халдвар, эдгэрэлт болон нас баралтыг хуримтлуулан бүртгэсэн байгаа тул тухайн өдөрт харгалзах халдвартай хүний тоог гаргахын тулд дараах үйлдэл хийв.

Infected <- infected - recovered - deaths

Ингээд Infected хувьсагч манай үндсэн өгөгдөл болно.

Өгөгдөл цэвэрлэх

Дээрх өгөгдөл 2020 оны 2 сарын 21-нээс эхэлж буй бөгөөд тэр үед өвчлөл илрээгүй байсан улсын мэдээг 0 гэж оруулсан байгаа тул Infected хувьсагчийн эхэнд байгаа тэгүүдийг зайлуулна.

if (any(Infected == 0)) {
  Infected <- Infected[-{1:max(which(Infected == 0))}]
}

SIR загварыг тодорхойлох дифференциал тэгшитгэлийн систем

Дифференциал тэгшитгэлийн системийг бодохдоо deSolve багцын ode() функц ашиглах тул тус функцийн шаардлагад нийцэх хэлбэрээр уг системийг оруулах ёстой. Үүний тулд дараах код бичсэн.

SIR <- function(time, state, parameters) {
  par <- as.list(c(state, parameters))
  with(par, {
    dS <- -beta / N * I * S
    dI <- beta / N * I * S - gamma * I
    dR <- gamma * I
    list(c(dS, dI, dR))
  })
}

Дээрх функцийн time болон state хувьсагчаар ode() функц утга дамжуулах ба харин parameters нь SIR загварын \(\beta\) болон \(\gamma\) параметрүүдийг дамжуулахаар нэмж бичсэн аргумент юм.

Загварын параметрийн үнэлгээ

Одоо SIR загварын \(\beta\) болон \(\gamma\) параметрүүдийг үнэлье. Үүний тулд эхлээд дээр тодорхойлсон \(RSS\) хэмжигдэхүүнд харгалзах функц зарлах хэрэгтэй.

RSS <- function(parameters) {
  names(parameters) <- c("beta", "gamma")
  out <- deSolve::ode(y = init, times = time, func = SIR, parms = parameters)
  fit <- out[,3]
  sum((Infected - fit)^2)
}

Дээрх функц өгөгдөл ба загвар хоёрын зөрүүг олж гаргаж өгнө. Мэдээж энэ зөрүү загварын SIR загварын \(\beta\) болон \(\gamma\) параметрүүдээс хамаарах тул эдгээр параметрүүдийн утгыг уг функцийн аргумент болгон дамжуулна.

Ингээд загварын параметрээ үнэлье. Үүний тулд дээр хэлж байсанчлан RSS() функцийг минимумчилна. Кодын эхэнд хугацааг зааж өгч байна. Бид хугацааг өдөр өдрөөр тооцсон. Яагаад гэвэл COVID-19 вирусын тархалтын мэдээ өдөр өдрөөр өгөгдсөн байгаа билээ. Удаах мөрд SIR загвар буюу бидний ашиглаж буй динамик системийн анхны нөхцөл харагдаж байна. Үүнийг хугацааны эхэнд байсан халдвартай хүний тоогоор авна. Харин эдгэрсэн өвчтөн байхгүй байх тул R = 0 гэж өгнө. Гурав дугаар мөрд RSS() функцийг минимумчлах тушаал харагдаж байна. Үүнд R програмын optim() функц ашиглав. Үүний цаана Nelder-Mead-ийн алгоритм ажиллана. Уг алгоритм нь олон хувьсагчтай, зааглалтгүй, шугаман бус оптимизацын бодлогод зориулсан арга юм. Мөн тус алгоритм градиент буюу уламжлал ашигладаггүй. Тодруулбал энэхүү алгоритм нь симплексэд тулгуурласан шууд хайлтын арга юм.

time <- 1:length(Infected)
init <- c(S = N - Infected[1], I = Infected[1], R = 0)
fit <- optim(par = c(0.5, 0.5), fn = RSS)
parameter <- setNames(fit$par, c("beta", "gamma"))

Прогноз

SIR загварын \(\beta\) болон \(\gamma\) параметрүүдийг үнэлсэн тул одоо загвараа ажиллуулж COVID-19 буюу шинэ Корона вирусын цар тахлын дэгдэлтийг прогнозлох боломжтой. Гэхдээ загвар бодит байдалтай хэр таарч тохирох вэ, параметрийн үнэлгээ хэр сайн болсон бэ гэдгийг анхаарах шаардлагатай. Иймд прогнозыг өнгөрсөн үеийн түүх буюу одоо байгаа мэдээлэлтэйгээ харьцуулах шаардлагатай. Ийнхүү прогноз гаргаад халдварын ажиглагдсан болон үнэлэгдсэн утгуудыг харьцуулсан диаграмм байгуулахын тулд дараах кодыг бичив.

t <- 1:{max(time) + 90} # прогноз хийх хугацаа, өдрөөр
forecast <- as.data.frame(deSolve::ode(y = init, times = t, func = SIR, parms = parameter))

plot(
  x = 1:length(infected),
  y = infected - recovered - deaths,
  main = paste0("COVID-19\n", country), xlab = "Day", ylab = "Infected",
  type = "l", lwd = 1, lty = 1, bty = "l", col = c("blue"),
  xlim = c(1,max(t)), ylim = c(0,max(Infected, forecast$I)))
lines(x = forecast$I, col = "red")
legend(
  x = 0, y = max(Infected, forecast$I), legend = c("Observed", "Estimated"),
  pch = 1, col = c("blue", "red"), bty = "n")

Ингээд жишээ болгон авч үзэж буй Япон болон Өмнөд Солонгос улсуудын хувьд гарсан үр дүнг дор харуулав.

SIR загвар Японы өгөгдөлтэй маш сайн тохирч байна. Хэрэв нөхцөл байдал энэ хэвээрээ байх аваас ойрын үед Япон дахь дэгдэлт намжих ажээ.

Харин Өмнөд Солонгос улсын хувьд загвар бодит байдлаас гажуудсан байна. Бодит байдлыг харуулсан цэнхэр динамикийг харвал тус улс дахь дэгдэлт нэлээн даамжрах хэдий ч намжихаар төлөвтэй ажээ. Дэгдэлт ийнхүү даамжирч буй явдал нь ЭМЯ-наас байсхийгээд л мэдээлээд байгаачлан ганц нэгхэн хүнээс болж томхон кластер халдвар гараад байгаатай холбоотой биз ээ. Мөн ийнхүү загвараас хазайж байгаа байдал нь ЭМЯ-ны мэдээнд олонтаа анхааруулж хэлээд буй давалгаа гээч нь ч байж болох юм.

Эцэст нь "Өөртөө биш бусдад илүү хайртай хүн маск зүүдэг" гэдгийг сануулъя. Тийнхүү би таны амь, амьдарлыг хайрлаж байна. Та ч бас минийхийг хайрлаарай.