如何在贝叶斯线性混合效应回归中分析截距与斜率差异？

Context & Model Setup #

I’m analyzing Arctic greenup dates using a Bayesian linear mixed-effects model (via blme) to address two core questions:

1. How do meteorological factors affect greenup dates within each of the 16 pan-Arctic regions?
2. Are there significant differences in baseline greenup dates (intercepts) across regions, and if so, which regions differ?

Scaled Data Preparation #

library(readr)
library(blme)
data <- read_csv("data.csv")
data.sc <- data.frame(scale(data))
head(data.sc)

Output:

regions      year    greenup          V1         V2          V3        V4          V5       V6            V7         V8       V9        V10
1          1 -1.687957  0.9336382  0.85187965 0.71761508  3.16360255 1.6670476 -0.68640856 -0.655334953 -1.799065  1.5577978 3.0125643 -0.5009276
2          1 -1.601395  0.2370278 -0.30385233 0.59926804  0.04203437 1.3281044 -0.62391808 -0.592095987 -1.590912  0.6735797 1.2016789 -0.5955858
3          2 -1.514833  0.8030237 -0.05341494 0.05770376 -0.03230812 0.4116001 -0.34697632 -0.311837094 -1.428521  0.2459790 0.7906408 -0.7097304
4          2 -1.428272 -0.1548155 -0.52504749 0.21499928 -0.53588809 1.3991222 -0.70552873 -0.674684140 -1.491065 -0.5336937 0.1615725 -0.6405863
5          3 -1.341710  0.4111804 -0.55526813 0.41495915 -0.28220429 1.0430792 -0.69965688 -0.668741960 -1.440828 -0.1372223 0.3667727 -0.8735882
6          3 -1.255148  1.7717476  0.37503107 0.49582150  3.25563398 1.1052932 -0.06417199 -0.005461234 -1.859388  1.6666640 3.6155112 -0.6626455

Bayesian Mixed Model (to avoid singular boundary warnings) #

reg <- blmer(greenup ~ V1+V2+V3+V4+V5+V6+V7+V8+(1+V1+V2+V3+V4+V5+V6+V7+V8|regions), data = data.sc, REML = F)
ranef(reg)

Output:

(Intercept)           V1          V2          V3          V4           V5          V6        V7          V8
1   -1.5088397 -0.678780774 -0.28176494  0.03297057  0.14762373  0.135628050  0.03870133 -0.16881483  0.01049209
2    1.4277155  0.002705333  0.63878199  0.04527960 -0.29932021 -0.092542944 -0.05871342  0.12192872  0.16659908
3   -1.1971171 -0.643673842 -0.36725321 -0.00437876  0.09999908 -0.195088268  0.02960333 -0.07320157  0.21971584
4   -1.8725315 -0.785756089 -0.38010681 -0.02575994  0.28474568  0.840094435 -0.13247112 -0.05102267 -0.77963082
5    1.0132713 -0.429509626 -0.31425681  0.29430628 -0.19690372  0.282548814  0.39303350 -0.12090616  0.69670371
6   -2.4881956 -1.692458471 -0.09449177  0.01172166  0.06284188 -0.844244673  0.79118303 -0.34123995  1.04264827
7    3.1058386  0.826832058  1.11178582  0.02043532 -0.47408031 -0.338428604 -0.33133031  0.31261663  0.10124955
8    1.9322901  3.317059011 -0.30085324 -0.22707802  0.28574031 -0.569026041 -0.42553220  0.17865665 -0.47420030
9   -1.2204600 -2.438705399  1.87732120 -0.07449666 -0.38793703  0.188712028 -0.38356492  0.11440041 -0.52707496
10   2.0227413  1.395666163  0.67593119 -0.14830203 -0.16997007 -0.385349601 -0.34510500  0.26727338 -0.24816029
11  -1.9223249 -0.767426078 -0.23096931 -0.10032745  0.14049704  0.003901228  0.01478041 -0.07175487 -0.12908571
12   0.0149983 -0.492942263 -0.40391390  0.16250771 -0.06522357  0.566016764  0.35380727 -0.11917404  0.19111680
13  -0.7919515 -0.731003237 -0.34630359  0.17115953  0.10735851  0.940255206 -0.07730108 -0.11185924 -0.37309296
14   4.7216459  4.388771043 -1.08519982  0.01923765  0.10673146 -0.205195431 -0.03391026  0.21661697  0.12766651
15   0.1163438 -0.623755530 -0.09248300  0.15112399 -0.08205915  0.193081476  0.03723803 -0.03449855  0.24479648
16  -3.3534245 -0.647022299 -0.40622380 -0.32839946  0.43995636 -0.520362438  0.12958139 -0.11902088 -0.26974330

1. Analyzing Random Intercept (Greenup Date) Differences Between Regions #

Step 1: Test if Regional Intercept Variation is Significant #

Compare your full model to a fixed-effects-only model using Bayes factors to assess if regional intercept variation is meaningful:

# Fixed-effects-only model
reg_null <- blmer(greenup ~ V1+V2+V3+V4+V5+V6+V7+V8, data = data.sc, REML = F)

# Calculate Bayes Factor
library(bayestestR)
bayes_factor(reg, reg_null)

A Bayes factor >3 indicates strong evidence in favor of the full model (i.e., significant regional intercept variation).

Step 2: Extract Full Regional Intercepts #

Each region’s intercept is the sum of the global fixed intercept and its random intercept deviation:

# Get fixed intercept
fixed_intercept <- fixef(reg)[1]

# Get random intercepts
random_intercepts <- ranef(reg)$regions[, "(Intercept)"]

# Calculate full intercept per region
regional_intercepts <- fixed_intercept + random_intercepts
regional_intercepts <- data.frame(region = names(regional_intercepts), intercept = regional_intercepts)

Step 3: Post-Hoc Pairwise Comparisons #

Use posterior credible intervals to identify regions with significantly different intercepts:

# Extract posterior samples
posterior <- as.data.frame(reg)

# Create pairwise difference matrix
n_regions <- 16
diff_matrix <- matrix(NA, nrow = n_regions, ncol = n_regions)
rownames(diff_matrix) <- colnames(diff_matrix) <- 1:n_regions

for(i in 1:n_regions) {
  for(j in 1:n_regions) {
    if(i != j) {
      # Compute posterior samples of intercept difference
      diff_samples <- posterior[, "(Intercept)"] + posterior[, paste0("regions.(Intercept)", i)] - 
                      (posterior[, "(Intercept)"] + posterior[, paste0("regions.(Intercept)", j)])
      # Check if 95% CI excludes 0
      ci <- quantile(diff_samples, c(0.025, 0.975))
      diff_matrix[i,j] <- ifelse(ci[1] > 0 | ci[2] < 0, "Significant", "Not Significant")
    }
  }
}

print(diff_matrix)

Visualize Regional Intercepts #

library(ggplot2)
ggplot(regional_intercepts, aes(x = region, y = intercept)) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = intercept - 1.96*sqrt(vcov(reg)[paste0("regions.(Intercept)", region), paste0("regions.(Intercept)", region)]),
                    ymax = intercept + 1.96*sqrt(vcov(reg)[paste0("regions.(Intercept)", region), paste0("regions.(Intercept)", region)])),
                width = 0.2) +
  labs(x = "Region", y = "Scaled Greenup Date Intercept", title = "Regional Baseline Greenup Dates (95% CIs)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

2. Exploring Region-Specific Effects of Meteorological Variables #

Step 1: Calculate Region-Specific Slopes #

Each region’s slope for a variable is the sum of the global fixed slope and its random slope deviation:

# Get fixed slopes (exclude intercept)
fixed_slopes <- fixef(reg)[-1]

# Get random slopes (exclude intercept)
random_slopes <- ranef(reg)$regions[, -1]

# Compute full slopes per region
regional_slopes <- t(apply(random_slopes, 1, function(x) fixed_slopes + x))
regional_slopes <- as.data.frame(regional_slopes)
regional_slopes$region <- rownames(regional_slopes)

Step 2: Visualize Region-Specific Effects #

library(tidyr)
# Reshape data for plotting
slopes_long <- pivot_longer(regional_slopes, cols = V1:V8, names_to = "Variable", values_to = "Slope")

ggplot(slopes_long, aes(x = region, y = Slope, color = Variable)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  facet_wrap(~Variable, scales = "free_y") +
  labs(x = "Region", y = "Scaled Slope", title = "Region-Specific Meteorological Effects") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Step 3: Test if Slope Variation is Significant #

Compare your full model to one with only random intercepts (no random slopes) to assess if region-specific slopes are meaningful:

# Model with only random intercepts
reg_no_rand_slopes <- blmer(greenup ~ V1+V2+V3+V4+V5+V6+V7+V8+(1|regions), data = data.sc, REML = F)

# Calculate Bayes Factor
bayes_factor(reg, reg_no_rand_slopes)

A Bayes factor >3 indicates strong evidence that slopes vary by region.

Step 4: Identify Variables with Significant Region-Specific Effects #

For individual variables, compare models with/