Hello everyone,

I am currently fitting a Gamma GLM on some artificial data. The response variable is the biomass values of fish, and the explanatory variables are Enforcement (an ordered factor), coral cover (continuous, ranging from 0 to 1), and habitat (a factor with two levels).

I am encountering an issue where the estimates from the model indicate the opposite direction of the conditional estimates. The predicted values match my expectations (as shown in the figure), where higher enforcement levels correspond to higher biomass values, but the estimates are negative.

```
> summary(reef_glm)
Family: gamma
Links: mu = log; shape = identity
Formula: biomass_fish ~ coral_cover:habitat + enforcement + habitat
Data: reef (Number of observations: 387)
Draws: 3 chains, each with iter = 10000; warmup = 4000; thin = 1;
total post-warmup draws = 18000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.33 0.26 0.81 1.84 1.00 12547 12273
enforcement.L 0.83 0.07 0.70 0.96 1.00 15739 12768
enforcement.Q 0.20 0.07 0.06 0.34 1.00 14039 13048
enforcement.C -0.44 0.07 -0.57 -0.31 1.00 14452 12850
enforcementE4 -0.22 0.08 -0.37 -0.07 1.00 15555 11882
habitatShallow -1.41 0.37 -2.13 -0.70 1.00 9978 10821
habitatDeep:coral_cover -0.96 0.38 -1.69 -0.22 1.00 12544 12160
habitatShallow:coral_cover -0.28 0.42 -1.10 0.56 1.00 12482 12522
```

My question is twofold:

- Why are the conditional estimates accurate, but the estimates are showing the opposite direction?
- Should I report the estimates in my research paper, or is it sufficient to report the conditional estimates?

Thank you for your help!

Here is the data:

artificial_reef.csv (14.3 KB)

And below is the code:

```
# GLM for reef data
# packages ----------------------------------------------------------------
library(tidyverse) # data wrangling; ggplot
library(brms) # brm() and the rest of Bayesian functions
# data --------------------------------------------------------------------
reef <- read_csv("artificial_reef.csv") |>
mutate(enforcement = factor(enforcement,
levels = c("No", "Low", "Medium_South",
"Medium_North", "High"),
ordered = TRUE)) |>
glimpse()
# exploration -------------------------------------------------------------
# boxplot
reef |>
ggplot(aes(enforcement, biomass_fish)) +
geom_boxplot(outlier.colour = NA) +
ylim(NA, 6) +
theme_minimal()
# data summary
reef |>
group_by(enforcement) |>
summarise(means = mean(biomass_fish),
medians = median(biomass_fish)) |>
# do other Enforcement levels have lower biomass
mutate(means_dif = means - first(means),
medians_dif = medians - first(medians)) |>
select(enforcement, means, means_dif, medians, medians_dif)
# model -------------------------------------------------------------------
# GLM (linear model)
reef_glm <- brm(biomass_fish ~ coral_cover:habitat + enforcement + habitat,
family = Gamma(link = "log"),
data = reef,
cores = 3,
chains = 3, iter = 10000, warmup = 4000,
file = "brms_reef_glm.RDS")
summary(reef_glm)
conditional_effects(reef_glm, "enforcement")
```

- R version 4.2.2 (2022-10-31 ucrt)
- Platform: x86_64-w64-mingw32/x64 (64-bit)
- Running under: Windows 10 x64 (build 22621)
- brms Version: brms_2.18.0