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