predict, pp_check, fitted do not compute predictions correctly for nonlinear model

[femiguez]

I have an example where if I compute predictions manually I get the correct values, but it seems that predict, pp_check and fitted do not. Here is a reproducible example:

nonlinear model

quadp3 <- function(x, a, b, lc){
  ans <- (x < -0.5 * b/(-exp(lc))) * (a + b * x - exp(lc) * x * x) + (x >= -0.5 * b/(-exp(lc))) * (a + -(b * b)/(4 * -exp(lc)))
  ans
}

#Simple dataset
x <- seq(0, 0.4, length.out = 30)
y <- quadp3(x, 6, 60, 4.7) + rnorm(length(x), sd = 0.5)
dat <- data.frame(x = x, y = y)

ggplot(dat, aes(x, y)) + 
  geom_point()

# brms fit
prs1 <- prior(normal(6, 3), nlpar = "a", coef = "Intercept") +
  prior(normal(60, 20), nlpar = "b", coef = "Intercept") + 
  prior(normal(4.7, 1), nlpar = "lc", coef = "Intercept")

bf1 <- bf(y ~ (x < -0.5 * b/(-exp(lc))) * (a + b * x - exp(lc) * x * x) + 
            (x >= -0.5 * b/(-exp(lc))) * (a + -(b * b)/(4 * -exp(lc))), 
          a + b + lc ~ 1,
          nl = TRUE)

brm1 <- brm(formula = bf1,
            prior = prs1,
            iter = 3000, 
            chains = 3,
            cores = 3,
            seed = 9,
            data = dat)
fixef(brm1)

This works fine and produces sensible estimates, but predict, pp_check and fitted do not work as expected.

post_samp <- posterior_samples(brm1)

# Computing predictions manually
yprds <- matrix(ncol = 12, nrow = length(x))

k <- 1
for(i in seq(1, nrow(post_samp), by = 400)){
       yprds[,k] <- quadp3(x, post_samp[i,1], post_samp[i,2], post_samp[i,3])
       k <- k + 1
}

psum <- posterior_summary(t(yprds))

ggplot(dat, aes(x, y)) + 
  geom_point() + 
  geom_line(aes(y = fttd, color = "fitted")) +
  geom_line(aes(y = psum[,1], color = "manual prediction"))

Fitted is flat while the "manual" prediction behaves as expected.

I apologize if this is known, I might be missing something or it has been reported elsewhere.

[paul-buerkner]

Thanks. Will investigate.

[paul-buerkner]

Found it. In some formula handling, white spaces where removed which accidentally transformed (x < -0.5 * b/(-exp(lc))) into (x <- 0.5 * b/(-exp(lc))), that is, combining < and -. I will fix this today.

[paul-buerkner]

Should now be fixed. Thanks again for reporting this issue!

[femiguez]

Thanks so much! Will test it! :)