Disrepancies in output of conditional_effects() between Windows and Linux for the same model

[Maxime]

Hello all,

I’ve been scratching my head with a strange issue I am facing with brms.

I produced plots from a brms model under Windows using conditional_effects(). I migrated my project to Linux and ensured that dependencies were the same using renv. I used the same package and R versions (see attached files).

We can see the differences in this figure :

• The left part shows the predictions under Windows for two models
• The right part shows the predictions under Linux Pop OS for the same two models

I did not refit the models under the different operating systems, I just imported the .rds files and computed the predictions (see R code below).

I’d love to know why this occurs and how we can correct it. This is for an upcoming scientific publication and I want the process to be as reproducible as possible. Running the code in any OS should result in the same outputs since they are from the same model, no?

Please find attached two files in the thread so that people can go and try reproducing my problem.

• sessionInfo() under Windows
• sessionInfo() under Linux (Pop OS)

Here is a link to the model files to reproduce the results. You can download GAMM-II.rds and GAMM-V.rds : OSF | Model outputs from: Prey movement shapes the acquisition of predator expertise in a virtual bi-trophic system

Here is the link to the raw data : GitHub - quantitative-ecologist/predator-expertise: Data and code for Fraser Franco et al. XXXX

Here is the code to generate the plots.

# ==========================================================================

#                           Plot GAMM models

# ==========================================================================



# ==========================================================================
# 1. Prepare script
# ==========================================================================



# Load libraries and model -------------------------------------------------

options(mc.cores = parallel::detectCores())

library(parallel)
library(brms)
library(data.table)
library(ggplot2)
library(ggpubr)
library(viridis)

path <- file.path(getwd(), "outputs", "outputs_models")

mod1 <- readRDS(file.path(path, "GAMM-II.rds"))
mod2 <- readRDS(file.path(path, "GAMM-V.rds"))



# Load the data ------------------------------------------------------------

data <- fread(
  "./data/FraserFrancoetal2023-data.csv",
  select = c("predator_id",
             "game_duration",
             "pred_speed",
             "prey_avg_speed",
             "prey_avg_rank",
             "cumul_xp_pred",
             "total_xp_pred",
             "hunting_success")
)

data[, predator_id := as.factor(predator_id)]



# Post processing preparations for custom family ---------------------------

expose_functions(mod1, vectorize = TRUE)

# Define the log likelihood function
log_lik_beta_binomial2 <- function(i, prep) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  phi <- brms::get_dpar(prep, "phi", i = i)
  trials <- prep$data$vint1[i]
  y <- prep$data$Y[i]
  beta_binomial2_lpmf(y, mu, phi, trials)
}

# Define function for posterior_predict
posterior_predict_beta_binomial2 <- function(i, prep, ...) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  phi <- brms::get_dpar(prep, "phi", i = i)
  trials <- prep$data$vint1[i]
  beta_binomial2_rng(mu, phi, trials)
}

# Define function for posterior_epred
posterior_epred_beta_binomial2 <- function(prep) {
  mu <- brms::get_dpar(prep, "mu")
  trials <- prep$data$vint1
  trials <- matrix(trials, nrow = nrow(mu), ncol = ncol(mu), byrow = TRUE)
  mu * trials
}

# ==========================================================================
# ==========================================================================





# ==========================================================================
# 2. Prepare the data for the plots
# ==========================================================================



# Prepare the plots --------------------------------------------------------

# Group-level smooths model 1
tab1_a <- conditional_effects(
  mod1, method = "fitted",
  effects = "Zcumul_xp:predator_id",
  robust = TRUE, re_formula = NULL
)
tab1_a <- data.table(tab1_a[[1]])

# Cumulative XP global trend model 1
tab1_b <- conditional_effects(
  mod1, method = "fitted",
  robust = TRUE, re_formula = NULL,
  effects = "Zcumul_xp",
  conditions = data.frame(predator_id = NA)
)
tab1_b <- data.table(tab1_b[[1]])

# Group-level smooths model 2
tab2_a <- conditional_effects(
  mod2, method = "fitted",
  effects = "Zcumul_xp:predator_id",
  robust = TRUE, re_formula = NULL
)
tab2_a <- data.table(tab2_a[[1]])

# Cumulative XP global trend model 2
tab2_b <- conditional_effects(
  mod2, method = "fitted",
  robust = TRUE, re_formula = NULL,
  effects = "Zcumul_xp",
  conditions = data.frame(predator_id = NA)
)
tab2_b <- data.table(tab2_b[[1]])



# Transform values --------------------------------------------------------

# Back transform x-axis values
sequence <- (seq(0, 500, 100) - mean(data$cumul_xp_pred))
standev <- sd(data$cumul_xp_pred)
scaled_breaks <- sequence / standev

# List the tables
tables <- list(
  tab1_a, tab1_b,
  tab2_a, tab2_b
)
names(tables) <- c(
  "tab1_a", "tab1_b",
  "tab2_a", "tab2_b"
)

# Function to apply transformation
# Computes non standardized cumulative XP
func <- function(x) {
  x[, cumul_xp := (Zcumul_xp * standev) + mean(data$cumul_xp_pred)]
}

# Apply function
lapply(tables, func)



# Cut fitted values based on player XP ------------------------------------

# Extract player IDs with their total XP from the original data
xp <- unique(data[, .(predator_id, total_xp_pred)])

# Merge the two tables adding the total XP
tab1_a <- merge(tab1_a, xp, by = "predator_id")
tab2_a <- merge(tab2_a, xp, by = "predator_id")

# Cut all matches where fitted values are above total XP
tab1_a <- tab1_a[cumul_xp <= total_xp_pred, ]
tab2_a <- tab2_a[cumul_xp <= total_xp_pred, ]



# Setup a custom theme for the plot ----------------------------------------

custom_theme <- theme(
  # axis values size
  axis.text = element_text(face = "plain",
                           size = 14,
                           color = "black"),
  # axis ticks lenght
  axis.ticks.length = unit(.15, "cm"),
  # axis ticks width
  axis.ticks = element_line(linewidth = 0.90,
                            color = "black"),
  # axis titles size
  axis.title = element_text(size = 16,
                            face = "plain",
                            color = "black"),
  axis.line = element_line(linewidth = 0.95,
                           color = "black"),
  legend.position = "none",
  panel.grid = element_blank(),
  panel.background = element_blank()
)

# ==========================================================================
# ==========================================================================





# ==========================================================================
# 3. Individual variance plots
# ==========================================================================


# Model 1 rank only --------------------------------------------------------
options(bitmapType = "cairo")
plot1_a <- ggplot(
  tab1_a,
  aes(x = Zcumul_xp,
      y = estimate__ / 4,
      color = predator_id)
) +
  geom_line(linewidth = 1, alpha = 0.5) +
  scale_color_viridis(discrete = TRUE, option = "D") + #B
  ylab("Hunting success\n") +
  ggtitle("Prey rank") +
  scale_y_continuous(breaks = seq(0, 1, 0.25),
                     limits = c(0, 1)) +
  scale_x_continuous(breaks = scaled_breaks,
                     labels = seq(0, 500, 100)) +
  xlab("\nCumulative experience") +
  custom_theme



# Model 2 rank + speed -----------------------------------------------------

plot2_a <- ggplot(
  tab2_a,
  aes(x = Zcumul_xp,
      y = estimate__ / 4,
      color = predator_id)
) +
  geom_line(linewidth = 1, alpha = 0.5) +
  scale_color_viridis(discrete = TRUE, option = "D") + #B
  ylab("Hunting success\n") +
  ggtitle("Prey rank + prey speed") +
  scale_y_continuous(breaks = seq(0, 1, 0.25),
                     limits = c(0, 1)) +
  scale_x_continuous(breaks = scaled_breaks,
                     labels = seq(0, 500, 100)) +
  xlab("\nCumulative experience") +
  custom_theme

# ==========================================================================
# ==========================================================================





# ==========================================================================
# 4. Global trend plots
# ==========================================================================


# Model 1 rank only --------------------------------------------------------

plot1_b <- ggplot(
  tab1_b,
  aes(x = Zcumul_xp,
      y = estimate__ / 4)
) +
  geom_vline(
    xintercept = min(tab1_b$Zcumul_xp),
    lty = "dashed",
    color = "#440154"
  ) +
  geom_text(
    aes(label = paste("y =", round(min(tab1_b$estimate__ / 4), digits = 2)),
        y = min(tab1_b$estimate__ / 4),
        x = min(tab1_b$Zcumul_xp) + 0.6),
    color = "#440154",
    size = 5
  ) +
  geom_vline(
    xintercept = tab1_b[which.max(estimate__), Zcumul_xp],
    lty = "dashed",
    color = "#440154"
  ) +
  geom_text(
    aes(label = paste("y =", format(round(max(tab1_b$estimate__ / 4), digits = 2), nsmall = 2)),
        y = max(tab1_b$estimate__ / 4) + 0.10,
        x = tab1_b[which.max(estimate__), Zcumul_xp] + 0.55),
    color = "#440154",
    size = 5
  ) +
  geom_ribbon(
    aes(
      ymin = lower__ / 4,
      ymax = upper__ / 4
    ),
    fill = "gray"
  ) +
  geom_line(linewidth = 1) +
  ylab("Hunting success\n") +
  ggtitle("Prey rank") +
  scale_y_continuous(
    breaks = seq(0, 1, 0.25),
    limits = c(0, 1)
  ) +
  scale_x_continuous(
    breaks = scaled_breaks,
    labels = seq(0, 500, 100)
  ) +
  xlab("\nCumulative experience") +
  custom_theme



# Model 2 rank + speed -----------------------------------------------------

# Predicted values from the GAMM
predicted_values <- tab2_b$estimate__

# Define the x-axis range
x <- tab2_b$Zcumul_xp

# Fit a spline to the predicted values
spline_fit <- splinefun(x, predicted_values)

# Calculate the first derivative using finite differences
dx <- mean(diff(x))
derivatives <- diff(predicted_values) / dx

# Establish the treshold at a value very close to 0 since its optimization
threshold <- 0.0088
# Find the point where the slope is close to zero
stabilized_point <- x[which(abs(derivatives) <= threshold)][1]

plot2_b <- ggplot(
  tab2_b,
  aes(x = Zcumul_xp,
      y = estimate__ / 4)
) +
  geom_vline(
    xintercept = min(tab2_b$Zcumul_xp),
    lty = "dashed",
    color = "#440154"
  ) +
  geom_text(
    aes(label = paste("y =", round(min(tab2_b$estimate__ / 4), digits = 2)),
        y = min(tab2_b$estimate__ / 4),
        x = min(tab2_b$Zcumul_xp) + 0.6),
    color = "#440154",
    size = 5
  ) +
  geom_vline(
    xintercept = tab2_b[Zcumul_xp == stabilized_point]$Zcumul_xp,
    lty = "dashed",
    color = "#440154"
  ) +
  geom_text(
    aes(label = paste(
      "y =",
      format(
        round(tab2_b[Zcumul_xp == stabilized_point]$estimate__ / 4, digits = 2),
        nsmall = 2
      )
    ),
    y = (tab2_b[Zcumul_xp == stabilized_point]$estimate__ / 4) + 0.10,
    x = tab2_b[Zcumul_xp == stabilized_point]$Zcumul_xp + 0.3),
    color = "#440154",
    size = 5
  ) +
  geom_ribbon(
    aes(
      ymin = lower__ / 4,
      ymax = upper__ / 4
    ),
    fill = "gray"
  ) +
  geom_line(linewidth = 1) +
  ylab("Hunting success\n") +
  ggtitle("Prey rank + prey speed") +
  scale_y_continuous(
    breaks = seq(0, 1, 0.25),
    limits = c(0, 1)
  ) +
  scale_x_continuous(
    breaks = scaled_breaks,
    labels = seq(0, 500, 100)
  ) +
  xlab("\nCumulative experience") +
  custom_theme

# ==========================================================================
# ==========================================================================





# ==========================================================================
# 3. Combine plots into 1 figure
# ==========================================================================

# Prepare figure ------------------------------------------------------------

# Arrange paneled figure
figure <- ggarrange(
  NULL, plot1_b, NULL, plot1_a,
  NULL, plot2_b, NULL, plot2_a,
  ncol = 4, nrow = 2,
  labels = c(
    "(A)", "", "(B)", "",
    "(C)", "", "(D)", ""
  ),
  widths = c(
    0.15, 1.5, 0.15, 1.5,
    0.15, 1.5, 0.15, 1.5
  )
)



# Export the figure -----------------------------------------------------

path <- file.path(getwd(), "outputs", "outputs_figures")

ggexport(
  figure,
  filename = file.path(path, "figure1.png"),
  width = 2700,
  height = 2200,
  res = 300
)

# ==========================================================================
# ==========================================================================

Thank you very much for your help!

I’m hopeful that the community will help me find answers.

Best

Maxime
session_infoLinux.txt (3.3 KB)
session_infoWindows.txt (2.8 KB)

[paul.buerkner]

this is a known issue that has been fixed in newer versions of brms. please install the latest cran or github version of brms and refit your model. this should fix the issue you see.

[Maxime]

Thank you very much for the swift response.

I’ll refit the models and get back on the post once I reach the solution.

Best

[Maxime]

Hello @paul.buerkner ,

Sorry for the late response. It seems that fitting the model on brms version 2.22.0 (on a computer cluster) and then computing the predictions on version 2.22.0 on my personal computer still results in different outputs on Pop OS compared to Windows.

The results make more sense biologically now (predictions computed under Linux), but the effect of my covariate is now negligible on the predictions.

This my code :

# ==========================================================================

#                           Plot GAMM models

# ==========================================================================



# ==========================================================================
# 1. Prepare script
# ==========================================================================



# Load libraries and model -------------------------------------------------

options(mc.cores = parallel::detectCores())

library(parallel)
library(brms)
library(rstan)
library(data.table)
library(ggplot2)
library(ggpubr)
library(viridis)

path <- file.path(getwd(), "outputs", "outputs_models")
path <- file.path(getwd(), "outputs", "outputs_models2024")

mod1 <- readRDS(file.path(path, "GAMM-II.rds"))
mod2 <- readRDS(file.path(path, "GAMM-V.rds"))



# Load the data ------------------------------------------------------------

data <- fread(
  "./data/FraserFrancoetal2023-data.csv",
  select = c("predator_id",
             "game_duration",
             "pred_speed",
             "prey_avg_speed",
             "prey_avg_rank",
             "cumul_xp_pred",
             "total_xp_pred",
             "hunting_success")
)

data[, predator_id := as.factor(predator_id)]



# Post processing preparations for custom family ---------------------------

expose_functions(mod1, vectorize = TRUE)

# Define the log likelihood function
log_lik_beta_binomial2 <- function(i, prep) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  phi <- brms::get_dpar(prep, "phi", i = i)
  trials <- prep$data$vint1[i]
  y <- prep$data$Y[i]
  beta_binomial2_lpmf(y, mu, phi, trials)
}

# Define function for posterior_predict
posterior_predict_beta_binomial2 <- function(i, prep, ...) {
  mu <- brms::get_dpar(prep, "mu", i = i)
  phi <- brms::get_dpar(prep, "phi", i = i)
  trials <- prep$data$vint1[i]
  beta_binomial2_rng(mu, phi, trials)
}

# Define function for posterior_epred
posterior_epred_beta_binomial2 <- function(prep) {
  mu <- brms::get_dpar(prep, "mu")
  trials <- prep$data$vint1
  trials <- matrix(trials, nrow = nrow(mu), ncol = ncol(mu), byrow = TRUE)
  mu * trials
}

# ==========================================================================
# ==========================================================================





# ==========================================================================
# 2. Prepare the data for the plots
# ==========================================================================



# Prepare the plots --------------------------------------------------------

# Group-level smooths model 1
tab1_a <- conditional_effects(
  mod1, method = "fitted",
  effects = "Zcumul_xp:predator_id",
  robust = TRUE, re_formula = NULL
)
tab1_a <- data.table(tab1_a[[1]])

# Cumulative XP global trend model 1
tab1_b <- conditional_effects(
  mod1, method = "fitted",
  robust = TRUE, re_formula = NULL,
  effects = "Zcumul_xp",
  conditions = data.frame(predator_id = NA)
)
tab1_b <- data.table(tab1_b[[1]])

# Group-level smooths model 2
tab2_a <- conditional_effects(
  mod2, method = "fitted",
  effects = "Zcumul_xp:predator_id",
  robust = TRUE, re_formula = NULL
)
tab2_a <- data.table(tab2_a[[1]])

# Cumulative XP global trend model 2
tab2_b <- conditional_effects(
  mod2, method = "fitted",
  robust = TRUE, re_formula = NULL,
  effects = "Zcumul_xp",
  conditions = data.frame(predator_id = NA)
)
tab2_b <- data.table(tab2_b[[1]])



# Transform values --------------------------------------------------------

# Back transform x-axis values
sequence <- (seq(0, 500, 100) - mean(data$cumul_xp_pred))
standev <- sd(data$cumul_xp_pred)
scaled_breaks <- sequence / standev

# List the tables
tables <- list(
  tab1_a, tab1_b,
  tab2_a, tab2_b
)
names(tables) <- c(
  "tab1_a", "tab1_b",
  "tab2_a", "tab2_b"
)

# Function to apply transformation
# Computes non standardized cumulative XP
func <- function(x) {
  x[, cumul_xp := (Zcumul_xp * standev) + mean(data$cumul_xp_pred)]
}

# Apply function
lapply(tables, func)



# Cut fitted values based on player XP ------------------------------------

# Extract player IDs with their total XP from the original data
xp <- unique(data[, .(predator_id, total_xp_pred)])

# Merge the two tables adding the total XP
tab1_a <- merge(tab1_a, xp, by = "predator_id")
tab2_a <- merge(tab2_a, xp, by = "predator_id")

# Cut all matches where fitted values are above total XP
tab1_a <- tab1_a[cumul_xp <= total_xp_pred, ]
tab2_a <- tab2_a[cumul_xp <= total_xp_pred, ]



# Setup a custom theme for the plot ----------------------------------------

custom_theme <- theme(
  # axis values size
  axis.text = element_text(face = "plain",
                           size = 14,
                           color = "black"),
  # axis ticks lenght
  axis.ticks.length = unit(.15, "cm"),
  # axis ticks width
  axis.ticks = element_line(linewidth = 0.90,
                            color = "black"),
  # axis titles size
  axis.title = element_text(size = 16,
                            face = "plain",
                            color = "black"),
  axis.line = element_line(linewidth = 0.95,
                           color = "black"),
  legend.position = "none",
  panel.grid = element_blank(),
  panel.background = element_blank()
)

# ==========================================================================
# ==========================================================================





# ==========================================================================
# 3. Individual variance plots
# ==========================================================================


# Model 1 rank only --------------------------------------------------------

plot1_a <- ggplot(
  tab1_a,
  aes(x = Zcumul_xp,
      y = estimate__ / 4,
      color = predator_id)
) +
  #geom_line(linewidth = 1, alpha = 0.5) +
  geom_line(linewidth = 1) +
  scale_color_viridis(discrete = TRUE, option = "D") + #B
  ylab("Hunting success\n") +
  ggtitle("Prey rank") +
  scale_y_continuous(breaks = seq(0, 1, 0.25),
                     limits = c(0, 1)) +
  scale_x_continuous(breaks = scaled_breaks,
                     labels = seq(0, 500, 100)) +
  xlab("\nCumulative experience") +
  custom_theme



# Model 2 rank + speed -----------------------------------------------------

plot2_a <- ggplot(
  tab2_a,
  aes(x = Zcumul_xp,
      y = estimate__ / 4,
      color = predator_id)
) +
  #geom_line(linewidth = 1, alpha = 0.5) +
  geom_line(linewidth = 1) +
  scale_color_viridis(discrete = TRUE, option = "D") + #B
  ylab("Hunting success\n") +
  ggtitle("Prey rank + prey speed") +
  scale_y_continuous(breaks = seq(0, 1, 0.25),
                     limits = c(0, 1)) +
  scale_x_continuous(breaks = scaled_breaks,
                     labels = seq(0, 500, 100)) +
  xlab("\nCumulative experience") +
  custom_theme

# ==========================================================================
# ==========================================================================





# ==========================================================================
# 4. Global trend plots
# ==========================================================================


# Model 1 rank only --------------------------------------------------------

plot1_b <- ggplot(
  tab1_b,
  aes(x = Zcumul_xp,
      y = estimate__ / 4)
) +
  geom_vline(
    xintercept = min(tab1_b$Zcumul_xp),
    lty = "dashed",
    color = "#440154"
  ) +
  geom_text(
    aes(label = paste("y =", round(min(tab1_b$estimate__ / 4), digits = 2)),
        y = min(tab1_b$estimate__ / 4),
        x = min(tab1_b$Zcumul_xp) + 0.6),
    color = "#440154",
    size = 5
  ) +
  geom_vline(
    xintercept = tab1_b[which.max(estimate__), Zcumul_xp],
    lty = "dashed",
    color = "#440154"
  ) +
  geom_text(
    aes(label = paste("y =", format(round(max(tab1_b$estimate__ / 4), digits = 2), nsmall = 2)),
        y = max(tab1_b$estimate__ / 4) + 0.10,
        x = tab1_b[which.max(estimate__), Zcumul_xp] + 0.55),
    color = "#440154",
    size = 5
  ) +
  geom_ribbon(
    aes(
      ymin = lower__ / 4,
      ymax = upper__ / 4
    ),
    fill = "gray"
  ) +
  geom_line(linewidth = 1) +
  ylab("Hunting success\n") +
  ggtitle("Prey rank") +
  scale_y_continuous(
    breaks = seq(0, 1, 0.25),
    limits = c(0, 1)
  ) +
  scale_x_continuous(
    breaks = scaled_breaks,
    labels = seq(0, 500, 100)
  ) +
  xlab("\nCumulative experience") +
  custom_theme



# Model 2 rank + speed -----------------------------------------------------

# Predicted values from the GAMM
predicted_values <- tab2_b$estimate__

# Define the x-axis range
x <- tab2_b$Zcumul_xp

# Fit a spline to the predicted values
spline_fit <- splinefun(x, predicted_values)

# Calculate the first derivative using finite differences
dx <- mean(diff(x))
derivatives <- diff(predicted_values) / dx

# Establish the treshold at a value very close to 0 since its optimization
threshold <- 0.0088
# Find the point where the slope is close to zero
stabilized_point <- x[which(abs(derivatives) <= threshold)][1]

plot2_b <- ggplot(
  tab2_b,
  aes(x = Zcumul_xp,
      y = estimate__ / 4)
) +
  geom_vline(
    xintercept = min(tab2_b$Zcumul_xp),
    lty = "dashed",
    color = "#440154"
  ) +
  geom_text(
    aes(label = paste("y =", round(min(tab2_b$estimate__ / 4), digits = 2)),
        y = min(tab2_b$estimate__ / 4),
        x = min(tab2_b$Zcumul_xp) + 0.6),
    color = "#440154",
    size = 5
  ) +
  geom_vline(
    xintercept = tab2_b[Zcumul_xp == stabilized_point]$Zcumul_xp,
    lty = "dashed",
    color = "#440154"
  ) +
  geom_text(
    aes(label = paste(
      "y =",
      format(
        round(tab2_b[Zcumul_xp == stabilized_point]$estimate__ / 4, digits = 2),
        nsmall = 2
      )
    ),
    y = (tab2_b[Zcumul_xp == stabilized_point]$estimate__ / 4) + 0.10,
    x = tab2_b[Zcumul_xp == stabilized_point]$Zcumul_xp + 0.3),
    color = "#440154",
    size = 5
  ) +
  geom_ribbon(
    aes(
      ymin = lower__ / 4,
      ymax = upper__ / 4
    ),
    fill = "gray"
  ) +
  geom_line(linewidth = 1) +
  ylab("Hunting success\n") +
  ggtitle("Prey rank + prey speed") +
  scale_y_continuous(
    breaks = seq(0, 1, 0.25),
    limits = c(0, 1)
  ) +
  scale_x_continuous(
    breaks = scaled_breaks,
    labels = seq(0, 500, 100)
  ) +
  xlab("\nCumulative experience") +
  custom_theme

# ==========================================================================
# ==========================================================================





# ==========================================================================
# 3. Combine plots into 1 figure
# ==========================================================================

# Prepare figure ------------------------------------------------------------

# Arrange paneled figure
figure <- ggarrange(
  NULL, plot1_b, NULL, plot1_a,
  NULL, plot2_b, NULL, plot2_a,
  ncol = 4, nrow = 2,
  labels = c(
    "(A)", "", "(B)", "",
    "(C)", "", "(D)", ""
  ),
  widths = c(
    0.15, 1.5, 0.15, 1.5,
    0.15, 1.5, 0.15, 1.5
  )
)



# Export the figure -----------------------------------------------------

path <- file.path(getwd(), "outputs", "outputs_figures")

ggexport(
  figure,
  filename = file.path(path, "figure1New.png"),
  width = 2700,
  height = 2200,
  res = 300
)

# ==========================================================================
# ==========================================================================

Another thing I noted is that expose_functions() now returned an error which did not occur before :

Error in expose_stan_functions(self$functions, global, verbose) : 
  'stanmodel' is not a valid object

At this point I’m a bit worried that I cannot trust the results of neither analyses that I did. What do you suggest I do?

[Maxime]

When I try to run the predictions under Windows, the results also changed and are not consistend with those under Linux.

1. As we see on Linux, the results show no effect of controlling for the prey’s speed.
2. Compared to the results on Linux, the trends are completely different.

Under Windows, expose_functions (which used to work) now return:

expose_functions(modA2, vectorize = TRUE)
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
 error in 'model7d0294e3545_User_defined_functions' at line 10, column 2
  -------------------------------------------------
     8:    *   an integer sequence from start to end
     9:    */
    10:   array[] int sequence(int start, int end) {
         ^
    11:     array[end - start + 1] int seq;
  -------------------------------------------------

PARSER EXPECTED: "}"
Error in stanc(model_code = mc, model_name = "User-defined functions",  :
  failed to parse Stan model 'User-defined functions' due to the above error.

[Maxime]

So one last thing.

I forgot to update brms on Windows before computing the predictions.

1. Results are now the same as those under Linux
2. I get the same error message for expose_functions()

Is it required that I used this to compute the predictions using conditional_effects() since I am using the beta_binomial family or not?

Thank you!

Maxime

[frank_hezemans]

The specific error message refers to this helper function, which is typically pasted at the top of the Stan code generated by brms (in the functions block).

Just a thought: Perhaps your version of rstan / Stan is out of date? The array[] int instead of int[] syntax is required as of Stan 2.32+ (August 2023).

Otherwise brms::expose_functions is basically just a wrapper around rstan::expose_stan_functions (or the expose_functions method for cmdstanr).

And regarding your last question: yes, as far as I understand, expose_functions is required here, because conditional_effects needs to generate draws from the (expected value of) the posterior predictive distribution of your custom model.

[paul.buerkner]

Hi,

(1) glad this is resolved!

(2) Since you have updated to the latest version of brms now, you can use the built-in beta-binomial family of brms via family = beta_binomial(). This means that you can remove all the custom Stan code etc. that you copied from the brms_customfamilies vignette, including the expose_functions part.

Paul

[Maxime]

Thanks a lot @frank_hezemans and @paul.buerkner for your input.

I have one last question which may not really be answerable.

Do you have an hypothesis as to why:

• on brms 2.17, results indicate a clear effect of prey speed (under Windows)
• but this effect is no longer there on brms 2.22 (under Windows and Linux)?

By that, I mean that the predictions (as well as variance parameters) changed when I controlled for prey speed before (first figure on thread), but they do not anymore (second figure on thread). This changes the whole interpretation of my results.

Best,

Maxime

[paul.buerkner]

There was a bug that could potentially lead to wrong spline predictions if the predictions were made under a different OS than where the model was fitted (see Spline model prediction inconsistencies when changing BLAS/LAPACK implementation · Issue #1465 · paul-buerkner/brms · GitHub). This bug was quite hard to find and we only fixed it relatively recently in brms 2.20.

[Maxime]

I see!

Many thanks for the help. Much appreciated.

Best,

Maxime