1

我正在尝试使用 deSolve 拟合 ODE 模型来提取参数值 alpha。我有来自配对实验的两个物种的数据，我想获得一个物种在第二个物种上的相互作用强度。我正在使用强制函数 approxfun() 来插入第二个物种的时间序列并将其传递到模型中。我的问题是我不断得到的参数值与我最初的猜测相同。

t <- c(0,4,5,7,10,12,14,17,19,21,24,26,28,31,33,35,38,40)
N1 = c(4.000000,772.727273,4139.175258,3585.635359,22.909507, 0.000000,
       148.936170,45.454545,4.95049,35.353535,131.979695,46.082949,
       61.452514,87.804878,155.778895,80.000000,76.081788,120.418848)
N2 = c(0.00000,0.00000,0.00000,20.61856,97.36541,109.43912,195.74468,
       227.27273,69.30693,116.16162,157.36041,153.60983,195.53073,82.92683,
       50.25126,20.00000,19.02045,109.94764)

N2data <- as.data.frame(list(t, N2))
InterpN2 <- approxfun(N2data, method = "constant", rule = 2)

parms <- c(alpha = 0.9)

state = c(N1 = BTdataT$final.dens[1])
r = 1.5
K = 7000


logistic <- function(t, state, parameters) {
  with(as.list(c(state,parameters)),{

    dNdt = r*N1*(1-(N1/K))-(alpha*N1*InterpN2(t)) 

    list(c(dNdt), InterpN2(t)) 
   }) 
 }

Out <- ode(times = t, func = logistic, y = state, parms = parms)
plot(Out)

 y.obs  <- BTdataT$final.dens

 shootfun <- function(alpha){
  Out <- ode(times = times, func = logistic, y = state, parms = parms)

  y.pred <- Out[,2] 

  sigma  <- sqrt(sum((y.obs-y.pred)^2)/length(y.obs))

  nloglike <- -sum(dnorm(y.obs, mean = y.pred, sd = sigma, log = TRUE))
  return(nloglike)
}

shoot <- mle2(shootfun,start=list(alpha = parms[1]),
              method="Nelder-Mead")
summary(shoot)

我的预期结果将是从该数据中提取的实际 alpha 值，而不是收敛到初始条件的值。我对此很陌生，目前还不太了解。非常感谢你的帮助！！

4

0 回答 0