向量化 mle2 模型的应用

Question

我已经编写了一个模型，我通过 mle2 包使用 ML 拟合数据。但是，我有一个大型样本数据框，我想将模型拟合到每个复制品，然后在数据框中检索模型的所有系数。

我尝试使用 plyr 包中的 ddply 函数没有成功。

我在尝试时收到以下错误消息：

Error in output[[var]][rng] <- df[[var]] : 
  incompatible types (from S4 to logical) in subassignment type fix

有什么想法吗？

这是我正在做的一个例子。

这是我的数据框。我在Pond 5...n 上测量了day 1....n。测量由 143 个通量 (flux.cor) 组成，这是我正在建模的变量。

     Pond Obs                Date     Time   Temp       DO   pH U day month    PAR
932    5 932 2011-06-16 17:31:00 17:31:00 294.05 334.3750 8.47 2   1     1 685.08
933    5 933 2011-06-16 17:41:00 17:41:00 294.05 339.0625 8.47 2   1     1 808.44
934    5 934 2011-06-16 17:51:00 17:51:00 294.02 340.6250 8.46 2   1     1 752.78
935    5 935 2011-06-16 18:01:00 18:01:00 294.00 340.6250 8.45 2   1     1 684.14
936    5 936 2011-06-16 18:11:00 18:11:00 293.94 340.9375 8.50 2   1     1 625.86
937    5 937 2011-06-16 18:21:00 18:21:00 293.88 341.5625 8.48 2   1     1 597.06
    day.night Treat            H  pOH           OH   DO.cor   sd.DO    av.DO   DO.sat
932         1     A 3.388442e-09 5.53 2.951209e-06 342.1406 2.63078 342.1406 274.0811
933         1     A 3.388442e-09 5.53 2.951209e-06 339.0625 2.63078 342.1406 274.0811
934         1     A 3.467369e-09 5.54 2.884032e-06 340.6250 2.63078 342.1406 274.2432
935         1     A 3.548134e-09 5.55 2.818383e-06 340.6250 2.63078 342.1406 274.3513
936         1     A 3.162278e-09 5.50 3.162278e-06 340.9375 2.63078 342.1406 274.6763
937         1     A 3.311311e-09 5.52 3.019952e-06 341.5625 2.63078 342.1406 275.0020
      DO_flux      NEP.hr  flux.cor  sd.flux    av.flux
932 -3.078125 -3.09222602 -3.078125 2.104482 -0.1070312
933  1.562500  1.54903673  1.562500 2.104482 -0.1070312
934  0.000000 -0.01375489  0.000000 2.104482 -0.1070312
935  0.312500  0.29876654  0.312500 2.104482 -0.1070312
936  0.625000  0.61126617  0.625000 2.104482 -0.1070312

这是我的模型：

    # function that generates predictions of O2 flux given GPP R and gas exchange
flux.pred <- function(GPP24, PAR, R24, Temp, U, DO, DOsat){
    # calculates Schmidt coefficient from water temperature
    Sc<-function(Temp){
        S<-0.0476*(Temp)^2 + 3.7818*(Temp)^2 - 120.1*Temp + 1800.6
        }
    # calculates piston velocity k (m h-1) from wind speed at 10m (m s-1)
    k600<-function(U){
        k.600<-(2.07 + 0.215*((U)^1.7))/100 
        }
    # calculates piston velocity k (m h-1) from wind speed at 10m (m s-1)
    k<-function(Temp,U){
        k<-k600(U)*((Sc(Temp)/600)^-0.5)
        }
    # physical gas flux (mg O2 m-2 10mins-1)
    D<-function(Temp,U,DO,DOsat){
        d<-(k(Temp,U)/6)*(DO-DOsat)
    }   

  # main function to generate predictions   
flux<-(GPP24/sum(YSI$PAR[YSI$PAR>40]))*(ifelse(YSI$PAR>40, YSI$PAR, 0))-(R24/144)+D(YSI$Temp,YSI$U,YSI$DO,YSI$DO.sat)
return(flux)
}

返回通量的预测。

然后我构建我的似然函数：

   # likelihood function
ll<-function(GPP24, PAR, R24, Temp, U, DO.cor, DO.sat){
    pred = (flux.pred(GPP24, PAR, R24, Temp, U, DO.cor, DOsat))
    pred = pred[-144]
    obs = YSI$flux.cor[-144]
    return(-sum(dnorm(obs, mean=pred, sd=sqrt(var(obs-pred)))))
    } 

并应用它

ll.fit

一天对一个池塘效果很好，但我想做的是自动将它应用于所有池塘。

我尝试了 ddply（如上所述）

metabolism<-ddply(YSI, .(Pond,Treat,day,month), summarise,
mle = mle2(ll,start=list(GPP24=100, R24=100)))

但没有成功。我也尝试使用 for 循环提取系数，但这也不起作用。

for(i in 1:length(unique(YSI$day))){
GPP<-numeric(length=length(unique(YSI$day)))
GPP[i]<-mle2(ll,start=list(GPP24=100, R24=100))
    }

我们将不胜感激地接受任何帮助。

Answer 1

您的函数至少存在一个问题：在您的函数flux.pred 或您的函数中没有一个参数可以实际指定所使用的数据。你硬编码了它。那么到底是怎么回事 *ply 应该猜测它需要将 YSI$... 更改为子集？

除此之外，正如@hadley 指出的那样，ddply 不适合您。 dlply 可能，或者您可能只使用 by() 或 lapply(split()) 的经典方法。

想象一下你做了一个函数

flux.pred <- function(data, GPP24, R24){
    # calculates Schmidt coefficient from water temperature
    Sc<-function(data$Temp){
        S<-0.0476*(data$Temp)^2 ...
    ...
    }   

还有一个函数

ll<-function(GPP24, R24, data ){
    pred = (flux.pred(data, GPP24, R24 ))
    pred = pred[-144] # check this
    obs = data$flux.cor[-144] # check this
    return(-sum(dnorm(obs, mean=pred, sd=sqrt(var(obs-pred)))))
    } 

然后你应该能够做例如：

dlply(data, .(Pond,Treat,day,month), .fun=function(i){
    mle2(ll,start=list(GPP24=100, R24=100, data=i))
})

数据参数的传递取决于您在 mle2 中用于优化的内容。在您的情况下，您使用默认优化器，即optim。有关详细信息，请参阅?optim。参数data=i 将从mle2 传递到optim 到ll。

我无法检查的是 optim 的行为方式。甚至可能是您的功能并没有真正按您的意愿工作。通常你应该有一个类似的函数：

ll <- function(par, data){
    GPP24 <- par[1]
    R24 <- par[2]
    ...
}

为了优化工作。但如果你说它像你写的那样有效，我相信你。确保它确实如此。我不相信...

附带说明：既不使用 by() / lapply(split()) 也不使用 dlply() 与矢量化相同。相反，所有这些结构都是内在循环。关于使用它们的原因，请阅读：Is R's apply family more than syntactic sugar?