Как написать функцию, где коэффициент стремится к 0 в R?

Question

Я хотел бы написать функцию, которая сглаживает коэффициент скорости роста до 0 за 60 дней. До сих пор мне удалось написать следующий код:

corona <- data.frame(Cases = c(3, 16, 79, 157, 229, 322, 400, 650, 888, 1128, 1694, 2036, 2502, 3089, 3858), Date = seq(as.Date("2020/02/20"), as.Date("2020/03/05"), by = "days"))

library(dplyr)
corona_entire <- corona %>% mutate(Growth = (Cases - lag(Cases))/lag(Cases)*100)  

mean(corona_entire$Growth[12:15])

ff = function(x) x*(1.2285823)^60

ff(3858)

Однако в моей функции скорость роста (0,2285823) постоянна в течение 60 периодов. Я хотел бы сказать R, чтобы эта скорость роста была равна 0, когда мы приближаемся к 60. Мне нужно написать в основном функцию сходимости для скорости роста.

Есть идеи, как я могу ее кодировать?

Спасибо!

Answer 1

В дополнение к моему комментарию выше, мне не ясно, что вы пытаетесь сделать. Если вы хотите смоделировать ставку Growth, вам нужно соответствовать какой-либо форме модели.

Для начала, как насчет экспоненциальной модели вида y = y0 * exp(k * time)?

В этом случае мы можем линеаризовать модель (и данные), взяв журнал, а затем использовать lm для оценки коэффициентов модели log(y0) и k.

df <- corona_entire %>% mutate(Time = as.integer(Date - min(Date)))
fit <- lm(log(Growth) ~ Time, weights = df$Growth, data = df)

Здесь я использовал взвешенные наименьшие квадраты подходить, взвешивая каждую точку по ее скорости Growth.

Затем мы можем построить точки и кривую наилучшего соответствия:

f <- function(x, fit) exp(coef(fit)[1])*exp(coef(fit)[2] * x)
ggplot(df, aes(Time, Growth)) +
    geom_point() +
    stat_function(fun = f, args = list(fit = fit)) +
    labs(x = sprintf("Days since %s", min(df$Date)))

Не совсем подходит, но это должно дать вам некоторые идеи. Вы, вероятно, хотите подобрать более подходящую нелинейную модель скорости роста и оценить параметры, используя nls.

---

Обновление

Поскольку вы действительно хотите предсказать Cases, давайте переформулируем нашу модель.

Мы начнем снова с экспоненциальной модели вида Cases ~ y0 * exp(k * Time)

ggplot(df, aes(Time, Cases)) +
    geom_point()
fit1 <- lm(log(Cases) ~ Time, data = df)
f1 <- function(x, fit) exp(coef(fit)[1])*exp(coef(fit)[2] * x)
ggplot(df, aes(Time, Cases)) +
    geom_point() +
    stat_function(fun = f1, args = list(fit = fit1)) +
    labs(x = sprintf("Days since %s", min(df$Date)))

Не совсем подходит! Результаты указывают на субэкспоненциальный рост. Простая модель субэкспоненциального роста в эпидемиологии - это модель вида Cases ~ (r / m * Time + A)^m, см., Например, Chowell et al., Phys. Life Rev. 18, 66 (2016) .

Итак, давайте подгоним модель, на этот раз с помощью нелинейной процедуры наименьших квадратов nls.

fit2 <- nls(
    Cases ~ (r / m * Time + A)^m,
    data = df,
    start = list(r = 4, m = 3, A = 1))
f2 <- function(x, r, m, A) (r / m * x + A)^m
ggplot(df, aes(Time, Cases)) +
    geom_point() +
    stat_function(
        fun = f2, 
        args = list(
            r = coef(fit2)[1],
            m = coef(fit2)[2],
            A = coef(fit2)[3])) +
    labs(x = sprintf("Days since %s", min(df$Date)))

Выглядит как приличная посадка. Вы можете проверить качество подгонки и нелинейных оценок наименьших квадратов для коэффициентов, если вы наберете summary(fit2)

summary(fit2)
#
#Formula: Cases ~ (r/m * Time + A)^m
#
#Parameters:
#  Estimate Std. Error t value Pr(>|t|)
#r   2.3308     0.6543   3.562  0.00391 **
#m   3.3316     0.4202   7.929 4.12e-06 ***
#A   2.1101     0.3126   6.750 2.04e-05 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 51.41 on 12 degrees of freedom
#
#Number of iterations to convergence: 6
#Achieved convergence tolerance: 6.514e-07
#

Answer 2

Глядя на данные, выглядит как квадратичная c кривая - лучший вариант для модели Cases как функции days

corona$days = as.numeric(corona$Date - corona$Date[1], "days") + 1
mod = lm(Cases ~ poly(days, 2, raw = TRUE), corona)
summary(mod)

#Call:
#lm(formula = Cases ~ poly(days, 2, raw = TRUE), data = corona)

#Residuals:
#    Min      1Q  Median      3Q     Max 
#-140.48  -50.63  -24.30   65.89  148.04 

#Coefficients:
#                           Estimate Std. Error t value Pr(>|t|)    
#(Intercept)                 264.912     84.071   3.151  0.00836 ** 
#poly(days, 2, raw = TRUE)1 -158.269     24.179  -6.546 2.75e-05 ***
#poly(days, 2, raw = TRUE)2   25.863      1.469  17.600 6.17e-10 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Residual standard error: 94.38 on 12 degrees of freedom
#Multiple R-squared:  0.9949,   Adjusted R-squared:  0.9941 
#F-statistic:  1181 on 2 and 12 DF,  p-value: 1.668e-14


plot(corona$days, corona$Cases)
lines(predict(mod, data.frame(days = corona$days)))

# Growth Rate 
d = predict(mod, data.frame(days = 59:60))
diff(d)/d[1]
#         2 
#0.03606188 

Answer 3

Если вы просто хотите, чтобы линейное падение темпов роста к 1 за 60 дней, вы можете сделать это:

ff = function(initial_n, initial_rate = 1.2285823, days = 60, time_to_stasis = 60)
{
  daily_rate <- seq(initial_rate, 1, length.out = time_to_stasis)
  result <- numeric(days)
  result[1] <- initial_n
  for(i in seq(days - 1)) result[i + 1] <- floor(daily_rate[i] * result[i])
  return(result)
}

Таким образом, вы получите ежедневное число, подобное этому:

ff(3858)
#>  [1]    3858    4739    5803    7084    8620   10456   12643   15239   18309   21926
#> [11]   26173   31141   36932   43656   51436   60403   70699   82477   95897  111129
#> [21]  128350  147743  169494  193790  220818  250760  283791  320073  359754  402961
#> [31]  449796  500332  554607  612621  674331  739644  808418  880454  955498 1033237
#> [41] 1113297 1195248 1278600 1362812 1447290 1531398 1614460 1695773 1774611 1850239
#> [51] 1921922 1988936 2050581 2106192 2155151 2196899 2230944 2256873 2274360 2283171

и вы можете настроить параметры так, как вам нравится.

Вы можете использовать его для построения таких проекций:

plot(1:60, ff(3858))

Я не уверен, насколько это биологически правдоподобно.