Estimation of parameters when adding a weibull 'growth spurt' to a increasing trend


#1

Hello,

I want to see if i can separate out a growth spurt from a linear growth. In the
following code I add some gompertz growth to a sqrt growth curve

library(data.table);
library(ggplot2)
library(brms)
addGrowth <- function(x,w,a,b,l){
  f <-sapply(x-w[1], function(s) if(s<0) 90 else s)
  log(a) - gompertz(f, log(a),b,l)
}
set.seed(3032)
z <- data.table(x=1:150)
dur <- c(50, 90);
asymp <- 1.15;b <-30;l <- 0.30
noise <- rnorm(nrow(z),0,sd=0.1)
z$y <-   z[,sqrt((x))+noise ]
z$ydup <- z[, exp(log(sqrt(x)) + addGrowth(x, dur, asymp, b, l)  ) +noise ]
ggplot(z,aes(x,ydup))+geom_line()+labs(x="x",y="y")+
  geom_line(data=z[x %between% dur,],aes(x,y=y),col='red',linetype=2)+
  geom_vline(data=data.table(dur=dur),aes(xintercept=dur),linetype = 2,color='grey')+
  ggtitle(glue("Asymptote(a)=Log({asymp}), Offset(b)={b} Rate(l)={l}",asymp=round(asymp,3)))

I tried to model this using brms. Notice i set gr1_id to 90 in regions
outside [50,90]. Thats because the growth additive effect is near zero for those values.
I’ve also set lb=0 for the parameters b and l since they need to be
positive .


z <- z[, gr1_id := sapply(x-dur[1], function(s) if(s<0) 90 else s)]
m1 <- brm(
  bf(
    log(ydup) ~  slope*x + (log(a) - (log(a) * exp(-b * exp(-l * gr1_id)))),
    slope ~  1,
    a    ~ 1,
    b    ~ 1,
    l    ~ 1,
    nl = TRUE),
  data=z, family=gaussian(),
  prior = c(
    prior("normal(0,10)",nlpar='slope'),
    prior("normal(0, 10)",lb=0, nlpar = "a"),
    prior("normal(30,5)",lb=0, nlpar = "b"),
    prior("normal(0.3, 1)",lb=0, nlpar = "l")
  ),
  control = list(adapt_delta = 0.999))

 Family: gaussian
  Links: mu = identity; sigma = identity
Formula: log(ydup) ~ slope * x + (log(a) - (log(a) * exp(-b * exp(-l * gr1_id))))
         slope ~ 1
         a ~ 1
         b ~ 1
         l ~ 1
   Data: z (Number of observations: 150)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects:
                Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
slope_Intercept     0.01      0.00     0.01     0.01       2243 1.00
a_Intercept         4.37      0.19     4.01     4.74       1937 1.00
b_Intercept        27.53      5.39    16.49    38.11       1117 1.00
l_Intercept         0.03      0.00     0.02     0.03       1007 1.00

Family Specific Parameters:
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     0.21      0.01     0.18     0.23       2066 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Though b is estimated decently, the others ) are far removed.
My first reaction is that i miss specified the model, have i? How would i correct
it?

Thanks much

  • Operating System:
    Linux
  • brms Version:
    2.7

#2

I also want to point out that if i replace log(a) with A and set a normal(0,10) prior onA and estimate A, the estimation is much faster and the parameter estimate is the same.

They are all still wrong however … (except you could say b is somewhat close, but i’m surprised its not closer)


#3

I think the data generating process didn’t match the model (the usefulness of the bayesian workflow!). I changed it to a simpler model

addGrowth <- function(x,w,a,b,l){
  f <-sapply(x-w[1], function(s) if(s<0) 90 else s)
  log(a) - gompertz(f, log(a),b,l)
}
set.seed(3032)
z <- data.table(x=1:150)
dur <- c(50, 90);
asymp <- 1.35;b <-30;l <- 0.30
noise <- rnorm(nrow(z),0,sd=0.1)
z$y <-   z[, 10 +noise ]
z$ydup <- z[, 10 + addGrowth(x, dur, asymp, b, l)   +noise ]

and model it thus

z <- z[, gr1_id := sapply(x-dur[1], function(s) if(s<0) 90 else s)]
m1 <- brm(
  bf(
    ydup ~  slope + (A - (A * exp(-b * exp(-l * gr1_id)))),
    slope ~  1,
    A    ~ 1,
    b    ~ 1,
    l    ~ 1,
    nl = TRUE),
  data=z, family=gaussian(),
  prior = c(
    prior("normal(0,10)",nlpar='slope'),
    prior("normal(0, 10)",nlpar = "A"),
    prior("normal(30,5)",lb=0, nlpar = "b"),
    prior("normal(0.3, 3)",lb=0, nlpar = "l")
  ),
  control = list(adapt_delta = 0.999))

to get

Population-Level Effects:
                Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
slope_Intercept    10.01      0.01     9.99    10.02       3574 1.00
A_Intercept         0.28      0.03     0.22     0.35       2646 1.00
b_Intercept        30.15      4.90    20.65    39.79       3056 1.00
l_Intercept         0.31      0.04     0.24     0.40       2102 1.00

Family Specific Parameters:
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     0.10      0.01     0.09     0.11       3136 1.00

all of which is spot on.


#4

Glad you found a way to get it working. I am sorry that I did not reply earlier, but some questions tend to take so much time to dive into in order to be of real help.


#5

No worries at all! I’m so thankful for BRMS. I was getting into STAN and for an R user familiar with standard tools of frequentist linear modelling BRMS is such a wonderful introduction to the bayesian workflow.