我正在尝试将幼鹅的生长数据拟合到 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")