Model 10.16 about statistical_rethinking_with_brms_ggplot2_and_the_tidyverse HOT 8 CLOSED

Sure thing — this is just from the stancode function in brms.

// generated with brms 2.12.0
functions {
}
data {
  int<lower=1> N;  // number of observations
  int<lower=2> ncat;  // number of categories
  int Y[N];  // response variable
  int<lower=1> K_b;  // number of population-level effects
  matrix[N, K_b] X_b;  // population-level design matrix
  // covariate vectors for non-linear functions
  int C_mu2_1[N];
  // covariate vectors for non-linear functions
  int C_mu3_1[N];
  int prior_only;  // should the likelihood be ignored?
}
transformed data {
}
parameters {
  vector[K_b] b_b;  // population-level effects
}
transformed parameters {
}
model {
  // initialize linear predictor term
  vector[N] nlp_b = X_b * b_b;
  // initialize non-linear predictor term
  vector[N] mu2;
  // initialize non-linear predictor term
  vector[N] mu3;
  // linear predictor matrix
  vector[ncat] mu[N];
  for (n in 1:N) {
    // compute non-linear predictor values
    mu2[n] = nlp_b[n] * C_mu2_1[n];
  }
  for (n in 1:N) {
    // compute non-linear predictor values
    mu3[n] = nlp_b[n] * C_mu3_1[n];
  }
  for (n in 1:N) {
    mu[n] = [0, mu2[n], mu3[n]]';
  }
  // priors including all constants
  target += normal_lpdf(b_b | 0, 1);
  // likelihood including all constants
  if (!prior_only) {
    for (n in 1:N) {
      target += categorical_logit_lpmf(Y[n] | mu[n]);
    }
  }
}
generated quantities {
}

from statistical_rethinking_with_brms_ggplot2_and_the_tidyverse.

This requires the new custom_family feature of brms. Maybe you find some inspiration in the related vignette: https://github.com/paul-buerkner/brms/blob/master/vignettes/brms_customfamilies.Rmd

from statistical_rethinking_with_brms_ggplot2_and_the_tidyverse.

Hi @ASKurz ,

First off, thanks so much for all the work you've put into your bookdown projects. I've found them tremendously helpful!

I'm still a bit naive when it comes to Bayesian statistics, but having now gone through McElreath's book, it also has been bugging me that I can't figure out how to replicate his model in brms, which I love to use whenever possible. (Thank you, @paul-buerkner for designing such a wonderful package and answering so many questions!) I was wondering whether some sort of non-linear model like the one specified below might work?

nl_model <-
  brm(data = dat,
      family = categorical(link = logit),
      bf(career ~ b * income, b ~ 1, nl = TRUE),
      prior = prior(normal(0, 5), nlpar = "b"),
      iter = 2500, warmup = 500, cores = 2, chains = 2,
      seed = 10)

That is, "income" here is just a variable that's equal to "career", as in the textbook. We could've just as well specified the model as "career ~ b * career"; the point is to use the nonlinear modeling feature of brms to force it to fit just a single "b" parameter. When I do this, I get an estimate of around b = .72, which is somewhat close to the right answer, but clearly biased upward. Interestingly, McElreath's own code gives a biased result in the opposite direction, but that doesn't use HMC.

Looking over the Stan code of the model I built, I couldn't figure out how it might differ from McElreath's model. I suspect something is wrong — especially if @paul-buerkner is saying you need to use a custom_family — but I can't figure out why it would be wrong.

All of this is to say that perhaps this is a workaround to the problem, but I'm not sure it's correct :)

from statistical_rethinking_with_brms_ggplot2_and_the_tidyverse.

Thanks for your kind words and your interest in the problem. Though I’m not fluent in Stan code, it might help others if you shared the Stan code for your model.

from statistical_rethinking_with_brms_ggplot2_and_the_tidyverse.

Hi All,

I've been struggling with McElreath's R code 10.57 for quite some time, and have (reluctantly) come to the conclusion that it's not quite right. In particular, to have just one "b" parameter, the map() functions' linear models need to be "centered" on whichever career is set to zero. Doing so preserves the constant relationships among the careers. (That is, the distances between the scores in the linear space.) If you don't do this, then McElreath's statement, "In the example above, I chose the the first event type, the first career. If you choose another, you'll get different estimates, but the same predictions" doesn't hold. (It does hold, though, if you use N-1 parameters, vs. just one. In that approach, you're establishing independent distances from zero to the other scores.)

I am fairly new at all this (both R and Bayesian statistical inferencing), and I don't know bmrs, so apologies if this is off-topic. But, it seems to me that "replicating" McElreath's code isn't going to work.

Here's code that shows that we do get the same results (probabilities), as long as we "center" around the zeroed career.

In the code below, you can set refCareer to 1, 2, or 3, and you get the same probabilities from your model, regardless. (I have another code block that uses N-1 [=2, in this case] parameters, and that doesn't require "centering." Happy to share if anyone is interested.)

library(rethinking)
library(reshape2)
library(ggplot2)

refCareer = 1  # Set to 1, 2, or 3
set.seed(1000)

## R code 10.56
# simulate career choices among 500 individuals
N <- 10000             # number of individuals
income <- 1:3        # expected income of each career
score <- 0.5*income  # scores for each career, based on income
# next line converts scores to probabilities
p <- softmax(score[1],score[2],score[3])

# now simulate choice
# outcome career holds event type values, not counts
career <- rep(NA,N)  # empty vector of choices for each individual
# sample chosen career for each individual
for ( i in 1:N ) career[i] <- sample( 1:3 , size=1 , prob=p )

# fit the model, using dcategorical and softmax link
if (refCareer == 1) {
  m10.16 <- map(
    alist(
      career ~ dcategorical( softmax(0,s2,s3) ),
      s2 <- b*(2-refCareer),    # linear model for event type 2
      s3 <- b*(3-refCareer),    # linear model for event type 3
      b ~ dnorm(0,5), 
    ) ,
    data=list(career=career) )
  # start=list(b=0) )
} else if (refCareer == 2) {
  m10.16 <- map(
    alist(
      career ~ dcategorical( softmax(s1,0,s3) ),
      s1 <- b*(1-refCareer),    # linear model for event type 1
      s3 <- b*(3-refCareer),    # linear model for event type 3
      b ~ dnorm(0,5)
    ) ,
    data=list(career=career) )
} else { # refCareer == 3
  m10.16 <- map(
    alist(
      career ~ dcategorical( softmax(s1,s2,0) ),
      s1 <- b*(1-refCareer),    # linear model for event type 1
      s2 <- b*(2-refCareer),    # linear model for event type 2
      b ~ dnorm(0,5)
    ) ,
    data=list(career=career) )
}

b = coef(m10.16)[1]

# Fn to convert linear params to probabilities (same as softmax())
softmaxFn = function(xVec) { # multinomial equiv to logisticFn
  pVec = exp(xVec) / sum(exp(xVec))
  return(pVec)
}

# Fn to convert probabilities to linear params
c = 0  # Any value works; doesn't matter! 
dcategoricalFn = function(p, c) {
  linearParams = log(p) + c
  return(linearParams)
}

careerProbFn = function(b) {
  s1 = b*1
  s2 = b*2
  s3 = b*3
  sVec = c(s1, s2, s3)
  careerProb = softmaxFn(sVec)
  return(careerProb)
}

careerProbMu = careerProbFn(b)
df = data.frame(income = c(1, 2, 3), genProb = p, 
                careerProb = careerProbMu)
dfMelt = melt(df, id.vars = "income")
plt = ggplot(data = dfMelt, aes(x = income, y = value, color = variable)) +
  geom_point() +
  geom_line() +
  ylim(0.0, 1.0) + 
  theme_minimal()
plot(plt)

from statistical_rethinking_with_brms_ggplot2_and_the_tidyverse.

Hey @ssickels,

Thanks for sharing your thoughts and work, on this. Updates to my translation of Statistical Rethinking (1st ed) are going to be on the backburner for a while, though I do hope to release a minor update by the end of the year. That update should contain revisions to this section and I will review your comment in detail when working on that section.

Heads up, though: @adambear91 and I hammered out this issue quite a bit in an email exchange, a few months ago. I [tentatively] think we solved the problem. Our efforts were focused in the corresponding problems in the second edition of Statistical Rethinking. In short, the solution had to do with the non-linear syntax. You can check out our attempts in Section 11.3 of my ebook.

from statistical_rethinking_with_brms_ggplot2_and_the_tidyverse.

Hi @ASKurz,

Many thanks for the note and for the link to your newest work. (So impressive!) And looks like I need to get the 2nd ed!

I'm still playing around with the (1st edition's) m10.56. I've now got it "working" with McElreath's map2stan (as well as with his map(), as in the code above), albeit with scary warning messages. But despite the warnings, it looks like it gets the right model. I don't (yet) know tidyverse -- nor stan (or bmrs). But it looks to me like your model has three learned parameters (b, and a[1] and a[2]). Whereas (I believe) my model gets the good results with just one (b). I'm not sure that my approach (shifting all the predictor variables to maintain their distances from the reference category at 0) is fully "general," though. I'm trying to get it to work when when the income predictor values aren't evenly spaced (as in your code, and, I presume, the 2nd ed's example). But no luck yet. If I make any particular progress on these fronts, I'll let you know. And, I'll dig into tidyverse, which has been on my "to-do list" for a while now. And into your code (and bmrs). Cheers -- and thanks again!

Edit: I got it working with unevenly-spaced incomes (with map(), and map2stan(), albeit with the scary warnings). So I am getting the idea that it's a good general solution, with just a single learned param (that is, just b).

from statistical_rethinking_with_brms_ggplot2_and_the_tidyverse.

This is the solution.

b10.16 <-
  brm(data = list(career = career), 
      family = categorical(link = logit),
      bf(career ~ 1,
         nlf(mu2 ~ b * 2),
         nlf(mu3 ~ b * 3),
         b ~ 1),
      prior(normal(0, 5), class = b, nlpar = b),
      iter = 2500, warmup = 500, cores = 2, chains = 2,
      seed = 10)

To confirm, run McElreath's code.

m10.16 <- map(
  alist(career ~ dcategorical( softmax(0,s2,s3) ), 
        s2 <- b*2, # linear model for event type 2 
        s3 <- b*3, # linear model for event type 3 
        b ~ dnorm(0,5)
  ), data=list(career=career) )

Both return b values of about 0.34. When you transform that value out of the log-odds metric, you get about 0.59. This is an approximation of one of the original values used to generate the data. It's an odd model, but we finally have the brms solution. I'll flesh this all out in the next release.

from statistical_rethinking_with_brms_ggplot2_and_the_tidyverse.