nlstools - nlsBoot 不工作

nlstools - nlsBoot not working

我正在尝试通过将 nlsBoot 函数应用于我创建的函数来使用 nlstools 包中的函数。从我的函数输出中使用 nlsBoot 时出现错误。但是,如果我直接将 nls 函数与它工作的数据一起使用。发生了什么事,有解决办法吗?

# Nonlinear function to generate data
NEE <- function(GPmax, alpha, resp, PAR) {
  ((alpha * PAR * GPmax)/((alpha * PAR)+ GPmax)) - resp

}

#some data
plot <- rep(c(1,2), each = 2000)
PAR <- 1:2000
dat <- data.table(plot, PAR)
dat[, GPP := (NEE(12, 0.73, -2, PAR) + rnorm(length(PAR), sd=2))]

library(nlstools)

# Function I created
model.fun <- function(df){
  fit <- nls(GPP ~ ((alpha * PAR * GPmax)/((alpha * PAR)+ GPmax)) - resp, 
             start = list(GPmax = 12, alpha = 0.73, resp = -2), data = df)

  return(list(b = summary(nlsBoot(fit))))
}

models <- dat[, list(model.fun(.SD)) , by = .(plot)]

# Error
# Error in data2[, var1] <- fitted1 + sample(scale(resid1, scale = FALSE),  : 
# object of type 'closure' is not subsettable

# Using nls directly outside of the function I created for plot 1.
mod1 <- nls(GPP ~ ((alpha * PAR * GPmax)/((alpha * PAR)+ GPmax)) - resp, 
           start = list(GPmax = 12, alpha = 0.73, resp = -2), data = dat[plot==1])

# Bootstrap of residuals
summary(nlsBoot(mod1,niter=5))

mod2 <- nls(GPP ~ ((alpha * PAR * GPmax)/((alpha * PAR)+ GPmax)) - resp, 
               start = list(GPmax = 12, alpha = 0.73, resp = -2), data = dat[plot==2])

    # Bootstrap of residuals
    summary(nlsBoot(mod2,niter=5))

# works

我会自己编写引导程序。这给了你更多的控制权,并且避免了这种范围界定问题:

library(boot)

model.fun <- function(dt){
  fit <- nls(GPP ~ ((alpha * PAR * GPmax)/((alpha * PAR)+ GPmax)) - resp, 
             start = list(GPmax = 12, alpha = 0.73, resp = -2), data = dt)
  #this copy should be avoided if you have big data, but I don't have time right now:
  df <- copy(dt)
  df[, fitted := fitted(fit)]
  df[, resid := residuals(fit)]
  fun <- function(df, inds) {
    df[, bootGPP := fitted + resid[inds]]
    tryCatch(coef(nls(bootGPP ~ ((alpha * PAR * GPmax)/((alpha * PAR)+ GPmax)) - resp, 
                      start = list(GPmax = 12, alpha = 0.73, resp = -2), data = df)),
             error = function(e) c("GPmax" = NA, "alpha" = NA, "resp" = NA))


  }
  b <- boot(df, fun, R = 1000)
  res0 <- b$t0
  res1 <- apply(b$t, 2, sd, na.rm = TRUE)
  res2 <- res0 - colMeans(b$t, na.rm = TRUE)
  return(as.list(setNames(c(sum(!is.na(b$t[,1])), res0, res1, res2), c("n", "GPmax", "alpha", "resp","SEGPmax", 
                                                 "SEalpha", "SEresp","BiasGPmax", "Biasalpha", "Biasresp"))))
}

set.seed(42)
models <- dat[, model.fun(.SD) , by = .(plot)]
#   plot    n    GPmax     alpha       resp  SEGPmax   SEalpha   SEresp   BiasGPmax   Biasalpha    Biasresp
#1:    1 1000 12.60382 0.9308744 -1.3579906 1.249449 0.2729928 1.263640 -0.11642169 -0.04376221 -0.11526536
#2:    2 1000 13.58702 0.8660450 -0.5081954 1.109599 0.2085150 1.125234 -0.06664517 -0.02241303 -0.06447617