基于贝叶斯框架的Wilcoxon检验需求：配对前后小样本比例数据分析

Great question—switching to Bayesian methods here is a smart move given your small, non-normal paired proportional data. Let’s break down how to approach this step by step, with practical, actionable code and interpretation tips.

Why This Is a Better Fit Than NHST #

• Your data is paired proportions (0-1) with n=10: Frequentist methods like t-tests rely on asymptotic normality, which doesn’t hold here, and Wilcoxon tests only focus on ranks (not the magnitude of proportional differences).
• Bayesian methods let you quantify uncertainty directly via posterior distributions, rather than just a single p-value. You’ll also get intuitive metrics like the probability that post-intervention scores are higher than pre-intervention scores.

Step 1: Model Choice for Paired Proportions #

Since your data is bounded between 0 and 1, the difference between paired observations (d_i = post_i - pre_i) ranges from -1 to 1. A robust, easy-to-implement approach is to use a scaled Beta distribution (since Beta is designed for [0,1] data, we can shift/scale it to fit the [-1,1] difference range). Alternatively, a Bayesian analog to the Wilcoxon test works if you care more about rank differences than magnitude.

Step 2: Practical Implementation with R (brms Package) #

We’ll use brms because it’s user-friendly and integrates well with tidy workflows.

1. Prepare Your Data #

First, calculate the paired difference and scale it to the [0,1] range (required for the Beta distribution):

# Assume your data has columns: subject, pre, post
data$diff <- data$post - data$pre
data$diff_scaled <- (data$diff + 1) / 2  # Shifts [-1,1] to [0,1]

2. Fit the Bayesian Model #

We’ll use a non-informative prior for the intercept (safe when you don’t have prior knowledge):

library(brms)

model <- brm(
  formula = diff_scaled ~ 1,
  data = data,
  family = Beta(link = "logit"),
  prior = prior(normal(0, 1), class = Intercept),  # Non-informative prior
  sample_prior = "yes",
  chains = 4,  # Standard for reliable posterior sampling
  iter = 2000,
  warmup = 1000  # Discard first half of samples as burn-in
)

3. Interpret the Results #

• Summary Output: Run summary(model) to get the posterior mean, standard deviation, and 95% credible interval (CI) for the intercept. To convert this back to your original difference scale:
  1. Apply the inverse logit to the intercept to get the mean of diff_scaled.
  2. Multiply by 2 and subtract 1 to shift back to [-1,1].
• Posterior Distribution: Use bayesplot::mcmc_dens(posterior_samples(model)) to visualize the posterior of the scaled difference, then transform it to see the actual difference distribution.
• Probability of Effect: Calculate how many posterior samples correspond to a positive difference (post > pre):

  post_samples <- posterior_samples(model)
  post_samples$diff_original <- (plogis(post_samples$b_Intercept) * 2) - 1
  prob_positive <- mean(post_samples$diff_original > 0)
  cat("Probability post > pre:", prob_positive)
  

  This is a far more intuitive metric than a p-value—it tells you directly how likely it is that the intervention had a positive effect.

Step 3: Bayesian Analog to Wilcoxon Signed-Rank Test #

If you want to focus on rank differences (like the Wilcoxon test), use the bayesfactor package to compute a Bayesian signed-rank test:

library(bayesfactor)

bf <- bayesfactor_signed_rank(data$post, data$pre)
print(bf)

The Bayes Factor (BF) quantifies evidence for the alternative hypothesis (post ≠ pre) vs the null. A BF >3 is moderate evidence, >10 is strong evidence for a difference.

Key Takeaways #

• Uncertainty: You’ll see exactly how much variation there is in your effect size estimate, not just a binary "significant" or "not".
• Small Sample Robustness: Bayesian methods don’t require large sample sizes to produce reliable results.
• Flexibility: If you have prior knowledge about the effect size (e.g., from similar studies), you can update the priors to incorporate that information.

内容的提问来源于stack exchange，提问作者Adrian Santos