0
$\begingroup$

I measured the change in three traits (y1, y2, y3) over time (x). The three traits each reached a maximum at different times and then declined. I am interested in the phase until the max and would like to fit Gompertz curves in the following form: $$ y = \alpha \exp(\beta(1−\exp(−\gamma x))) $$

as the estimated parameters have biological meaning and will be used for further calculations. One approach I found in an (old) paper is a two-phase (until and after the max) fitting of Gompertz curves.

Major problems:

I could not determine until which time point of the data I would use for the fitting (i.e. to discard the decline phase):

  1. Within the same group, y1, y2 and y3 reached max at different times (graph below: y2 reached max earlier than y1 and y3).
  2. The decline is greater in some groups than others (graph below: the decline phase is more obvious in treatment a than in treatment b).

The graph shows the Gompertz curves I fit to all data (y1) and data until x = 55 (y2 and y3) for treatments a and b.

enter image description here

Attempt:

A solution I think of is to first fit broken-line regression to find the max point, and then fit the Gompertz curve to the data until that point. However, it is difficult to fit Gompertz curves to less than four time points of observations (e.g. for y3 in treatment b, the estimated upper asymptote of the Gompertz is higher biased).

Is there any better approach to fit the Gompertz curves? Any suggestions or comments would be greatly appreciated. Thanks a lot in advance.


Edit:

Sample data are pasted.

There are 7 time points each with ~10 observations, which also differed in time in a microscopic scale. In the dataset below, 3-4 observations per time are provided. After reading the comments , I am abandoning the idea of fitting Gompertz curves to those data (I did not expect the declines at different times).

I would like to compare 1) time to max among treatments and among traits, 2) max y among treatments 3) growth rate and 4) decline rate.

For 2), as they reached max y at different times, would it be a fair comparison if I just take the max y regardless of x? For 3) and 4), the same problem with time remains. As the $\Delta x$ are different, would it be more suitable to compare the instantaneous initial growth rate (at x = 0) and final decline rate (at x = 100) rathre than compare the averages from the initial to max point (or from max to final point)?

        x    y1    y2    y3
1   3.088 7.922 0.368 2.575
2   2.994 8.061 0.368 2.636
3   2.860 7.962 0.361 2.560
4   2.741 7.915 0.353 2.499
5  21.430 8.383 0.396 2.955
6  21.285 8.127 0.391 2.811
7  21.561 8.262 0.386 2.863
8  35.046 8.690 0.385 3.056
9  35.142 8.513 0.376 2.984
10 35.256 8.386 0.377 2.998
11 53.666 8.743 0.386 3.122
12 53.767 8.664 0.359 2.997
13 53.545 8.542 0.365 2.966
14 68.315 8.490 0.347 2.744
15 68.083 8.820 0.359 2.983
16 68.215 8.205 0.347 2.756
17 86.383 8.567 0.346 2.823
18 86.275 8.709 0.346 2.860
19 86.505 8.707 0.331 2.823
20 99.480 8.465 0.328 2.705
21 99.399 8.545 0.331 2.747
22 99.553 8.523 0.325 2.657
$\endgroup$
5
  • $\begingroup$ With two back-to-back Gompertz curves and an unspecified join-point, there's at least 7 parameters ... fitted to 7 times. Least squares curves will simply fit the sample mean at each point. $\endgroup$
    – Glen_b
    Nov 13, 2022 at 22:35
  • $\begingroup$ I would like to fit only one Gompertz curve from start to the max. The estimates are needed for other calculation. But yes there are few time points for 3 parameters. $\endgroup$
    – vetna
    Nov 17, 2022 at 5:52
  • $\begingroup$ This doesn't help in any way, in fact it may often be worse (consider if the peak precedes the third observation). I didn't wish to compound the problem by making the number of observations available be uncertain as well. $\endgroup$
    – Glen_b
    Nov 17, 2022 at 6:59
  • $\begingroup$ "One approach I found in an (old) paper is a two-phase (until and after the max) fitting of Gompertz curves." Could you give a reference to that paper. The approach seems a bit arbitrary and the solution to your problem might be to use a more modern approach. In any case, your data has only very little variation and growth is very difficult to express. You may use something like modeling an initial growth rate and the peak height, but because there are so few time points during which the curve increases it is difficult to fit any complicated growth mechanism. $\endgroup$ Nov 17, 2022 at 9:10
  • $\begingroup$ Thanks for the comments. I have added some detail in the post. I did not expect the early decline at different times. Here I will first compare the max and the growth/decline rates. It seems segmented regression with 2 linear lines would be an easy solution, but least squares curve should give a better fit (?) The old paper I referred to is this one (also where the reparameterization comes from; two stage-fitting shown in figure 3; although I am not using Gompertz here, I would also like to know a better approach): spo.nmfs.noaa.gov/sites/default/files/pdf-content/1976/743/… $\endgroup$
    – vetna
    Nov 17, 2022 at 16:20

3 Answers 3

1
$\begingroup$

After the comments on my first answer and the data added later to the original question I am afraid that the whole has to be reconsidered. That is why a different answer is posted in order to avoid confusion.

They are two key points :

  • The Gompertz model cannot be used to model the whole data.

  • The Gompertz model might be convenient for only a part of data at low values of $x$. But they are not enough values of $x$ sufficiently different one to another. One cannot fit the Gompertz equation for signifiant computed values of the parameters. I think that any variant of method to select a set of points at low $x$ is a dead end.

Of course one cannot take for granted that the Gompertz model is convenient in all cases even if they are a lot of publications. $$y = \alpha \exp(\beta(1−\exp(−\gamma x))) \tag 1$$ They are other competitive models (Makeham, Gavrilov, ...).

For example (not exclusively) I suggest to try the model $(2)$ of the Makeham kind. Note that those two models become equivalent at low values of $x$. $$y = \alpha \exp\Big(-\delta x+\beta\big(1−\exp(−\gamma x)\big)\Big) \tag 2$$ which is equvalent to : $$\ln(y)=A+B\:x+C\:e^{\:p\:x} \quad \begin{cases} A=\ln(\alpha)+\beta\\ B=-\delta\\ C=-\beta\\ p=-\gamma \end{cases} \tag 3$$ This function can be fitted to the whole data thanks to non-linear regression that is an iterative method starting from "guessed" values of parameters.

A not iterative method not requiring "guessed" parameters is shown below.

enter image description here

This is an application of the theory from https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales

RESULTS :

enter image description here

enter image description here

enter image description here

NOTE :

As a matter of fact the method of "Regression with Integral Equation" is a linear regression but not wrt the model equation itself (Eq.3) but wrt a linear integral equation (below) to which the model equation is solution. In the present case :

enter image description here

The discret values of the integral S(x) are first computed by numerical calculus. Then p,c2,c1,c0 are obtained by linear regression. Only the numerical value of p is used latter to compute the discret values of exp(px). Then Eq.(3) becomes linear wrt A,B,C which are finally obtained by linear regression.

$\endgroup$
4
  • $\begingroup$ The 'not iterative method' should only be used the get starting values right? It doesn't find a least squares solution. $\endgroup$ Nov 18, 2022 at 11:36
  • $\begingroup$ @Sextus Empiricus. This depends on the wanted criteria of fitting. The above "not iterative" method is a least squares but not directly wrt the function to be fitted but to a linear integral equation which solution is the fitted function. Often no criteria of fitting is clearly specified when a problem is raised. This is a pity because this opens the floodgates to arbitrariness. On the other hand if no criteria of fitting is specified any one can be chosen to answer. Then why not using the "not iterative" method ? Of course an iterative method is required to satisfy a specific criteria. $\endgroup$
    – JJacquelin
    Nov 18, 2022 at 13:47
  • $\begingroup$ Thank you very much for your helpful and detailed solution. I cannot express my gratitude enough! The Gompertz-Makeham approach solves several of my problems at once. It fitted the data nicely and I am now looking for references which I could cite in my field. The "not iterative method" also saved a lot of time. It is a quick way to generate a fit. Thank you! $\endgroup$
    – vetna
    Mar 27, 2023 at 14:44
  • $\begingroup$ Don't mention it. I added a note at the end of my answer to explain the application of the method of Regression with Integral Equation in the present case. For more information about the theory and applications refer to the document quoted in my answer. $\endgroup$
    – JJacquelin
    Mar 28, 2023 at 10:33
1
$\begingroup$

Obviously your data has to be modeled with a function which is increasing at low value of $x$ and decreasing at high values of $x$.

On the other hand the function $$y = \alpha \exp(\beta(1−\exp(−\gamma x)))\tag 1$$ is always increasing. Thus the fitting of such a function is wrong. That way the estimated parameters will have no physical meaning.

It is strange that the experimental measurements doesn't agree even roughly with the proposed theoretical law. You should clarify this point. But this is not the subject to be discussed here.

If you maintain that the function $(1)$ is a correct model only for low values of $x$ and that your measurements are false at high values of $x$ then I understand why you want to eliminate the points above a maximum.

I think that isn't a good way to obtain significant result. Nevertheless you could first choose another model without physical meaning but using a function which has a maximum (for example the quadratic function). Fitting this "toy" model will give a maximum $(x_m\:,\:y_m)$. Then you could fit the function $(1)$ with only the points $x\leq x_m$.

I would have liked to test this proposed method with your data. But this was not possible with only your graphs and without numerical data.

$\endgroup$
3
  • $\begingroup$ It seems you have posted the same answer twice. $\endgroup$ Nov 16, 2022 at 10:14
  • $\begingroup$ I cannot understand why double post occurs. I delate one. $\endgroup$
    – JJacquelin
    Nov 16, 2022 at 10:19
  • $\begingroup$ Thanks very much for your reply! I have posted some data (unsure if they are sufficient?) above. Function (1) is from literature and is indeed fitted to early phase (low x). I have tried segmented regression, which fitted two lines and gave me the x where the slope changes. I will try the quadratic function to see if the results are similar. My main concern is after removing the data, it would be difficults to fit Gompertz as the upper asymptote is gone... $\endgroup$
    – vetna
    Nov 17, 2022 at 6:09
0
$\begingroup$

Here are two approaches:

1) loess Since the edit indicates that we do not need any specific parametric form we can use a nonparametric loess fit. (Alternately a nonparametric gam fit, not shown, as in this answer could be tried.) Using the data in the Note at the end and using R code, we define a prediction function, pred, that accepts a loess fit lo and the x values and returns the corresponding y values. We then run it for each of the y columns plotting the result and the maximum. The rates are the average growth and decline rates before and after the maximum is reached.

span <- 0.6
pred <- function(lo, xx) predict(lo, data.frame(x = xx))
rngx <- range(DF$x)
xx <- seq(rngx[1], rngx[2], length = 200)
lo <- list()
result <- matrix(NA, 3, 4, 
  dimnames = list(names(DF)[-1], c("xAtMax", "ymax", "growth", "decline")))

opar <- par(mfrow = c(3, 1), mar = c(2, 4, 1, 1))
for(yi in names(DF)[-1]) {
  plot(DF$x, DF[[yi]], xlab = "", ylab = yi)
  fo <- reformulate("x", yi)  # generate formula needed by loess
  lo[[yi]] <- loess(fo, DF, span = span)
  lines(pred(lo[[yi]], xx) ~ xx, DF, col = "red")
  opt <- optimize(pred, rngx, lo = lo[[yi]], maximum = TRUE)
  with(opt, abline(h = objective, v = maximum, lty = 2))
  result[yi, "xAtMax"] <- opt$maximum
  result[yi, "ymax"] <- opt$objective
  result[yi, "growth"] <- (opt$objective - pred(lo[[yi]], rngx[1])) / 
    (opt$maximum - rngx[1])
  result[yi, "decline"] <- (opt$objective - pred(lo[[yi]], rngx[2])) / 
    (opt$maximum - rngx[2])
}
par(opar)
result

giving:

     xAtMax      ymax      growth      decline
y1 51.26490 8.6533897 0.014244016 -0.002362361
y2 21.94926 0.3911388 0.001521732 -0.000803318
y3 46.50930 3.0538271 0.011206357 -0.006279221

(continued after chart) screenshot

2) Two regression lines We can use two regression lines and take the minimum of them as our curve fit. This follows my answer here. It performs linear regression fits on the first 10 points and separately on the last 10 points and then uses those fits as starting values for a nonlinear regression that takes the minimum of two lines fitting all the data at once. In this case we can solve for the intersection of the lines and read the growth and decline right off the slopes of the lines.

result <- matrix(NA, 3, 4, 
  dimnames = list(names(DF)[-1], c("xAtMax", "ymax", "growth", "decline")))
opar <- par(mfrow = c(3, 1), mar = c(2, 4, 1, 1))
fm <- list()
for(yi in names(DF)[-1]) {
  y <- DF[[yi]]
  fm_1 <- lm(y ~ x, DF, subset = 1:10)
  fm_2 <- lm(y ~ x, DF, subset = seq(to = nrow(DF), length = 10))
  st <- list(a = coef(fm_1)[[1]], b = coef(fm_1)[[2]],
       c = coef(fm_2)[[1]], d = coef(fm_2)[[2]])
  fm[[yi]] <- nls(y ~ pmin(a + b * x, c + d * x), DF, start = st)
  co <- coef(fm[[yi]])
  X <- with(as.list(co), (a - c) / (d - b))
  Y <- with(as.list(co), a + b * X)
  plot(DF$x, DF[[yi]], xlab = "", ylab = yi)
  abline(a = co[[1]], b = co[[2]], col = "red")
  abline(a = co[[3]], b = co[[4]], col = "red")
  abline(h = Y, v = X, lty = 2)
  result[yi, "xAtMax"] <- X
  result[yi, "ymax"] <- Y
  result[yi, "growth"] <- co[["b"]]
  result[yi, "decline"] <- co[["d"]]
}
par(opar)
result

giving

     xAtMax      ymax      growth       decline
y1 41.81144 8.6327074 0.017346330 -0.0014532939
y2 21.89605 0.3917683 0.001544078 -0.0008119577
y3 32.38736 3.0602008 0.016732095 -0.0049408889

screenshot

Note

Lines <- "x    y1    y2    y3
1   3.088 7.922 0.368 2.575
2   2.994 8.061 0.368 2.636
3   2.860 7.962 0.361 2.560
4   2.741 7.915 0.353 2.499
5  21.430 8.383 0.396 2.955
6  21.285 8.127 0.391 2.811
7  21.561 8.262 0.386 2.863
8  35.046 8.690 0.385 3.056
9  35.142 8.513 0.376 2.984
10 35.256 8.386 0.377 2.998
11 53.666 8.743 0.386 3.122
12 53.767 8.664 0.359 2.997
13 53.545 8.542 0.365 2.966
14 68.315 8.490 0.347 2.744
15 68.083 8.820 0.359 2.983
16 68.215 8.205 0.347 2.756
17 86.383 8.567 0.346 2.823
18 86.275 8.709 0.346 2.860
19 86.505 8.707 0.331 2.823
20 99.480 8.465 0.328 2.705
21 99.399 8.545 0.331 2.747
22 99.553 8.523 0.325 2.657"
DF <- read.table(text = Lines)
$\endgroup$
1
  • $\begingroup$ Thank you for the clear examples and code. This was originally what I asked for. $\endgroup$
    – vetna
    Mar 27, 2023 at 14:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.