在NLME中拟合数据的技巧？ [英] Tricks for fitting data in nlme?
问题描述
当我在nlme中放入数据时，我第一次尝试都没有成功，在 nlme（fit.model）之后，我习惯于看到诸如此类的东西：
When I fit data in nlme, I never succeed on the first try, and after nlme(fit.model) I am accustomed to seeing things such as:
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的单位轴（例如，从几年到几天，或者从几天到成长度的天数）。
1)Change the units of the x-axis (e.g. from years to days, or days to growing degree days).
2）在我的数据集中进行ax = 0，y = 0的测量
2)Make a x=0, y=0 measurement in my dataset
3）添加 random = pdDiag（）
4）弄乱什么是随机的，什么是固定的
4)Mess with what is random and what is fixed
5）切入我的数据集并尝试在不同时间拟合不同的部分
5)Chop up my dataset and try to fit different parts at different times
6）实现一个非常简单的拟合，然后使用 update 使模型正确
6)Achieve a very simple fit, then use update to make the model proper
好像 工作。还有其他人要添加到此列表吗？什么可以帮助您使nlme处理数据？
Eventually something seems to work. Does anyone else have something to add to this list? What helps you get nlme to work with your data?
我知道这个问题可能已经解决，但是如果有任何建议如何改写它，以便于接受所以，我将不胜感激。
I realize this question will probably be closed, but if there are any suggestions on how to reword it to be acceptable to SO, I would appreciate the input.
这里是一个示例，其中我尝试了其中一些方法，但到目前为止还没有成功：
Here is an example where I have tried some of these things, but have not had success so far:
数据： 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）：
And here is the function (SSbgf) for the curve:
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"))
推荐答案
另一个窍门是增加pnls值
Another trick is to increase the pnls tolerance.
所需的代码为：
control = nlmeControl(pnlsTol = x, msVerbose = TRUE)
pnls公差的起始值为0.001，所以我喜欢以0.01或0.02开始。只需将x替换为您的数字就可以了。
The starting value for pnls tolerance is 0.001, so I like to start with 0.01 or 0.02. Just replace x with your number and you should be set.
这篇关于在NLME中拟合数据的技巧？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
