在数据框中的每个组上运行模型并将关键指标输出到新的数据框中

Question

我有一个以下model，我目前正在为一种数据状态运行。来自这个博客的代码：https://blog.ephorie.de/epidemiology-how-contagious-is-novel-coronavirus-2019-ncov

它根据现有的每日 COVID-19 感染数据对未来的感染进行建模。我有一个美国所有县的数据集，并希望逐个县运行相同的分析，并将关键变量输出到每个县的数据集中。这是我为一个状态运行的代码。

library(deSolve)
library(tidyverse)

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)

cp <- data.frame(date,fips,Infected,day,N)

还有model：

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

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

Opt <- optim(c(0.5, 0.5), RSS, method = "L-BFGS-B", lower = c(0, 0), upper = c(1, 1))
Opt_par <- setNames(Opt$par, c("beta", "gamma"))
t <- 1:365 # time in days
fit <- data.frame(ode(y = init, times = t, func = SIR, parms = Opt_par))
R0 <- setNames(Opt_par["beta"] / Opt_par["gamma"], "R0")
infections <- fit[fit$I == max(fit$I), "I", drop = FALSE]
deaths <- max(fit$I) * 0.02

这是我拥有的县级数据的一小部分数据集：

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

我想为每个 fips（县代码）运行上述模型。作为输出，我想要一个dataframe，其中包括： time_stamp、fips、max(Infected)、max(day)、N、beta、gamma、R0、infections、deaths

我想要一个dplyr 管道解决方案并尝试使用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)))

错误：

ode(y = c(S, I, R), times = Day, func = SIR, parms = parameters) 中的错误： 找不到对象“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(.))