calculate_posterior <- function(prior_mean, prior_se,
data_mean, data_se,
bias_mean, bias_se) {
# Precisions for prior and data measurement
prior_var <- prior_se^2
prior_precision <- 1 / prior_var
data_var <- data_se^2
data_precision <- 1 / data_var # Precision of the measurement itself
# Precision related to bias belief
bias_var <- bias_se^2
# Handle case where bias_se is 0 (certainty about bias) to avoid division by zero
if (bias_var == 0) {
bias_precision <- Inf
} else {
bias_precision <- 1 / bias_var # How certain are we about the bias value
}
# Calculate data's effective precision for estimating tau
# Uncertainty about bias adds to the data's measurement variance
total_data_uncertainty_variance <- bias_var + data_var
if (total_data_uncertainty_variance == 0) { # Should only happen if bias_se=0 and data_se=0
effective_data_precision <- Inf
} else {
effective_data_precision <- 1 / total_data_uncertainty_variance
}
# Calculate posterior precision (certainty about tau)
posterior_precision <- prior_precision + effective_data_precision
posterior_precision <- prior_precision + effective_data_precision
# Calculate posterior variance and SE for tau
if (posterior_precision == Inf) { # Handle certainty case
posterior_var <- 0
} else {
posterior_var <- 1 / posterior_precision
}
posterior_se <- sqrt(posterior_var)
# Calculate posterior mean for tau
# Weighted average using prior precision and data's *effective* precision
# Note: Data point is corrected by expected bias (data_mean - bias_mean)
# Handle edge case of infinite precisions
if (is.infinite(prior_precision) && is.infinite(effective_data_precision)) {
# This case is ill-defined unless prior_mean equals (data_mean - bias_mean)
# For simplicity, let's prioritize data if both are infinitely precise.
posterior_mean <- data_mean - bias_mean
} else if (is.infinite(prior_precision)) {
posterior_mean <- prior_mean
} else if (is.infinite(effective_data_precision)) {
posterior_mean <- data_mean - bias_mean
} else if (posterior_precision == 0) { # Both prior and effective data precision are 0 (infinite variance)
posterior_mean <- NA # Undefined / Non-informative
} else {
posterior_mean <- (prior_precision * prior_mean +
effective_data_precision * (data_mean - bias_mean)) / posterior_precision
}
# Return results
list(prior_mean = prior_mean, prior_se = prior_se, prior_precision = prior_precision,
data_mean = data_mean, data_se = data_se, data_precision = data_precision,
bias_mean = bias_mean, bias_se = bias_se, bias_precision = bias_precision,
effective_data_precision = effective_data_precision,
posterior_mean = posterior_mean, posterior_se = posterior_se, posterior_precision = posterior_precision)
}
# --- Example cases ---
# Base info from Fig 11.1
prior_mean_base <- 0
prior_se_base <- 2
data_mean_base <- 10
data_se_base <- 1
# Case A: Certainty of no bias (Bias SE = 0)
results_A <- calculate_posterior(prior_mean_base, prior_se_base, data_mean_base, data_se_base, bias_mean = 0, bias_se = 0)
# Case B: Moderate uncertainty about bias (Bias SE = 2)
results_B <- calculate_posterior(prior_mean_base, prior_se_base, data_mean_base, data_se_base, bias_mean = 0, bias_se = 2)
# Case C: High uncertainty about bias (Bias SE = 10)
results_C <- calculate_posterior(prior_mean_base, prior_se_base, data_mean_base, data_se_base, bias_mean = 0, bias_se = 10)
# --- Print results ---
print_results <- function(results, case_label) {
cat("--- Case:", case_label, "---\n")
cat(sprintf("Bias Beliefs: Mean=%.1f, SE=%.1f (Bias Var=%.1f, Bias Precision=%.4f)\n",
results$bias_mean, results$bias_se, results$bias_se^2, results$bias_precision))
cat(sprintf("Data Effective Precision: %.4f\n", results$effective_data_precision))
cat(sprintf("Posterior: Mean=%.2f, SE=%.2f (Posterior Precision=%.4f)\n\n",
results$posterior_mean, results$posterior_se, results$posterior_precision))
}
print_results(results_A, "A: Certain No Bias")
