I want to compare a robust linear model with no change point to a piecewise linear model with a random change point (using information criteria). I feel fairly confident I have specified the robust linear model correctly, but I am not confident I have specified the piecewise model and its priors correctly. I have received errors after running the piecewise model (with more iterations) that suggest I should adjust alpha and maximum tree depth. I welcome advice on that, but I mostly want feedback on model specification and priors. I am completely new to Bayesian analysis and piecewise models and my code is naively based off of this post: Piecewise Linear Mixed Models With a Random Change Point - #9 by lindeloev.
# Load packages
library("brms")
# Data
dat <-
structure(list(log_q = c(0, -0.25, 0.96, 0.22, 0.37, 1.2, 0.99,
0.87, 0.17, 0.52, 1.21, 0.6, 0.91, -0.04, 0.99, 0.59, 1.21, -0.07,
0.38, 1.35, 1.29, 0.76, -0.09, 1.62, 1.01, 1.36, -0.17, -0.19,
1, -0.34, 0.7, 0.97, 1.58, 0.98, 1.7, -0.03, 0.65, 0.87, 1.11,
1.19, 0.39, 1.16, 0.14, 0.37, 1.02, -0.17, -0.13, 0.73, 1.29,
1.16, 0.97, 0.24, 1.26, 0.9, 1.08, 0.63, 0.38, 1.55, -0.06, 0.33,
1.14, 1.1, -0.19, -0.34, 0.33, 0.91, 0.85, -0.4, 0.03, 0.66,
-0.43, 0.98, 0.08, 0.68, 1.1), log_c = c(0.84, 0.58, 1.2, 0.99,
0.96, 2.09, 1.79, 1.47, 0.53, 1.07, 2.23, 0.66, 1.2, 0.38, 1.34,
0.88, 2.59, 0.97, 0.86, 2.51, 2.03, 1.52, 0.79, 2.5, 1.72, 2.05,
0.31, 0.2, 1.81, 1.51, 1.8, 1.58, 2.36, 1.41, 2.28, 0.26, 1.62,
1.63, 1.81, 1.87, 1.78, 1.79, 0.54, 0.29, 2.55, 0.86, 0.9, 0.98,
2.64, 1.68, 1.37, 0.82, 2, 1.76, 1.48, 1.06, 0.6, 2.03, 1, 0.69,
1.51, 2.43, 0.75, 0.32, 0.94, 1.37, 1.99, 0.68, 0.49, 1.59, 0.26,
1.9, 0.47, 1.28, 2.17)), row.names = c(NA, -75L), class = c("tbl_df",
"tbl", "data.frame"))
# Calculate necessary values we'll use below
sd_x <- sd(dat$log_q)
sd_y <- sd(dat$log_c)
m_x <- mean(dat$log_q)
m_y <- mean(dat$log_c)
min_x <- min(dat$log_q)
max_x <- max(dat$log_q)
# The following are for specifying the vague priors
stanvars <-
brms::stanvar(sd_y, name = "sd_y") +
brms::stanvar(10*abs(m_x*sd_y/sd_x), name = "Intercept_sigma") +
brms::stanvar(10*abs(sd_y/sd_x), name = "slope1_sigma") +
brms::stanvar(1/29, name = "one_over_twentynine") +
brms::stanvar(min_x, name = "min_x") +
brms::stanvar(max_x, name = "max_x")
# Robust linear regression (no random change point)
fit_lm <-
brms::brm(
data = dat,
family = student,
log_c ~ 1 + log_q,
prior = c(
prior(normal(0, Intercept_sigma), class = Intercept),
prior(normal(0, slope1_sigma), class = b),
prior(normal(0, sd_y), class = sigma),
prior(exponential(one_over_twentynine), class = nu)
),
stanvars = stanvars
)
# Piecewise linear model w/ random change point
fit_pw <-
brm(
data = dat,
family = gaussian(),
formula =
bf(log_c ~
Intercept +
slope1 * log_q * step(change - log_q) +
(slope1 * change + slope2 * (log_q - change)) * step(log_q - change),
Intercept + slope1 + slope2 ~ 1,
change ~ 1, # Should this be included in the above statement?
nl = TRUE),
prior = c(
prior(normal(0, Intercept_sigma), nlpar = Intercept),
prior(normal(0, slope1_sigma), nlpar = "slope1"),
prior(normal(0, slope1_sigma), nlpar = "slope2"),
prior(uniform(min_x, max_x), lb = min_x, ub = max_x, nlpar = "change") # Is this the best vague prior?
# Do I need to specify a prior for sigma? e.g.:
# prior(normal(0, sd_y), class = sigma),
),
stanvars = stanvars
)
- Operating System: macOS Big Sur Version 11.1
- brms Version: 2.16.3