r - 在 nlme 中拟合数据的技巧？

Question

当我在 nlme 中拟合数据时，我第一次尝试从未成功，之后 nlme(fit.model)我习惯于看到这样的事情:

Error in nlme.formula(model = mass ~ SSbgf(day, w.max, t.e, t.m), random = list( :
  step halving factor reduced below minimum in PNLS step

Error in MEestimate(nlmeSt, grpShrunk) : 
  Singularity in backsolve at level 0, block 1

所以我回去

1) 更改 x 轴的单位(例如，从年到天，或从天到生长度日)。

2)在我的数据集中进行 x=0, y=0 测量

3)添加random=pdDiag()
4)混淆随机的和固定的

5)切碎我的数据集并尝试在不同的时间拟合不同的部分

6)实现非常简单的拟合，然后使用update使模型正确

最终似乎有些东西奏效了。还有其他人要添加到此列表中吗？什么可以帮助您让 nlme 处理您的数据？

我意识到这个问题可能会被关闭，但如果有任何关于如何改写它以被 SO 接受的建议，我将不胜感激。

下面是一个例子，我尝试了其中的一些东西，但到目前为止还没有成功:

资料:https://www.dropbox.com/s/4inldx7617fip01/proots.csv .这已经只是整个系列的一部分。

编码:

roots<-read.table("proots.csv", header = TRUE)

#roots$day[roots$year == 2007] <- 0 #when I use a dataset with time=0, mass=0
roots$day[roots$year == 2008] <- 153
roots$day[roots$year == 2009] <- 518
roots$day[roots$year == 2010] <- 883
roots$day[roots$year == 2011] <- 1248
roots$day[roots$year == 2012] <- 1613
roots$day[roots$year == 2013] <- 1978

#or bigger time steps
roots$time[roots$year == 2008] <- 1
roots$time[roots$year == 2009] <- 2
roots$time[roots$year == 2010] <- 3
roots$time[roots$year == 2011] <- 4
roots$time[roots$year == 2012] <- 5
roots$time[roots$year == 2013] <- 6

roots$EU<- with(roots, factor(plot):factor(depth)) #EU is "experimental unit"
rootsG<-groupedData(mass ~ day | EU, data=roots)

#I will post the SSbgf function below -- run it first
fit.beta <- nlsList(mass ~ SSbgf(day, w.max, t.e, t.m), data = rootsG) 

fit.nlme.bgf<-nlme(fit.beta)
fit.nlme.bgf<-nlme(fit.beta, random=list(w.max + t.e + t.m ~1))
fit.nlme.bgf<-nlme(fit.beta, random=list(w.max ~ 1))
fit.nlme.bgf<-nlme(fit.beta, random=pdDiag(w.max ~1))

fit.nlme.bgf<-nlme(fit.beta, random=pdDiag(w.max + t.e + t.m ~1))
fit.nlme.bgf<-nlme(fit.beta, random=list(t.m ~1)) 
fit.nlme.bgf<-nlme(fit.beta, random=list(t.e ~1))

fit.nlme.bgf<-nlme(fit.beta, random=pdSymm(w.max ~1))
fit.nlme.bgf<-nlme(fit.beta, random=pdDiag(w.max ~1))

这是曲线的函数 (SSbgf):

bgfInit <- function(mCall, LHS, data){

  xy <- sortedXyData(mCall[["time"]], LHS, data)
  if(nrow(xy) < 4){
    stop("Too few distinct input values to fit a bgf")
  }
  w.max <- max(xy[,"y"])
  t.e <- NLSstClosestX(xy, w.max)
  t.m <- NLSstClosestX(xy, w.max/2)
  value <- c(w.max, t.e, t.m)
  names(value) <- mCall[c("w.max","t.e","t.m")]
  value

}


bgf <- function(time, w.max, t.e, t.m){

  .expr1 <- t.e / (t.e - t.m)
  .expr2 <- (time/t.e)^.expr1
  .expr3 <- (1 + (t.e - time)/(t.e - t.m))
  .value <- w.max * .expr3 * .expr2

  ## Derivative with respect to t.e
  .exp1 <- ((time/t.e)^(t.e/(t.e - t.m))) * ((t.e-time)/(t.e-t.m) + 1)
  .exp2 <- (log(time/t.e)*((1/(t.e-t.m) - (t.e/(t.e-t.m)^2) - (1/(t.e - t.m)))))*w.max
  .exp3 <- (time/t.e)^(t.e/(t.e-t.m))
  .exp4 <- w.max * ((1/(t.e-t.m)) - ((t.e - time)/(t.e-t.m)^2))
  .exp5 <- .exp1 * .exp2 + .exp3 * .exp4 

  ## Derivative with respect to t.m
  .ex1 <- t.e * (time/t.e)^((t.e/(t.e - t.m))) * log(time/t.e) * ((t.e - time)/(t.e -     
 t.m) + 1) * w.max
  .ex2 <- (t.e - time) * w.max * (time/t.e)^(t.e/(t.e-t.m))
  .ex3 <- (t.e - t.m)^2
  .ex4 <- .ex1 / .ex3 + .ex2 / .ex3

  .actualArgs <- as.list(match.call()[c("w.max", "t.e", "t.m")])

##  Gradient
  if (all(unlist(lapply(.actualArgs, is.name)))) {
    .grad <- array(0, c(length(.value), 3L), list(NULL, c("w.max", 
                                                      "t.e", "t.m")))
    .grad[, "w.max"] <- .expr3 * .expr2
    .grad[, "t.e"] <- .exp5
    .grad[, "t.m"] <- .ex4 
    dimnames(.grad) <- list(NULL, .actualArgs)
    attr(.value, "gradient") <- .grad
  }
    .value
}

SSbgf <- selfStart(bgf, initial = bgfInit, c("w.max", "t.e", "t.m"))

Answer 1

另一个技巧是增加 pnls 容差。

所需的代码是:

control = nlmeControl(pnlsTol = x, msVerbose = TRUE)

pnls 容差的起始值是 0.001，所以我喜欢从 0.01 或 0.02 开始。只需用您的号码替换 x 就可以了。