通过使用 optim() 最小化残差平方和来优化参数，将模型拟合到观察数据

问题描述

我试图将这个非常简单的 4 种线性 Lotka-Volterra 竞争模型拟合到观察到的数据中，但由于某种原因，当我尝试 optim() 函数时，关于 deSolve 的某些东西似乎失败了。

# Data
data <- data.frame(Cod = c(0.1966126,0.1989563,0.2567677,0.3158896,0.4225435,0.7219856,1.0570824,0.7266830,0.6286763,0.6389475),Herring = c(1.988372,2.788014,3.397138,2.557245,2.627013,3.045617,3.161002,3.531306,3.432021,3.617174),Sprat = c(2.030273,3.480469,3.009277,1.895996,2.457520,1.991211,2.350098,2.118164,1.693359,1.869141),Flounder = c(0.4758220,0.4425532,0.4185687,0.4967118,0.7102515,0.5733075,0.7404255,0.5996132,0.6235977,0.7187621))
# Model formulation
LLV <- function(time,state,parameters) {
  with(as.list(c(state,parameters)),{
    db1.dt = b1*(r1+a11*b1+a12*b2+a13*b3+a14*b4)
    db2.dt = b2*(r2+a22*b2+a21*b1+a23*b3+a24*b4)
    db3.dt = b3*(r3+a33*b3+a31*b1+a32*b2+a34*b4)
    db4.dt = b4*(r4+a44*b4+a41*b1+a42*b2+a43*b3)
    list(c(db1.dt,db2.dt,db3.dt,db4.dt))
  })
}
# Model input and simulation
# Model input
params <- c(r1 = -0.342085,r2 = 0.6855681,r3 = 2.757769,r4 = 0.9744113,a11 = -1.05973762,a12 = 0.09577309,a13 = -0.01915480,a14 = 1.36098939,a21 = 0.17533326,a22 = -0.32247342,a23 = 0.03111628,a24 = 0.30212711,a31 = 0.5303516,a32 = -0.4869761,a33 = -0.3194882,a34 = -1.5089027,a41 = 0.004418133,a42 = 0.163716414,a43 = -0.237873378,a44 = -1.519158802)
ini <- c(b1 = data[1,1],b2 = data[1,2],b3 = data[1,3],b4 = data[1,4])
tmax <- 10
t <- seq(1,tmax,0.1)
# Results and first parameter guess is more or less okay
results <- deSolve::ode(y = ini,times = t,func = LLV,parms = params)
matplot(data,pch = 1)
matplot(x = results[,y = results[,-1],type = "l",add = TRUE)

在这里，我继续编写一个函数，该函数最小化残差平方和，当包含在 optim() 中时，使用上述初始参数猜测应该会产生我正在寻找的结果。

min.RSS <- function(data,params) {
  output <- deSolve::ode(y = ini,parms = params)
  predictions <- exp(output[,-1])
  observations <- data
  return(sum((predictions-observations)^2))
}
result <- optim(par = params,fn = min.RSS,data = data)
fit <- deSolve::ode(y = ini,parms = result$par)
matplot(x = fit[,y = fit[,lwd = 3,add = TRUE)

任何有关如何解决此问题的想法将不胜感激。

解决方法

你更适合，但你应该非常小心这个问题。我有点疯狂并使用（开发中）fitode package 来解决这个问题。我拟合了模型并得到了更好的拟合，还尝试在我的最佳拟合周围使用 100 个随机变化的起点进行拟合。你的残差平方和是 1.19； fitode 在第一次尝试时达到 0.29，100 次拟合中最好的结果是 RSS=0.16。然而：这些拟合高度不稳定。该图显示了对未来 5 个时间步长的数据和预测的拟合（1）您的拟合（虚线）；（2）fitode初始拟合（虚线）； (3) 其他100个fitode拟合（最佳拟合0.05RSS以内的为实心，比这差的画得很轻）。

您可以看到样本外预测大多是疯狂的。你的拟合实际上比一些更好的拟合更稳定 - 在整个社区崩溃之前到达时间步骤 13 - 但底线是在这种情况下非常适合数据绝不保证一个明智的答案。看起来 100 个拟合中的一个在没有崩溃的情况下到达了预测时间序列的末尾（这似乎是基于观察到的时间序列的合理合理的“常识”预测）。

为了可靠地拟合这些数据，您要么需要一个参数少得多的模型，要么需要以先验形式提供的外部信息，要么正则化 - 某种方式进行惩罚拟合，这意味着 '摆动的确定性轨迹，或不合理的交互参数/增长率。p>

## remotes::install_github("parksw3/fitode")
library(fitode)

## data with tags for fitode
data2 <- setNames(data,paste0(names(data),"_obs"))
data2 <- data.frame(times=seq(nrow(data2)),data2)


## Model formulation (for fitode)
LV_model <- odemodel(
    name="4-species LV",model=list(
        Cod ~ Cod*(r1+a11*Cod+a12*Herring+a13*Sprat+a14*Flounder),Herring ~ Herring*(r2+a22*Herring+a21*Cod+a23*Sprat+a24*Flounder),Sprat ~ Sprat*(r3+a33*Sprat+a31*Cod+a32*Herring+a34*Flounder),Flounder ~ Flounder*(r4+a44*Flounder+a41*Cod+a42*Herring+a43*Sprat)
    ),observation=list(
        Cod_obs ~ ols(mean=Cod),Herring_obs ~ ols(mean=Herring),Sprat_obs ~ ols(mean=Sprat),Flounder_obs ~ ols(mean=Flounder)
    ),initial=list(
        Cod ~ data2$Cod_obs[1],Herring ~ data2$Herring_obs[1],Sprat ~ data2$Sprat_obs[1],Flounder ~ data2$Flounder_obs[1]
    ),link=setNames(rep("identity",length(pars)),pars),par= pars
)

## plot results
plotres <- function(p,ODEint="rk",lty=1,dt=0.1,tvec=seq(1,10,by=dt),...) {
    par(las=1,bty="l")
    res <- deSolve::ode(ini,tvec,LLV,p,method=ODEint)
    matplot(res[,1],res[,-1],type="l",lty=lty,...)
    return(invisible(res[,-1]))
}

f1 <- fitode(
    LV_model,data=data2,start=params,control=list(maxit=1e5,trace=1000)
)

## fitode with multistart

ranfit <- function(n,fit,range=0.5) {
    ## 
    rpars <- params*runif(length(params),1-range,1+range)
    newfit <- try(update(fit,start=rpars))
    return(newfit)
}

cl <- makeCluster(10)
clusterSetRNGStream(cl = cl,101)
clusterExport(cl,c("params","LV_model","data2"))
clusterEvalQ(cl,invisible(library(fitode)))
system.time(
    multifit <- parLapply(cl,1:100,ranfit,fit=f1,tvec=tvec)
)
saveRDS(multifit,file="SO65440448_multifit.rds")

ivec <- seq_along(multifit)
ivec <- ivec[sapply(multifit,function(x) !inherits(x,"try-error"))]
coef <- pred <- vector("list",length=length(ivec))
ll <- conv <- rep(NA,length(ivec))
for (i in seq_along(ivec)) {
    nf <- multifit[[ivec[i]]]
    coef[[i]] <- coef(nf)
    pp <- predict(nf,times=1:10)
    pred[[i]] <- cbind(times=pp[[1]][,do.call(cbind,lapply(pp,"[",-1)))
    ll[i] <- logLik(nf)
    conv[i] <- nf@mle2@details$convergence
}

par(las=1,bty="l")
matplot(pred[[1]][,pred[[1]][,type="n",ylim=c(0,6),xlab="time",ylab="density")
lthresh <- 0.05
for (i in 1:length(pred)) {
    good <- ll[i]>(max(ll)-lthresh)
    alpha <- if (good) 0.8 else 0.1
    lwd <- if (good) 2 else 1
    matlines(pred[[i]][,pred[[i]][,col=adjustcolor(palette()[1:4],alpha.f=alpha),lwd=lwd)
}
matpoints(data2[,data2[,pch=16,cex=3)
plotres(optimres$par,add=TRUE,lwd=3,lty=2,dt=1)
plotres(coef(f1),lty=3,dt=1)

对于那些感兴趣的人，我设法获得了一个涉及更改 ode 集成方法的解决方案。这是工作优化器：

# Optimising parameter fit
LVmse = function(parms) {
  out = as.matrix(deSolve::ode(ini,1:10,parms,method="rk")[,-1])
  RSS = sum((spp-out)^2,na.rm = TRUE) # Minimising residual sum of squares
  return(RSS)
}
optimres <- optim(par = params,fn = LVmse)

model-fitting ode optimization r r