Нахождение поворотной точки (точки перегиба) инфекций COVID-19 с использованием пакета перегиба - PullRequest
2 голосов
/ 17 февраля 2020

Я видел несколько SO вопросов о расчете точки перегиба. Я все еще не уверен, сделал ли я это правильно. На основании лабораторно подтвержденных кумулятивных данных о случаях в эпицентре текущей эпидемии c мы попытались определить точку перегиба. Я использовал пакет inflection и рассчитал точку перегиба как «08 февраля 2020». Я также попытался вычислить первую и вторую директивы как предполагаемое увеличение каждого и изменение скорости. Я немного разбираюсь в этом по математике, но просто следую примеру из разных SO сообщений. Мой вопрос: согласуются ли эти результаты из следующих графиков? Если нет, как улучшить мой код.

df<-structure(list(date = structure(c(18277, 18278, 18279, 18280, 
18281, 18282, 18283, 18284, 18285, 18286, 18287, 18288, 18289, 
18290, 18291, 18292, 18293, 18294, 18295, 18296, 18297, 18298, 
18299, 18300, 18301, 18302, 18303, 18304, 18305, 18306, 18307), 
class = "Date"), 
cases = c(45, 62, 121, 198, 258, 363, 425, 
        495, 572, 618, 698, 1590, 1905, 2261, 2639, 3125, 4109, 5142, 
        6384, 8351, 10117, 11618, 13603, 14982, 16903, 18454, 19558, 
       20630, 21960, 22961, 23621)), 
class = "data.frame", row.names = c(NA, -31L))
xlb_0<- structure(c(18281, 18285, 18289, 18293, 
                    18297, 18301, 18305,   
                    18309), class = "Date")
library(tidyverse)
# Smooth cumulative cases over time 
df$x = as.numeric(df$date)
fit_1<- loess(cases ~ x, span = 1/3, data = df) 
df$case_sm <-fit_1$fitted 

# use inflection to obtain inflection point
library(inflection)
guai_0 <- check_curve(df$x, df$case_sm)
check_curve(df$x, df$cases)
#> $ctype
#> [1] "convex_concave"
#> 
#> $index
#> [1] 0
guai_1 <- bese(df$x, df$cases, guai_0$index)
structure(guai_1$iplast, class = "Date")
#> [1] "2020-02-08"

# Plot cumulativew numbers of cases 
df %>% 
  ggplot(aes(x = date, y = cases ))+
  geom_line(aes(y = case_sm), color = "red") +
  geom_point() + 
  geom_vline(xintercept = guai_1$iplast) + 
  labs(y = "Cumulative lab confirmed infections")

# Daily new cases (first derivative) and changing rate (second derivative)
df$dt1 = c(0, diff(df$case_sm)/diff(df$x))
fit_2<- loess(dt1 ~ x, span = 1/3, data = df) 
df$change_sm <-fit_2$fitted 
df$dt2 <- c(NA, diff(df$change_sm)/diff(df$x)) 

df %>%  
  ggplot(aes(x = date, y = dt1))+
  geom_line(aes(y = dt1, 
                color = "Estimated number of new cases")) + 
  geom_point(aes(y = dt2*2, color = "Changing rate")) +
  geom_line(aes(y = dt2*2, color = "Changing rate"))+
  geom_vline(xintercept = guai_1$iplast) + 
  labs(y = "Estimatede number of new cases") +
  scale_x_date(breaks = xlb_0, 
               date_labels = "%b%d") + 
  theme(legend.title = element_blank())
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 row(s) containing missing values (geom_path).

Создано в 2020-02-17 пакетом Представлять (v0.3.0)

Ответы [ 2 ]

1 голос
/ 17 февраля 2020

Очень быстрый закулисный график, основанный на ваших данных, показывает

calc_d <- function(x) c(0, diff(x))
df %>%
    mutate(
        first_deriv_cases = calc_d(cases),
        second_deriv_cases = calc_d(calc_d(cases))) %>%
    pivot_longer(-date) %>%
    ggplot(aes(date, value)) +
        geom_line() +
        facet_wrap(~name, scale = "free_y", ncol = 1) +
        geom_smooth()

enter image description here

Итак, точка перегиба 8 февраля соответствует первой производной (т.е. функции плотности), имеющей максимум в этой точке.

1 голос
/ 17 февраля 2020

Я собирался написать комментарий, но я выдвигал ограничение на число символов.

Я не знаком с пакетом inflection, поэтому я не один, чтобы судить, является ли 2020-02-08 истинным перегибом , Однако я скажу, что это трудно ответить с помощью R, потому что R не обязательно хорош в вычислении производных. Если бы у вас было оценочное уравнение линии - тогда вы могли бы потенциально использовать его для построения первой и второй производных. Вычисление грубой дельты с помощью разницы в (Y_n+1-Y_n)/(X_n+1-X_n) никогда не является оптимальным, поскольку теоретически производный представляет собой дельту двух точек, бесконечно близко расположенных друг к другу. Вы принципиально не можете получить отличную оценку производной. Вы даже можете увидеть это, потому что вы вынуждены изменить эту оценку на n или n+1. Кроме того, можно ожидать, что точка перегиба x_0 будет локальным минимумом / максимумом в первой производной и равна нулю во второй производной. Так что я не думаю, что ваш второй заговор помогает. Но это может быть связано с подсчетом дельты.

Что бы я сделал, это сначала подгонил ваши данные к какой-то модели. В этом примере я собираюсь использовать пакет dr4pl, чтобы смоделировать ваши данные для модели logisti c с 4 параметрами. Поскольку функция четырехпараметрической модели хорошо известна, я могу написать, какой должна быть первая и вторая производные функции, затем построить эти значения, используя stat_function в пакете ggplot2.

library(ggplot2)
library(dr4pl)
df<-structure(list(date = structure(c(18277, 18278, 18279, 18280, 
                                      18281, 18282, 18283, 18284, 18285, 18286, 18287, 18288, 18289, 
                                      18290, 18291, 18292, 18293, 18294, 18295, 18296, 18297, 18298, 
                                      18299, 18300, 18301, 18302, 18303, 18304, 18305, 18306, 18307), 
                                    class = "Date"), 
                   cases = c(45, 62, 121, 198, 258, 363, 425, 
                             495, 572, 618, 698, 1590, 1905, 2261, 2639, 3125, 4109, 5142, 
                             6384, 8351, 10117, 11618, 13603, 14982, 16903, 18454, 19558, 
                             20630, 21960, 22961, 23621)), 
              class = "data.frame", row.names = c(NA, -31L))
xlb_0<- structure(c(18281, 18285, 18289, 18293, 
                    18297, 18301, 18305,   
                    18309), class = "Date")

df$dat_as_num <- as.numeric(df$date)

dr4pl_obj <- dr4pl(cases~dat_as_num, data = df, init.parm = c(30000, 18300, 2, 0))
#first derivative derivation
d1_dr4pl <- function(x, theta, scale = F){
  if (any(is.na(theta))) {
    stop("One of the parameter values is NA.")
  }
  if (theta[2] <= 0) {
    stop("An IC50 estimate should always be positive.")
  }
  f <- -theta[3]*((theta[4]-theta[1])/((1+(x/theta[2])^theta[3])^2))*((x/theta[2])^(theta[3]-1))
  if(scale) {
    f <- scales::rescale(x = f, to = c(theta[4],theta[1]))
  }
  return(f)
}
#Second derivative derivation
d2_dr4pl <- function(x, theta, scale = F){
  if (any(is.na(theta))) {
    stop("One of the parameter values is NA.")
  }
  if (theta[2] <= 0) {
    stop("An IC50 estimate should always be positive.")
  }
  f <- 2*((theta[3]*(x/theta[2])^(theta[3]-1))^2)*((theta[4]-theta[1])/((1+(x/theta[2])^(theta[3]))^3))-theta[3]*(theta[3]-1)*((x/theta[2])^(theta[3]-2))*((theta[4]-theta[1])/((1+(x/theta[2])^theta[3])^2))
  if(scale) {
    f <- scales::rescale(x = f, to = c(theta[4],theta[1]))
    f <- f - f[1]
  }
  return(f)
}

ggplot(df, aes(x = dat_as_num)) +
  geom_hline(yintercept = 0) +
  [![enter image description here][1]][1]geom_point(aes(y = cases), color = "grey", alpha = .6, size = 5) +
  stat_function(fun = d1_dr4pl, args = list(theta = dr4pl_obj$parameters, scale = T), color = "red") +
  stat_function(fun = d2_dr4pl, args = list(theta = dr4pl_obj$parameters, scale = T), color = "blue") +
  stat_function(fun = dr4pl::MeanResponse, args = list(theta = dr4pl_obj$parameters), color = "gold") +
  geom_vline(xintercept = dr4pl_obj$parameters[2], linetype = "dotted") +
  theme_classic()

dr4pl fitted curve in yellow,

Как вы можете видеть, точка перегиба, которая является значением IC50 (тета 2) модели logisti c с 4 параметрами, хорошо выстраивается, когда мы приближаемся к этому way.

summary(dr4pl_obj)
#$call
#dr4pl.formula(formula = cases ~ dat_as_num, data = df, init.parm = c(30000, 18300, 2, 0))
#
#$coefficients
#               Estimate       StdErr       2.5 %      97.5 %
#Upper limit 25750.61451 4.301008e-05 25750.59681 25750.63221
#Log10(IC50) 18298.75347 4.301008e-09 18298.67889 18298.82806
#Slope        5154.35449 4.301008e-05  5154.33678  5154.37219
#Lower limit    58.48732 4.301008e-05    58.46962    58.50503
#
#attr(,"class")
#[1] "summary.dr4pl"

Кроме того, используя dr4pl, он говорит, что значение IC50 примерно 18298.8, что позднее 2020-02-06. Не слишком далеко от значения inflection. Я уверен, что может быть лучше использовать модель, чем 4pl, но я знал, что это была единственная модель, которую я мог бы написать первое и второе производные для ответа на этот вопрос.

Я уверен другие языки кодирования более специализированы, когда дело доходит до производных, и могут даже рассчитать их для вас, если вы начнете с начальной функции. Я думаю, что один из этих языков - mathematica.

A disclaimer , я закончил масштабирование первого и второго производных так, чтобы они могли быть построены вместе. Их фактические значения намного больше, чем показано здесь.

...