# Disrepancies in output of conditional\_effects() between Windows and Linux for the same model

**URL:** <https://discourse.mc-stan.org/t/disrepancies-in-output-of-conditional-effects-between-windows-and-linux-for-the-same-model/36526>\
**Category:** brms\
**Tags:** conditional-effects, brms\
**Created:** [September 8, 2024, 2:54pm UTC](https://discourse.mc-stan.org/t/disrepancies-in-output-of-conditional-effects-between-windows-and-linux-for-the-same-model/36526 "2024-09-08T14:54:16Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Maxime](https://avatars.discourse-cdn.com/v4/letter/m/a9adbd/32.png) [@Maxime](https://discourse.mc-stan.org/u/Maxime)\
**Post date:** [September 8, 2024, 2:54pm UTC](https://discourse.mc-stan.org/t/disrepancies-in-output-of-conditional-effects-between-windows-and-linux-for-the-same-model/36526/1 "2024-09-08T14:54:16Z")

</div>

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 :

 ![Screenshot from 2024-09-08 10-37-57](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/7/7/77237803a2f7c42ba8642f23a84dbfe9342fa8f8.jpeg)

- 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](https://osf.io/hdv38/)

Here is the link to the raw data : [GitHub - quantitative-ecologist/predator-expertise: Data and code for Fraser Franco et al. XXXX](https://github.com/quantitative-ecologist/predator-expertise)

Here is the code to generate the plots.

```no-highlight
# ==========================================================================

# 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](https://discourse.mc-stan.org/uploads/short-url/A8d7PCFIrODiQchdCdyNROpMfDu.txt) (3.3 KB)  
[session\_infoWindows.txt](https://discourse.mc-stan.org/uploads/short-url/bT48EPZ3aFRRdxTWSDZ2bk2DkHt.txt) (2.8 KB)

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [September 8, 2024, 5:07pm UTC](https://discourse.mc-stan.org/t/disrepancies-in-output-of-conditional-effects-between-windows-and-linux-for-the-same-model/36526/2 "2024-09-08T17:07:57Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Maxime](https://avatars.discourse-cdn.com/v4/letter/m/a9adbd/32.png) [@Maxime](https://discourse.mc-stan.org/u/Maxime)\
**Post date:** [September 8, 2024, 5:25pm UTC](https://discourse.mc-stan.org/t/disrepancies-in-output-of-conditional-effects-between-windows-and-linux-for-the-same-model/36526/3 "2024-09-08T17:25:32Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Maxime](https://avatars.discourse-cdn.com/v4/letter/m/a9adbd/32.png) [@Maxime](https://discourse.mc-stan.org/u/Maxime)\
**Post date:** [February 9, 2025, 9:47pm UTC](https://discourse.mc-stan.org/t/disrepancies-in-output-of-conditional-effects-between-windows-and-linux-for-the-same-model/36526/4 "2025-02-09T21:47:40Z")

</div>

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 :

```no-highlight
# ==========================================================================

# 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 :

```no-highlight
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?

 ![figure1New](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/c/0/c0587f800e7726745b1366e8164d477e7ddce914.png)

---

<div class="post-metadata">

**Author:** ![Maxime](https://avatars.discourse-cdn.com/v4/letter/m/a9adbd/32.png) [@Maxime](https://discourse.mc-stan.org/u/Maxime)\
**Post date:** [February 9, 2025, 10:48pm UTC](https://discourse.mc-stan.org/t/disrepancies-in-output-of-conditional-effects-between-windows-and-linux-for-the-same-model/36526/5 "2025-02-09T22:48:41Z")

</div>

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:

```no-highlight
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.

```

 ![figure1New](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/b/1/b13595e98e1f6f8f60e6780b7227e6a5c31f887b.jpeg)

---

<div class="post-metadata">

**Author:** ![Maxime](https://avatars.discourse-cdn.com/v4/letter/m/a9adbd/32.png) [@Maxime](https://discourse.mc-stan.org/u/Maxime)\
**Post date:** [February 10, 2025, 12:15am UTC](https://discourse.mc-stan.org/t/disrepancies-in-output-of-conditional-effects-between-windows-and-linux-for-the-same-model/36526/6 "2025-02-10T00:15:32Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![frank\_hezemans](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/frank_hezemans/32/17981_2.png) [@frank\_hezemans](https://discourse.mc-stan.org/u/frank_hezemans)\
**Post date:** [February 10, 2025, 9:13am UTC](https://discourse.mc-stan.org/t/disrepancies-in-output-of-conditional-effects-between-windows-and-linux-for-the-same-model/36526/7 "2025-02-10T09:13:40Z")

</div>

The specific error message refers to [this helper function](https://github.com/paul-buerkner/brms/blob/master/inst/chunks/fun_sequence.stan), 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.

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [February 10, 2025, 9:56am UTC](https://discourse.mc-stan.org/t/disrepancies-in-output-of-conditional-effects-between-windows-and-linux-for-the-same-model/36526/8 "2025-02-10T09:56:32Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Maxime](https://avatars.discourse-cdn.com/v4/letter/m/a9adbd/32.png) [@Maxime](https://discourse.mc-stan.org/u/Maxime)\
**Post date:** [February 10, 2025, 10:28pm UTC](https://discourse.mc-stan.org/t/disrepancies-in-output-of-conditional-effects-between-windows-and-linux-for-the-same-model/36526/9 "2025-02-10T22:28:55Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [February 11, 2025, 6:09am UTC](https://discourse.mc-stan.org/t/disrepancies-in-output-of-conditional-effects-between-windows-and-linux-for-the-same-model/36526/10 "2025-02-11T06:09:36Z")

</div>

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](https://github.com/paul-buerkner/brms/issues/1465)). This bug was quite hard to find and we only fixed it relatively recently in brms 2.20.

---

<div class="post-metadata">

**Author:** ![Maxime](https://avatars.discourse-cdn.com/v4/letter/m/a9adbd/32.png) [@Maxime](https://discourse.mc-stan.org/u/Maxime)\
**Post date:** [February 11, 2025, 4:46pm UTC](https://discourse.mc-stan.org/t/disrepancies-in-output-of-conditional-effects-between-windows-and-linux-for-the-same-model/36526/11 "2025-02-11T16:46:17Z")

</div>

I see!

Many thanks for the help. Much appreciated.

Best,

Maxime