我用于求解和绘制 ODE 耦合系统的 R 代码有什么问题？

Question

I'm reasonably new to R and have this system of ODEs
\frac{dS}{dT} = -\beta(1-\mu)S(t)\frac{I(t))}{1-D(t))}
\frac{dI}{dT} = -\beta(1-\mu)S(t)\frac{I(t))}{1-D(t))}-\delta I(t))-\phi I(t))
\frac{dR}{dT}=\delta I(t))
\frac{dD}{dT}=\phi I(t))

我正在尝试创建 S、I、R 和 D 相对于时间的图，并利用我的 R 基础知识以及 deSolve 手册中类似示例的帮助编写了此代码。然而，我的代码生成的图并不是我所期望的。

（注意。我使用的参数值是任意的，我计划使用平方和方法将这些参数调整为我的数据集中的真实参数值，即来自墨西哥的 COVID-19 病例和死亡数据，而且我也命名错误我的代码中的参数字母，它们代表相同的值）

This is my code: 


```
install.packages("deSolve")
library("deSolve")

sird_equations <- function(time, variables, parameters) {
with(as.list(c(variables, parameters)), {
dS <- -beta * (1- mu) * S * (I/(1-D))
dI <- -beta * (1- mu) * S * (I/(1-D))-sigma * I - theta * I
dR <- sigma * I 
dD <- theta * I 
return(list(c(dS, dI, dR, dD)))
})
} 

parameter_values <- c( 
beta = 0.95, #infectious contact rate (/person/day)
sigma = 0.5, #recovery rate (/day)
theta = 0.043, #mortality rate (/day)
mu = 0.1 #adjustment for efficacy of countermeasures
)

inital_values <- c(
S = 126000000, #suseptible population (adjust based on ylur data)
I = 1, #infected on day 1 of data 
R = 0, #recovered on day 1
D = 0 #dead on day 1
)


time_values <- seq(0, 10) #days(adjust based on your data)
ls()

sird_values_1 <- ode(
y = inital_values,
times = time_values,
func = sird_equations,
parms = parameter_values 
)

sird_values_1 <- as.data.frame(sird_values_1)
with(sird_values_1, {
plot(time, S, type = "l", col = "blue",
xlab = "time (days)", ylab = "# of people")
lines(time, I, col = "black")
lines(time, R, col = "red")
lines(time, D, col = "orange")
})

legend("right", c("susceptible", "infectious", "recovered", "dead"),
col = c("blue", "black", "red", "orange"), lty = 1, bty = "n")
```

我期待这样的剧情 但是我得到的是这个

我尝试过调整参数值，但我所做的任何事情都对情节没有太大影响，我非常感谢支持理解我哪里出了问题。

Answer 1

这是一个基于

wikipedia

上描述的 SIRD 模型的模型，并进行

mu

调整：

sird_equations <- function(time, variables, parameters) {
  with(as.list(c(variables, parameters)), {
    dS <- - (beta * (1 - mu) * S * I ) / (S + I + R + D)
    dI <- (beta * (1 - mu) * S * I) / (S + I + R + D) - sigma * I - theta * I
    dR <- sigma * I 
    dD <- theta * I 
    return(list(c(dS, dI, dR, dD)))
  })
} 

我改变了时间步长：

time_values <- seq(0, 100, 0.01)

然后我还更改了绘图代码，以采用正确的 y 限制：

with(sird_values_1, {
  plot(
    time, S, 
    ylim = c(0, max(sird_values_1)),
    type = "l", col = "blue",  xlab = "time (days)", ylab = "# of people"
  )
  lines(time, I, col = "black")
  lines(time, R, col = "red")
  lines(time, D, col = "orange")
})

或使用

ggplot2

:

sird_values_1 |> 
  tidyr::pivot_longer(-time) |> 
  ggplot(aes(time, value, color = name)) +
  geom_line() +
  scale_color_manual(values = c(S = "blue", I = "black", R = "red", D = "orange")) +
  scale_y_continuous(labels = scales::label_number()) +
  labs(y = '# of people', x = 'time (days)', color = NULL) +
  theme_classic()