Advice/Peer Review: Compare Robust Linear Model to Piecewise Linear Model w/ Random Change Point

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"))

image

# 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
1 Like

I know this is a forum for stan/brms, but I just want to briefly mention another approach using mcp as a frontend to JAGS, if nothing else, just to verify the results you infer via brms.

# Fit models with and without change points
library(mcp)
model_full = list(
  log_c ~ 1 + log_q,  # Intercept and slope in the first segment
  ~ 0 + log_q  # Joined slope in the second segment
)
model_null = list(
  log_c ~ 1 + log_q  # Just one segment with intercept+slope
)

fit_full = mcp(model_full, dat)
fit_null = mcp(model_null, data)

# Inspect results
summary(fit_full)
plot_pars(fit_full)
plot(fit_full, q_fit = TRUE)

image

The density on x is the posterior for the change point and the dashed red lines is the 95% credible interval. mcp tries to mimic brms in many ways, e.g., the default priors. But it includes default priors for the change point as well (read more in the preprint and this guide on setting priors in mcp).

Comparing the models using LOO/WAIC (read more here):

loo::loo_compare(loo(fit_full), loo(fit_null))
loo::loo_compare(waic(fit_full), waic(fit_null))
2 Likes

I noticed both of your variables appear to have been log transformed. If your criterion c is naturally non-negative and continuous, you might consider using one of the likelihood functions designed for that purpose (e.g., exponential, gamma, log-normal, Weibull). For example, here’s a version of your model where log_q predicts c, using the Weibull likelihood.

# load
library(tidyverse)
library(brms)

# wrangle
dat <- dat %>% 
  mutate(q = exp(log_q),
         c = exp(log_c))


# fit
fit_weibull <-
    brm(
      data = dat,
      family = weibull,
      c ~ 1 + log_q,
      cores = 4, seed = 4
    )

# plot
conditional_effects(fit_weibull) %>% 
  plot(points = TRUE)

Screen Shot 2021-12-15 at 9.28.55 AM

The fit isn’t perfect, but this model’s a big win from a parsimony perspective.

2 Likes

Is there a way to limit the x-range of the change point in the full model? My data will only have 1 change point, if any. For example, I think pp below is how a colleague limits the range in a JAGS model string:

# Specify the data:
  data {
    Ntotal <- length(y)
    range <- max(x)-min(x)
    low <- min(x) + 0.15*range
    hi <- max(x) - 0.20*range 
  }
  # Specify the model:
  model {
    #Likelihood
    for (i in 1:Ntotal){
      y[i] ~ dnorm(yMean[i], tau.y[K[i]])
      yMean[i]<-beta0 + 
         (beta1 + delta*step(x[i]-x.change[i]) )*(x[i]-x.change[i])
      x.change[i] <- (1-pp)*low+pp*hi   #calculation of flow threshold
      K[i]<- step(x[i]-x.change[i])+1
      }
      beta0.true <- beta0 - x.change[1]*beta1
      beta1.star <- beta1 + delta
      
    # Priors Vague
    for ( K in 1:2 ) {
      tau.y[K] <- pow(sigma.y[K],-2)
      sigma.y[K] ~ dunif(0,4)
      }
    delta ~ dnorm(0,0.0001)
    pp ~ dbeta(1,1)	#this is uniform beta distrib. when a=1, b=1 the mean and mode are 0.5 per Kruschke, p.129
    beta0 ~ dnorm(0,0.0001)
    beta1 ~ dnorm(0,0.0001)
  }

For more back-and-forth regarding mcp , I think it’ll be better to raise a discussion in the mcp forum. I only dared respond here, because the discussions on how to set priors for change points in mcp may be readily translated to stan / brms .

But briefly (read more here), you can truncate priors in JAGS using:

  • An inherently truncated distribution, such as prior = list(cp_1 = "dunif(0, 0.5)") or
  • Imposing truncation on other distributions, e,g, prior = list(cp_1 = "dnorm(0.25, 0.5) T(0, 0.5)")
1 Like