 # 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)
f <-sapply(x-w, 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, 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 on`A` 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, 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, 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.