Brms predictions inconsistent with custom model (likely user error)

General
1

Hi all,

This is probably a bit embarassing, but I cannot for the life of me figure out why the posterior preditions from a brmsfit object are different (better) than my custom negative binomial model. The model example is a simple regression with two harmonic terms. I have tried generating posterior predictions with both neg_binomial_2_rng() with the inverse-link function and neg_binomial_2_log_rng() with the linear predictor. The posterior distributions for the parameters are almost identical, so the model (despite being extremely inefficient) is fine. It’s just the predictions that are confusing me a bit. The code for my custom model is below, and I am attaching my R code and a sample of the data. Pardon the clunkiness of the stan code, the real model I am working is much more complex, and I wanted to firgure out the issues of the toy example before moving to correct the real one, so I just stripped down the more complex versions so it was easier to fix the other model structures. I’m 99% positive it’s a user error issue, I just don’t know what it is.

Thanks!

data {
  int<lower=0> M_Y;                                    // # Months x # Years of data
  int<lower = 0> Trips[M_Y];                           // Response: Trips in month m in year y
  matrix[M_Y,2] X;                                     // Predictors of max effort: column 1 = sine, 2 = cosine
  vector[M_Y] Pop;                                     // Florida population
  vector[M_Y] CPUE;                                    // CPUE for each species s in year y
}

parameters {
  real<lower = 0, upper = 100> phi;                    // Dispersion parameters
  vector[2] beta;                                      // Effects of seasonality (harmonic terms)
  real eta;                                            // Minimum effort (intercept)
  real xi;                                             // Biomass
  real tau;                                            // Florida population effect

}

transformed parameters {
  vector[M_Y] mu;                                      // Mean structure for counts (log scale)
  vector[M_Y] alpha;                                  
  alpha = X[,1]*beta[1] + X[,2]*beta[2] +xi*CPUE + tau*Pop;
  mu = alpha+eta;
}

model {
    for (my in 1:M_Y){
      Trips[my] ~ neg_binomial_2_log(mu[my], phi);
    }
}

generated quantities {
  vector[M_Y] pred;
  vector[M_Y] log_lik;
  for (my in 1:M_Y){
      log_lik[my] = neg_binomial_2_log_lpmf(Trips[my] | mu[my], phi);
      pred[my] = neg_binomial_2_rng(exp(mu[my]), phi);
    }
}

R code:

#------------------------------------------------------------------------------#
##################################### Libraries ################################
#------------------------------------------------------------------------------#
## Syntax packages
suppressMessages(library(dplyr))
suppressMessages(library(tidyr))
suppressMessages(library(tidyverse))
suppressMessages(library(lubridate))

## Visualization packages
suppressMessages(library(ggplot2))
options(width = 80)

### Set ggplot theme for plotting
theme_set(theme_bw(base_family = 'serif')) 


## Modeling
suppressMessages(library(rstan))
suppressMessages(library(brms))

#------------------------------------------------------------------------------#
##################################### Housekeeping #############################
#------------------------------------------------------------------------------#
## Clear out old files in R
rm(list=ls(all=TRUE)) 

## Set working directory
setwd("C:\\Users\\ichal\\Documents\\Hyman gag grouper models")

## User-defined functions
`%nin%` <- Negate(`%in%`)

#------------------------------------------------------------------------------#
##################################### Formatting ###############################
#------------------------------------------------------------------------------#
## Load data
Data <- read.csv("Snapper.csv")
#Datalist for model
DataList <- list(
  "M_Y" = nrow(Data),
  "Pop" = Data$pop,
  "Trips" = ceiling(Data$VS),
  "X" = Data[,2:3],
  "CPUE" = Data$SSB
)

## My model
My_Model <- stan('Example.stan', data = DataList, chains = 1, iter = 2000, control = list(adapt_delta = 0.99)) 
BRM_Model <- brm(ceiling(VS)~sin1+cos1+pop + SSB, data = Data, family = negbinomial())

BRM_preds <- predict(BRM_Model, Data, type = 'response')
Data$Median_BRM <- BRM_preds[,1]
Data$Upper_BRM  <- BRM_preds[,4]
Data$Lower_BRM  <- BRM_preds[,3]


Model_pred <- extract(My_Model)
Data$Median_custom <- apply(Model_pred$pred, 2, function(x){quantile(x, 0.5)})
Data$Lower_custom <- apply(Model_pred$pred, 2, function(x){quantile(x, 0.025)})
Data$Upper_custom <- apply(Model_pred$pred, 2, function(x){quantile(x, 0.975)})
Data$Time <- 1:276

ggplot(Data)+
  geom_point(aes(x = Time, y = VS), col = 1)+
  geom_line(aes(x = Time, y = Median_BRM, color = 'a'))+
  geom_line(aes(x = Time, y = Median_custom, color = 'b'))+
  scale_colour_manual(name = 'Model', 
                      values = c('a'='forestgreen','b'=2), 
                      labels = c('brms','custom'))
2

Help.R (2.5 KB)
Snapper.csv (19.1 KB)
Example.stan (1.5 KB)

3

Howdy! I haven’t had a chance to run your example, but could you define what you mean by “different (better)”? Also, it doesn’t appear that you have set any priors for either the Stan or brms model. I believe when the user doesn’t set priors, that for Stan, they are all flat, but for brms, only the ‘population-level’ effects have a flat prior, while the Intercept and the shape parameter have custom default priors. So, your Stan and brms models are not quite the same. You can use stancode(BRM_Model) to see the Stan code from brms to compare to your Stan code.

4

Sure! Thanks for taking a look. What I meant by “different (better)” is that the fit of the posterior predictions from the brms model are more closely aligned with the obsereved data than the ones from my hand-written one, which underestimates the response variable. I don’t think it is a prior issue. Although I removed the priors for this example to make it as simple as possible, I usually add diffuse priors and had done so in previous iterations. Moreover, the fact that the posterior distributions of the fixed effects are nearly identical among the two models provides evidence that this is probably related to the differences between predict_brms() and the neg_binomial_2_rng().

Since the posterior distributions are the same, worst-case scenario is that I just throw the scans into rnbinom() to get posterior predictive values, but figuring out what is wrong with neg_binomial_2_rng() would save me a little time.

-Challen

5

image
image1408×753 10.8 KB

This is the plot that I got, using your above code and data. Did you get the same plot with a Stan model that used neg_binomial_2_log_rng?

6

Yes, exactly the same. Any idea what I’m doing wrong here?

7

You are comparing means from brms to medians from Stan.

Try this instead:

BRM_preds <- predict(BRM_Model, Data, type = 'response', summary=FALSE)
Data$Median_BRM <- apply(BRM_preds, 2, function(x){quantile(x, 0.5)})


Model_pred <- extract(My_Model)
Data$Median_custom <- apply(Model_pred$pred, 2, function(x){quantile(x, 0.5)})

Data$Time <- 1:276

ggplot(Data)+
  geom_point(aes(x = Time, y = VS), col = 1)+
  geom_line(aes(x = Time, y = Median_BRM, color = 'a'))+
  geom_line(aes(x = Time, y = Median_custom, color = 'b'))+
  scale_colour_manual(name = 'Model', 
                      values = c('a'='forestgreen','b'=2), 
                      labels = c('brms','custom'))

image
image1408×753 8.94 KB

2 Likes
10

Yep, that worked. I feel a little silly and should have read more carefully. Thanks!

2 Likes