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

Question

У меня есть следующая модель , которую я сейчас использую для одного состояния данных. Это код из этого блога: https://blog.ephorie.de/epidemiology-how-contagious-is-novel-coronavirus-2019-ncov

Он моделирует будущие инфекции на основе существующих ежедневных данных о заражении COVID-19. У меня есть набор данных для всех округов США, и я хочу запустить один и тот же округ анализа по округам и вывести ключевые переменные в набор данных для каждого округа. Вот код для одного штата data:

date <- '2020-03-24','2020-03-25','2020-03-26','2020-03-24','2020-03-25','2020-03-26'
fips <- 1001,1001,1001,1002,1002,1002
Infected <- 1,2,4,4,7,9
day <- 1,2,3,1,2,3
N <- 55601,55601,55601,2231,2231,2231

Я бы хотел запустить приведенную выше модель для каждого слияния (код округа). В качестве вывода я хочу иметь фрейм данных , который включает в себя: time_stamp, fips, max(Infected), max(day), N, beta, gamma, R0, infections, deaths

Я бы хотел использовать dplyr pipe и пытался использовать group_modify(), но постоянно получал ошибку. Можете ли вы помочь ??

Вот что у меня есть и ошибка, которую я получаю:

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))
    })
}

RSS <- function(parameters) {
    names(parameters) <- c("beta", "gamma")
    out <- ode(y = c(S, I, R)
               , times = Day, func = SIR, parms = parameters)
    fit <- out[ , 3]
    sum((Infected - fit)^2)
}

cp <- cp %>%
    group_by(fips) %>%
    mutate(S = N-Infected[1], I = Infected[1], R = 0) %>%
    group_modify(optim(c(0.5, 0.5), RSS, method = "L-BFGS-B", lower = c(0, 0), upper = c(1, 1)))

Ошибка:

Ошибка в оде (у = c (S , I, R), времена = день, веселье c = SIR, parms = параметры): объект 'S' не найден

Answer 1

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

library(deSolve)
library(tidyverse)

cp <- data.frame(
  date     = c('2020-03-24','2020-03-25','2020-03-26','2020-03-24','2020-03-25','2020-03-26'),
  fips     = c(1001,1001,1001,1002,1002,1002),
  Infected = c(1,2,4,4,7,9),
  day      = c(1,2,3,1,2,3),
  N        = c(55601,55601,55601,2231,2231,2231)
)

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

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

fit_model <- function(data) {
  init <- slice(data, 1) %>% select(S, I, R) %>% unlist()
  opt <- optim(c(0.5, 0.5), RSS,
               init = init,
               data = data,
               method = "L-BFGS-B",
               lower = c(0, 0),
               upper = c(1, 1))

  data.frame(alpha=opt$par[1], beta=opt$par[2], RSS=opt$value,
             conv= opt$convergence, ok = opt$convergence == 0)
}

cp %>%
  group_by(fips) %>%
  mutate(S = N[1]-Infected[1], I = Infected[1], R = 0) %>%
  group_modify( ~ fit_model(.))