Translating robust logistic regression from rstan to brms

rstan

#1

I am interested in fitting a “robust” logistic regression model from Doing Bayesian Data Analysis (2nd Edition). This model treats data as a mixture of two sources: (1) guessing and (2) a logistic function of predictors. For guessing, the y values are thought to come from y ~ Bernoulli(m=1/2).

Using some data I found in this vignette, I think I’ve been able to implement the model in a simple case using rstan:

library("rstan")

# Prepare data 
url <- "http://stat.columbia.edu/~gelman/arm/examples/arsenic/wells.dat"
wells <- read.table(url)
X <- model.matrix(~ arsenic, wells)
standata <- list(y = wells$switch, arsenic = wells$arsenic, N = nrow(wells))

# Model
guess.model <- '
data {                          
  int<lower=0> N;                
  int<lower=0,upper=1> y[N];     
  vector[N] arsenic;            
}
parameters {
  real alpha;                   
  real b_arsenic;                
  real<lower=0,upper=1> guess;
}
transformed parameters { 
  vector[N] mu;
  for ( i in 1:N) {
    mu[i] = .5*guess + (1-guess)*inv_logit(alpha + b_arsenic * arsenic[i]);
  }
}
model {
  y ~ bernoulli(mu);
  alpha ~ normal(0,1);         
  b_arsenic ~ normal(0,1);     
  guess ~ beta(1,9);
}
'

# Fit model
fit <- stan(model_code = guess.model, data = standata)
print(fit, pars = c("alpha","b_arsenic","guess"))

However, I’m unsure how to implement the guessing parameter in brms. Could anyone point me in the right direction?

Operating System: Windows 10
Interface Version: 3.2.1


#2

You can learn more about that in vignette("brms_nonlinear")


#3

Many thanks to jrrae for posing this question. I’ve been working through Kruschke’s book and was initially baffled by this model, too (chapter 21.3). In case any others are curious, I was able to follow Paul’s lead and fit the model like so.

fit <-
  brm(data = my_data, 
      family = binomial(link = "identity"),
      bf(male | trials(1) ~ a * .5 + (1 - a) * 1 / (1 + exp(-1 * (b0 + b1 * weight_z))),
         a + b0 + b1  ~ 1,
         nl = TRUE),
      prior = c(prior(normal(0, 4), nlpar = "b0"),
                prior(normal(0, 4), nlpar = "b1"),
                prior(beta(1, 9),   nlpar = "a")))

Much like Paul cautioned in his vignette("brms_nonlinear"), the model yielded a number of divergent transitions which might well be alleviated with better priors. But the results match Kruschke’s pretty well.

Note. My model is based on the data and parameters from Kruschke’s text, not those in the vignette jrrae cited, above.