基于MLE拟合年龄结构化水痘SIR模型的数值参数错误排查

我尝试用最大似然估计（MLE）将年龄结构化的确定性SIR模型拟合到水痘发病数据中。模型单独运行时表现正常且结果合理，但拟合执行过程中出现错误。

R代码 #

var.data <- read.csv("/Users/laurenadams/Documents/gitrepo/Varicella/Input CSV/Age Specific Incidence.csv")

source("/Users/laurenadams/Documents/gitrepo/Varicella/Version 1/Equations/No Vax.R")
source("/Users/laurenadams/Documents/gitrepo/Varicella/Version 1/Functions.R")
source("/Users/laurenadams/Documents/gitrepo/Varicella/Version 1/Force of Infection.R")
source("/Users/laurenadams/Documents/gitrepo/Varicella/Version 1/Parameters.R")

#set initial conditions
M_inits     <- P[1]
Sv0_inits   <- P[2] - 0.000049*P[2] - 0.5*P[2] - 0.000012*P[2] - 0.20*P[2] - 0.20*P[2]      
Sv_inits    <- P_new[-1]-0.000049*P_new[2:100] - 0.5*P[2:100] - 0.000012*P_new[2:100] - 0.20*P_new[2:100] - 0.20*P_new[2:100]
Iv0_inits   <- 0.000049*P[2]
Iv_inits    <- 0.000049*P_new[2:100] 
Sz0_inits   <- 0.5*P[2]-0.000012*P[2]   
Sz_inits    <- 0.5*P[2:100]-0.000012*P_new[2:100]
Iz0_inits   <- 0.000012*P[2]   
Iz_inits    <- 0.000012*P_new[2:100] 
Rz0_inits   <- 0.40*P[2]
Rz_inits    <- 0.40*P_new[2:100] 
B10_inits   <- 0    
B1_inits    <- 0*P_new[2:100]
VInc0_inits <- 0
VInc_inits  <- c(rep(0,99))
ZInc0_inits <- 0
ZInc_inits  <- c(rep(0,99))

inits_all <- c(M=M_inits, Sv0=Sv0_inits, Sv=Sv_inits,
           Iv0=Iv0_inits,Iv=Iv_inits,
           Sz0 = Sz0_inits,Sz=Sz_inits,
           Iz0=Iz0_inits,Iz=Iz_inits,
           Rz0=Rz0_inits,Rz=Rz_inits,
           B10=B10_inits, B1=B1_inits,
           VInc0=VInc0_inits, VInc=VInc_inits,
           ZInc0=ZInc0_inits, ZInc=ZInc_inits)


# mle.sir - Maximum likelihood estimation function for closed SIR model
mle.sir <- function(b) {
  t <- seq(0, 365)    
  beta <- exp(b)

  results <- as.data.frame(ode(y= inits_all, 
                           times=t, 
                           varicella_fun_novax,
                           parms=c(parms)), row.names=F)
  Y <- results[365,]
  Y <- Y[103:202]
 nll <- -sum(dpois(x=var.data$Case_Proportion, lambda=Y, log=TRUE))
 return(nll)
 #return(results)
}

# initial - Initial estimates of beta and gamma
initial <- list(b=beta)
# fit0 - preliminary fit of model to data using the initial estimates
fit0 <- mle2(mle.sir, start=initial)

报错信息 #

Error in dpois(x = var.data$Case_Proportion, lambda = Y, log = TRUE) : 
  Non-numeric argument to mathematical function

解决思路 #

• 修正Y的类型转换：results[365,]提取的是数据框的一行，属于数据框类型，而dpois要求lambda为数值向量。修改提取方式，直接将目标列转为数值：

  # 替换原Y的赋值代码
  Y <- as.numeric(results[365, 103:202])
  

• 验证观测数据类型：检查var.data$Case_Proportion是否为数值型，执行str(var.data$Case_Proportion)查看类型。如果是因子或字符型，转换为数值：

  var.data$Case_Proportion <- as.numeric(as.character(var.data$Case_Proportion))
  

• 检查维度匹配：确保Y和var.data$Case_Proportion的长度一致，可在函数中加入cat(length(Y), length(var.data$Case_Proportion), "\n")验证。

• 排查模型输出合理性：确认Y的所有取值为非负数（泊松分布要求lambda > 0），若存在负数或NA，需检查模型方程、初始条件或参数设置是否正确。

• 调试技巧：临时在mle.sir函数中加入str(Y)和str(var.data$Case_Proportion)，直接调用mle.sir(initial$b)运行函数，查看两者的类型和数值情况，定位具体问题。

内容的提问来源于stack exchange，提问作者Lauren Adams