Подгонка модели к наблюдаемым данным путем оптимизации параметров с использованием optim(), минимизирующей остаточную сумму квадратов

У вас лучше подходит, но вы должны быть очень осторожны с этой проблемой. Я немного сошел с ума и использовал (находящийся в разработке) пакет fitode, чтобы решить эту проблему. Я подогнал модель и получил гораздо лучшее соответствие, а также попытался подобрать 100 случайно меняющихся начальных точек вокруг моего наилучшего соответствия. Ваша остаточная сумма квадратов составила 1,19; fitode получил значение 0,29 с первой попытки, а лучшее из 100 совпадений было RSS=0,16. Однако эти соответствия крайне нестабильны. На этом графике показаны соответствия данным и прогнозы на 5 временных шагов в будущем для (1) вашего соответствия (пунктирные линии); 2 – начальная подгонка фитода (пунктир); (3) 100 других подгонок фитодов (те, что в пределах 0,05 RSS от наилучшего соответствия, выделены сплошным цветом, те, что хуже этого, нарисованы очень слабо).

Вы можете видеть, что прогнозы вне выборки в основном сумасшедшие. Ваша подгонка на самом деле более стабильна, чем некоторые из лучших подгонок — она достигает шага 13 по времени до того, как все сообщество выйдет из строя — но суть в том, что в данном случае это хорошее соответствие данным никоим образом не гарантирует разумного ответа. Похоже, что одна из 100 подгонок достигает конца временного ряда прогноза без коллапса (что кажется разумным прогнозом здравого смысла, основанным на наблюдаемом временном ряду).

Для надежной подгонки этих данных вам либо нужна модель с гораздо меньшим количеством параметров, либо внешняя информация, предоставленная в виде априорных значений, либо регуляризация — какой-то способ наложить штраф на подгонки, подразумевающие «волнистую» детерминированность. траектории или параметры взаимодействия/темпы роста, которые являются необоснованными.

## 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]][,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]][,1],pred[[1]][,-1],
        type="n",lty=1,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]][,1],pred[[i]][,-1],lty=1,
             col=adjustcolor(palette()[1:4],alpha.f=alpha),
             lwd=lwd)
}
matpoints(data2[,1],data2[,-1],pch=16,cex=3)
plotres(optimres$par,add=TRUE, lwd=3,lty=2,dt=1)
plotres(coef(f1),add=TRUE, lwd=3,lty=3,dt=1)

Для тех, кто заинтересован, мне удалось получить решение, которое включает изменение метода интеграции оды. Вот рабочий оптимизатор:

# Optimising parameter fit
LVmse = function(parms) {
  out = as.matrix(deSolve::ode(ini, 1:10, LLV, 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)