I am trying to fit growth data from juvenile geese to the Gompertz model, but I would like to fix it at the y-intercept,
which would correspond in this case to body mass of 72.1 for t = 0. Taking the parametization form used in SSgompertz in nls,
I transformed the equation $Asym * (e^{-b3 * b2^x}$ into $\frac{X_0}{exp(-b3)} * (e^{-b3 * b2^x}$ by setting x = 0,
equalling with the y-intercept and solving for Asym, which I consequently replaced in the model by the term of the other side: $\frac{X_0}{exp(-b3)}$
The problem is that with this altered model I cannot use the starting values obtained from SSgompertz anymore:
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)
Error in qr.default(.swts * attr(rhs, "gradient")) : NA/NaN/Inf in foreign function call (arg 1)
In addition: Warning message: In log(b3): NaNs produced
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, start = list(I = coef(fm13a)[2], K = coef(fm13a)[3]), na.action = na.exclude)
Is there maybe a way to get from parameter estimates from SSgompertz to suitable initual values for a Gompertz model with fixed y-intercept?
?SSgompertz says that the inflection point b2 is related to the value of the function at x = 0.
In Ricklefs bookchapter ‘Avian postnatal development’ (1983), he displays a different parameterization of the Gompertz model:
$M(x) = M(∞) * e^{-[log M(∞) - log M(0)] * e^(-K * x)}$
So the difference between asymptote M(∞) or Asym and y-intercept M(0) or X0 represents the inflection point b2:
$b2 = log(Asym) - log(X0) <--> X0 = Asym - e^{b2}$
I checked this relationship out with the coefficients from fm1:
coef(fm1) Asym b2 b3 1079.4598299 8.6815201 0.8710308 But
coef(fm1)[1] - exp(coef(fm1)[2])resulted in-4813.538which is not in the least near the true y-intercept near 0: SSgompertz(0, Asym = coef(fm1)1, b2 = coef(fm1)[2], b3 = coef(fm1)[3])
0.1831767
Here a plot of the fit with the parameters from SSgompertz: Fit of Gompertz model to body mass measurements
So I would greatly appreciate if someone can help me in finding initial values for a Gompertz model with fixed y-intercept!
Here are the data:
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")
