Adding fixed effects and random effects to a nonlinear Stan model via brms

brms-nonlinear.r

# The data set and model are described in the *brms* vignette
library(brms)

url <- paste0("https://raw.githubusercontent.com/mages/diesunddas/master/Data/ClarkTriangle.csv")
loss <- read.csv(url)

set.seed(42)

# Generated a random continuous feature
loss$ramey <- runif(nrow(loss))

# Prior stays the same for the models below for illustration
nlprior <- c(prior(normal(5000, 1000), nlpar = "ult"),
             prior(normal(1, 2), nlpar = "omega"),
             prior(normal(45, 10), nlpar = "theta"))


# Nonlinear Example from the brms Vignette
nlform <- bf(cum ~ ult * (1 - exp(-(dev/theta)^omega)),
              ult ~ 1 + (1 + ramey | AY),
              omega ~ 1,
              theta ~ 1,
              nl = TRUE)

fit_loss <- brm(formula = nlform, data = loss,
                family = gaussian(), prior = nlprior,
                control = list(adapt_delta = 0.9))

# The nonlinear parameter *ult* is modeled by year (AY).
# The slope on the *ramey* variable varies by year (AY) when modeling the *ult* parameter.
nlform2 <- bf(cum ~ ult * (1 - exp(-(dev/theta)^omega)),
             ult ~ 1 + (1 + ramey | AY),
             omega ~ 1,
             theta ~ 1,
             nl = TRUE)

fit_loss2 <- brm(formula = nlform2, data = loss,
                 family = gaussian(), prior = nlprior,
                 control = list(adapt_delta = 0.9))


# Same as above, but this time the "ramey" variable is modeled at the population level.
# Fixed effect.
nlform3 <- bf(cum ~ ult * (1 - exp(-(dev/theta)^omega)),
             ult ~ 1 + ramey + (1 | AY),
             omega ~ 1,
             theta ~ 1,
             nl = TRUE)


fit_loss3 <- brm(formula = nlform3, data = loss,
                 family = gaussian(), prior = nlprior,
                 control = list(adapt_delta = 0.9))


# To see the marginal effect of the "ramey" fixed effect in the 3rd model, run:
marginal_effects(fit_loss3)



# Same thing, but this time fixed and random effect on "ramey"
nlform4 <- bf(cum ~ ult * (1 - exp(-(dev/theta)^omega)),
              ult ~ 1 + ramey + (1 + ramey | AY),
              omega ~ 1,
              theta ~ 1,
              nl = TRUE)


fit_loss4 <- brm(formula = nlform4, data = loss,
                 family = gaussian(), prior = nlprior,
                 control = list(adapt_delta = 0.9))


# To see the marginal effect of the "ramey" fixed effect in the 4th model, run:
marginal_effects(fit_loss4)