3

我正在尝试将幼鹅的生长数据拟合到 Gompertz 模型,但我想将它固定在 y 截距处,在这种情况下,它对应于 72.1 的体重t = 0。采用 inSSgompertz中使用的参数化形式nls,我将方程 $Asym * (e^{-b3 * b2^x}$ 转换为 $\frac{X_0}{exp(-b3)} * (e^{-b3 * b2 ^x}$ 通过设置 x = 0,等于 y 截距并求解Asym,因此我在模型中将其替换为另一边的项:$\frac{X_0}{exp(-b3)}$

问题是,对于这个改变的模型,我不能再使用从 SSgompertz 获得的起始值:

library(nlme)
fbmgrp   <- groupedData(bm   ~ age | dummy, data = bmdata) # I put the data at the bottom of my question

# Defining the Gompertz model with fixed y-intercept
gobm   <- function(Age, b2, b3) ((72.1 / exp(-b2)) * (exp(-b2 * (b3^Age))))
gobm   <- deriv(body(gobm)[[2]], namevec = c("Age", "b2", "b3"), function.arg = gobm)

fm1 <- nls(bm ~ SSgompertz(age, Asym, b2, b3), data = bmgrp) 
# Using parameter estimates for b2 and b3 from SSgompertz in the customized Gompertz model leads to an error:
fm2 <- nls(bm ~ gobm(age, b2, b3), data = bmgrp, start = list(b2 = coef(fm1)[2], b3 = coef(fm1)[3]), na.action = na.exclude)

qr.default(.swts * attr(rhs, "gradient")) 中的错误:外部函数调用中的 NA/NaN/Inf (arg 1)
此外:警告消息:在 log(b3) 中:NaN 产生
5:qr。 default(.swts * > attr(rhs, "gradient")) 4: qr(.swts * attr(rhs, "gradient"))
3: assign("QR", qr(.swts * attr(rhs, "gradient) ")), envir = thisEnv)
2: (function (newPars) {setPars(newPars) assign("resid", .swts * (lhs - assign("rhs", getRHS(), envir = thisEnv))), envir = thisEnv) assign("dev", sum(resid^2), envir = thisEnv) assign("QR", qr(.swts * attr(rhs, "gradient")), envir = thisEnv) return(QR$rank < min(dim(QR$qr)))})(c(7.52734494474091, -0.123046281607752))
1: nls(bm ~ gobm(age, I, K), data = bmgrpd,开始 = 列表(I = coef(fm13a)[2],K = coef(fm13a)[3]),na.action = na.exclude)

是否有一种方法可以从 SSgompertz 的参数估计中获得具有固定 y 截距的 Gompertz 模型的合适初始值? ?SSgompertz说拐点b2与 x = 0 处的函数值有关。在 Ricklefs bookchapter 'Avian postnatal development' (1983) 中,他展示了 Gompertz 模型的不同参数化:$M(x) = M(∞) * e^{-[log M(∞) - log M(0)] * e^(-K * x)}$ 所以渐近线 M(∞) orAsym和 y 截距 M(0) or的差值X0表示拐点 b2:$b2 = log(Asym) - log(X0) <--> X0 = Asym - e^{b2}$ 我用 fm1 的系数检查了这个关系:

coef(fm1) Asym b2 b3 1079.4598299 8.6815201 0.8710308 但coef(fm1)[1] - exp(coef(fm1)[2])导致-4813.538其在 0 附近的真实 y 截距一点也不接近:SSgompertz(0, Asym = coef(fm1) 1 , b2 = coef(fm1)[2], b3 = coef(fm1)[3])

0.1831767

这是与 SSgompertz 的参数的拟合图: Gompertz 模型与体重测量值的拟合

因此,如果有人可以帮助我找到具有固定 y 截距的 Gompertz 模型的初始值,我将不胜感激!

以下是数据:

bmdata <- read.table(header = TRUE, stringsAsFactors = TRUE, text = "
pop year    sex age bm  dummy
RUS 2013    f   7   137 1
NL  2012    f   10  132 1
RUS 2013    f   10  225 1
RUS 2013    f   12  183 1
RUS 2013    f   19  362 1
RUS 2003    f   19  880 1
RUS 2005    f   19  750 1
RUS 2005    f   19  625 1
RUS 2004    f   20  670 1
RUS 2004    f   20  900 1
RUS 2004    f   20  740 1
NL  2012    f   21  260 1
RUS 2015    f   21  510 1
RUS 2014    f   22  660 1
RUS 2015    f   22  500 1
RUS 2003    f   22  776 1
RUS 2003    f   22  890 1
RUS 2003    f   22  1010    1
RUS 2004    f   22  612 1
RUS 2004    f   22  316 1
RUS 2004    f   22  622 1
RUS 2004    f   22  438 1
RUS 2005    f   22  684 1
RUS 2005    f   22  652 1
RUS 2014    f   23  695 1
RUS 2015    f   23  520 1
RUS 2015    f   23  450 1
RUS 2003    f   23  854 1
RUS 2003    f   23  897 1
RUS 2003    f   23  1180    1
RUS 2004    f   23  990 1
RUS 2004    f   23  1110    1
RUS 2004    f   23  1080    1
RUS 2004    f   23  500 1
RUS 2004    f   23  480 1
RUS 2005    f   23  820 1
NL  2012    f   24  190 1
RUS 2013    f   24  400 1
RUS 2014    f   24  740 1
RUS 2014    f   24  775 1
RUS 2014    f   24  840 1
RUS 2015    f   24  490 1
RUS 2004    f   24  870 1
RUS 2004    f   24  516 1
RUS 2004    f   24  672 1
RUS 2004    f   24  604 1
RUS 2005    f   24  820 1
RUS 2005    f   24  910 1
NL  2015    f   25  490 1
RUS 2013    f   25  620 1
RUS 2014    f   25  880 1
RUS 2014    f   25  860 1
RUS 2015    f   25  930 1
RUS 2003    f   25  980 1
RUS 2004    f   25  680 1
RUS 2004    f   25  1070    1
RUS 2004    f   25  940 1
RUS 2004    f   25  900 1
RUS 2004    f   25  910 1
RUS 2014    f   26  985 1
RUS 2014    f   26  882 1
RUS 2014    f   26  860 1
RUS 2014    f   26  945 1
NL  2004    f   26  480 1
RUS 2003    f   26  900 1
RUS 2004    f   26  720 1
RUS 2004    f   26  634 1
RUS 2014    f   27  800 1
RUS 2014    f   27  830 1
RUS 2014    f   27  940 1
RUS 2014    f   27  950 1
RUS 2014    f   27  980 1
RUS 2015    f   27  600 1
RUS 2003    f   27  1080    1
RUS 2003    f   27  762 1
RUS 2003    f   27  880 1
RUS 2003    f   27  1020    1
RUS 2003    f   27  950 1
RUS 2003    f   27  1050    1
RUS 2003    f   27  1140    1
RUS 2003    f   27  1080    1
RUS 2003    f   27  1060    1
RUS 2004    f   27  1300    1
RUS 2004    f   27  630 1
RUS 2005    f   27  990 1
NL  2012    f   28  237 1
RUS 2014    f   28  790 1
RUS 2014    f   28  1005    1
RUS 2014    f   28  855 1
RUS 2014    f   28  1005    1
RUS 2015    f   28  660 1
RUS 2003    f   28  966 1
RUS 2003    f   28  902 1
RUS 2003    f   28  740 1
RUS 2003    f   28  960 1
RUS 2003    f   28  1140    1
RUS 2004    f   28  744 1
RUS 2005    f   28  517 1
NL  2015    f   29  430 1
NL  2015    f   29  425 1
RUS 2014    f   29  1100    1
RUS 2014    f   29  1115    1
RUS 2014    f   29  1000    1
RUS 2014    f   29  1055    1
RUS 2014    f   29  1100    1
RUS 2014    f   29  1020    1
RUS 2014    f   29  1125    1
RUS 2014    f   29  1245    1
RUS 2014    f   29  985 1
RUS 2015    f   29  980 1
RUS 2004    f   29  612 1
RUS 2004    f   29  860 1
RUS 2005    f   29  510 1
RUS 2005    f   29  902 1
RUS 2005    f   29  910 1
RUS 2014    f   30  1000    1
RUS 2014    f   30  1000    1
RUS 2014    f   30  1100    1
RUS 2014    f   30  1200    1
RUS 2014    f   30  1030    1
RUS 2014    f   30  980 1
RUS 2014    f   30  1000    1
RUS 2014    f   30  890 1
RUS 2004    f   30  965 1
RUS 2004    f   30  720 1
RUS 2004    f   30  1090    1
RUS 2004    f   30  840 1
RUS 2005    f   30  1146    1
NL  2012    f   31  430 1
RUS 2014    f   31  1140    1
RUS 2014    f   31  1225    1
RUS 2015    f   31  1000    1
RUS 2015    f   31  1020    1
RUS 2015    f   31  1130    1
RUS 2015    f   31  730 1
RUS 2004    f   31  1024    1
RUS 2004    f   31  704 1
RUS 2005    f   31  910 1
RUS 2005    f   31  870 1
NL  2012    f   32  418 1
NL  2012    f   32  282 1
RUS 2014    f   32  1140    1
RUS 2014    f   32  1125    1
RUS 2014    f   32  1085    1
RUS 2015    f   32  1200    1
NL  2004    f   32  530 1
RUS 2003    f   32  1140    1
RUS 2004    f   32  640 1
RUS 2005    f   32  830 1
RUS 2005    f   32  1050    1
RUS 2005    f   32  1140    1
RUS 2005    f   32  1070    1
RUS 2005    f   32  1050    1
RUS 2005    f   32  930 1
RUS 2014    f   33  1090    1
RUS 2014    f   33  1310    1
RUS 2014    f   33  1105    1
RUS 2014    f   33  1060    1
RUS 2015    f   33  820 1
RUS 2015    f   33  750 1
RUS 2003    f   33  850 1
RUS 2004    f   33  476 1
RUS 2004    f   33  810 1
RUS 2004    f   33  1050    1
RUS 2004    f   33  935 1
NL  2015    f   34  560 1
NL  2012    f   34  438 1
NL  2012    f   34  440 1
RUS 2014    f   34  990 1
RUS 2014    f   34  1220    1
RUS 2014    f   34  1090    1
RUS 2014    f   34  1085    1
RUS 2015    f   34  820 1
RUS 2004    f   34  1042    1
RUS 2004    f   34  980 1
RUS 2004    f   34  930 1
RUS 2013    f   35  860 1
RUS 2014    f   35  1210    1
RUS 2014    f   35  1040    1
RUS 2014    f   35  1080    1
RUS 2014    f   35  1270    1
RUS 2004    f   35  1020    1
RUS 2005    f   35  1150    1
NL  2015    f   36  480 1
NL  2015    f   36  530 1
NL  2012    f   36  460 1
RUS 2015    f   36  720 1
RUS 2015    f   36  1060    1
RUS 2005    f   36  1000    1
RUS 2005    f   36  1170    1
RUS 2013    f   37  1030    1
RUS 2013    f   37  1050    1
RUS 2014    f   37  1095    1
RUS 2014    f   37  1025    1
RUS 2014    f   37  1250    1
RUS 2015    f   37  1060    1
RUS 2004    f   37  1120    1
RUS 2014    f   38  1140    1
RUS 2014    f   38  935 1
RUS 2014    f   38  1430    1
RUS 2014    f   38  1225    1
RUS 2014    f   38  1225    1
RUS 2014    f   38  1315    1
NL  2012    f   39  535 1
NL  2012    f   39  195 1
NL  2012    f   39  1090    1
RUS 2014    f   39  1325    1
RUS 2015    f   39  1110    1
RUS 2015    f   39  1030    1
RUS 2015    f   39  960 1
RUS 2015    f   39  980 1
RUS 2014    f   40  1230    1
RUS 2014    f   40  1300    1
RUS 2015    f   40  1040    1
RUS 2003    f   40  1180    1
NL  2012    f   41  500 1
RUS 2015    f   41  1190    1
NL  2004    f   41  1000    1
NL  2012    f   42  580 1
NL  2015    f   44  770 1
NL  2012    f   44  580 1
NL  2012    f   44  700 1
NL  2005    f   44  1050    1
NL  2012    f   45  680 1
NL  2004    f   47  1210    1
NL  2004    f   47  1250    1
NL  2012    f   48  480 1
NL  2012    f   49  605 1
NL  2012    f   49  407 1
NL  2012    f   50  1060    1
NL  2004    f   50  1320    1
NL  2004    f   50  940 1
NL  2012    f   51  1325    1
NL  2012    f   52  1220    1
NL  2012    f   53  880 1
NL  2012    f   53  880 1
NL  2012    f   53  1300    1
NL  2012    f   54  695 1
NL  2012    f   54  640 1
NL  2012    f   54  830 1
NL  2004    f   54  910 1
NL  2004    f   54  1110    1
NL  2012    f   56  540 1
NL  2012    f   56  1225    1
NL  2012    f   56  1280    1
NL  2015    f   57  1080    1
NL  2012    f   57  1180    1
NL  2015    f   58  980 1
NL  2015    f   58  1040    1
NL  2012    f   58  740 1
NL  2015    f   59  1040    1
NL  2015    f   59  1030    1
NL  2015    f   59  1140    1
NL  2012    f   59  1360    1
NL  2012    f   61  900 1")
4

0 回答 0