Publication Bias Models with MARS
mars authors
2026-07-17
Publication-Bias-Models.RmdThis vignette demonstrates publication-bias modeling with
publication_bias(method = "beta_selection") and
publication_bias(method = "beta_binomial") and
publication_bias(method = "logistic_selection") after
fitting mars models.
Current implementation note: - Univariate publication-bias fits
estimate a scalar tau2 inside the selection model. -
Multilevel and multivariate publication-bias fits use the fitted MARS
covariance blocks directly, rather than collapsing to diagonal sampling
variances. - Reported standard errors use a conservative rule: when both
Hessian-based and OPG-based SEs are finite, the larger value is
reported.
library(mars)
extract_full_data_coefficients <- function(fit) {
est_raw <- if(!is.null(fit$beta_r)) {
fit$beta_r
} else if(!is.null(fit$beta_fe)) {
fit$beta_fe
} else {
fit$est_values$mu
}
est_vec <- as.numeric(est_raw)
est_names <- names(est_raw)
if(is.matrix(est_raw) || is.data.frame(est_raw)) {
est_vec <- as.numeric(est_raw[, 1, drop = TRUE])
est_names <- rownames(est_raw)
}
if(length(est_names) != length(est_vec)) {
est_names <- colnames(fit$design_matrix)
}
se_vec <- rep(NA_real_, length(est_vec))
vc <- fit$varcov_beta
if(is.matrix(vc) && nrow(vc) == length(est_vec)) {
se_vec <- sqrt(pmax(diag(vc), 0))
if(!is.null(colnames(vc)) && length(colnames(vc)) == length(est_vec)) {
est_names <- colnames(vc)
}
}
data.frame(
term = est_names,
estimate = est_vec,
std_error = se_vec,
stringsAsFactors = FALSE
)
}
compare_pub_bias_3 <- function(x, y, z,
label_x = "beta_selection",
label_y = "beta_binomial",
label_z = "logistic_selection",
ref = NULL,
label_ref = "all_data") {
a <- x$coefficients
b <- y$coefficients
c <- z$coefficients
names(a)[names(a) == "estimate"] <- paste0("estimate_", label_x)
names(a)[names(a) == "std_error"] <- paste0("se_", label_x)
names(b)[names(b) == "estimate"] <- paste0("estimate_", label_y)
names(b)[names(b) == "std_error"] <- paste0("se_", label_y)
names(c)[names(c) == "estimate"] <- paste0("estimate_", label_z)
names(c)[names(c) == "std_error"] <- paste0("se_", label_z)
out <- merge(
a[, c("term", paste0("estimate_", label_x), paste0("se_", label_x))],
b[, c("term", paste0("estimate_", label_y), paste0("se_", label_y))],
by = "term",
all = TRUE
)
out <- merge(
out,
c[, c("term", paste0("estimate_", label_z), paste0("se_", label_z))],
by = "term",
all = TRUE
)
if(!is.null(ref)) {
r <- ref
names(r)[names(r) == "estimate"] <- paste0("estimate_", label_ref)
names(r)[names(r) == "std_error"] <- paste0("se_", label_ref)
out <- merge(
out,
r[, c("term", paste0("estimate_", label_ref), paste0("se_", label_ref))],
by = "term",
all = TRUE
)
}
out[order(out$term), ]
}
plot_pub_bias_compare_3 <- function(tbl,
label_x = "beta_selection",
label_y = "beta_binomial",
label_z = "logistic_selection",
label_ref = "all_data",
main = "Method Comparison") {
ref_est_name <- paste0("estimate_", label_ref)
ref_se_name <- paste0("se_", label_ref)
est_x <- tbl[[paste0("estimate_", label_x)]]
se_x <- tbl[[paste0("se_", label_x)]]
est_y <- tbl[[paste0("estimate_", label_y)]]
se_y <- tbl[[paste0("se_", label_y)]]
est_z <- tbl[[paste0("estimate_", label_z)]]
se_z <- tbl[[paste0("se_", label_z)]]
has_ref <- ref_est_name %in% names(tbl)
est_ref <- if(has_ref) tbl[[ref_est_name]] else rep(NA_real_, nrow(tbl))
se_ref <- if(has_ref && ref_se_name %in% names(tbl)) tbl[[ref_se_name]] else rep(NA_real_, nrow(tbl))
k <- nrow(tbl)
xpos <- seq_len(k)
offsets <- if(has_ref) c(-0.27, -0.09, 0.09, 0.27) else c(-0.18, 0, 0.18)
all_vals <- c(
est_x, est_y, est_z, est_ref,
est_x - 1.96 * se_x, est_x + 1.96 * se_x,
est_y - 1.96 * se_y, est_y + 1.96 * se_y,
est_z - 1.96 * se_z, est_z + 1.96 * se_z,
est_ref - 1.96 * se_ref, est_ref + 1.96 * se_ref
)
ylim <- range(all_vals[is.finite(all_vals)])
old_par <- par(no.readonly = TRUE)
on.exit(par(old_par), add = TRUE)
x_axis_cex <- max(0.55, min(0.8, 8 / max(nchar(tbl$term), 1L)))
mar <- old_par$mar
mar[1] <- max(mar[1], 8 + 3 * x_axis_cex)
par(mar = mar)
plot(xpos, est_x, type = "n", xaxt = "n", xlab = "", ylab = "Estimate", main = main, ylim = ylim)
abline(h = 0, lty = 2, col = "gray60")
axis(1, at = xpos, labels = tbl$term, las = 2, cex.axis = x_axis_cex)
title(xlab = "Coefficient Term", line = 6)
keep_x <- is.finite(est_x) & is.finite(se_x)
keep_y <- is.finite(est_y) & is.finite(se_y)
keep_z <- is.finite(est_z) & is.finite(se_z)
keep_ref <- has_ref & is.finite(est_ref) & is.finite(se_ref)
pt_x <- is.finite(est_x)
pt_y <- is.finite(est_y)
pt_z <- is.finite(est_z)
pt_ref <- has_ref & is.finite(est_ref)
miss_se_x <- is.finite(est_x) & !is.finite(se_x)
miss_se_y <- is.finite(est_y) & !is.finite(se_y)
miss_se_z <- is.finite(est_z) & !is.finite(se_z)
miss_se_ref <- has_ref & is.finite(est_ref) & !is.finite(se_ref)
x_x <- xpos + offsets[1]
x_y <- xpos + offsets[2]
x_z <- xpos + offsets[3]
x_ref <- if(has_ref) xpos + offsets[4] else xpos
if(any(keep_x)) {
segments(x_x[keep_x], est_x[keep_x] - 1.96 * se_x[keep_x],
x_x[keep_x], est_x[keep_x] + 1.96 * se_x[keep_x], col = "#1F78B4", lwd = 2)
}
if(any(keep_y)) {
segments(x_y[keep_y], est_y[keep_y] - 1.96 * se_y[keep_y],
x_y[keep_y], est_y[keep_y] + 1.96 * se_y[keep_y], col = "#E31A1C", lwd = 2)
}
if(any(keep_z)) {
segments(x_z[keep_z], est_z[keep_z] - 1.96 * se_z[keep_z],
x_z[keep_z], est_z[keep_z] + 1.96 * se_z[keep_z], col = "#33A02C", lwd = 2)
}
if(any(keep_ref)) {
segments(x_ref[keep_ref], est_ref[keep_ref] - 1.96 * se_ref[keep_ref],
x_ref[keep_ref], est_ref[keep_ref] + 1.96 * se_ref[keep_ref], col = "#4D4D4D", lwd = 2)
}
if(any(pt_x)) {
points(x_x[pt_x], est_x[pt_x], pch = 19, col = "#1F78B4")
}
if(any(pt_y)) {
points(x_y[pt_y], est_y[pt_y], pch = 17, col = "#E31A1C")
}
if(any(pt_z)) {
points(x_z[pt_z], est_z[pt_z], pch = 15, col = "#33A02C")
}
if(any(pt_ref)) {
points(x_ref[pt_ref], est_ref[pt_ref], pch = 18, col = "#4D4D4D")
}
if(any(miss_se_x)) {
points(x_x[miss_se_x], est_x[miss_se_x], pch = 1, col = "#1F78B4", cex = 1.2, lwd = 1.5)
}
if(any(miss_se_y)) {
points(x_y[miss_se_y], est_y[miss_se_y], pch = 2, col = "#E31A1C", cex = 1.2, lwd = 1.5)
}
if(any(miss_se_z)) {
points(x_z[miss_se_z], est_z[miss_se_z], pch = 0, col = "#33A02C", cex = 1.2, lwd = 1.5)
}
if(any(miss_se_ref)) {
points(x_ref[miss_se_ref], est_ref[miss_se_ref], pch = 5, col = "#4D4D4D", cex = 1.2, lwd = 1.5)
}
legend_labels <- c(label_x, label_y, label_z)
legend_cols <- c("#1F78B4", "#E31A1C", "#33A02C")
legend_pch <- c(19, 17, 15)
if(has_ref) {
legend_labels <- c(legend_labels, label_ref)
legend_cols <- c(legend_cols, "#4D4D4D")
legend_pch <- c(legend_pch, 18)
}
legend("topright", legend = legend_labels, col = legend_cols, pch = legend_pch, bty = "n")
if(any(miss_se_x) || any(miss_se_y) || any(miss_se_z) || any(miss_se_ref)) {
mtext("Open symbols indicate estimates with unavailable SE.", side = 3, line = 0.2, cex = 0.8)
}
}
extract_alt_interval <- function(x, term = "(Intercept)") {
if(!is.null(x$summary) &&
all(c("point_estimate", "ci_lower", "ci_upper") %in% colnames(x$summary)) &&
term %in% rownames(x$summary)) {
out <- x$summary[term, c("point_estimate", "ci_lower", "ci_upper"), drop = FALSE]
return(out)
}
if(!is.null(x$summary) &&
all(c("estimate", "ci_lower", "ci_upper") %in% colnames(x$summary)) &&
term %in% rownames(x$summary)) {
out <- x$summary[term, c("estimate", "ci_lower", "ci_upper"), drop = FALSE]
names(out) <- c("point_estimate", "ci_lower", "ci_upper")
return(out)
}
stop("Could not extract a point estimate + interval for term: ", term)
}Univariate Example
fit_uni <- mars(
formula = yi ~ 1,
data = bcg,
studyID = "trial",
variance = "vi",
varcov_type = "univariate",
structure = "univariate",
estimation_method = "MLE"
)
summary(fit_uni)
#> Results generated with MARS:v 0.5.3
#> Friday, July 17, 2026
#>
#> Model Type:
#> univariate
#>
#> Estimation Method:
#> Maximum Likelihood
#>
#> Model Formula:
#> yi ~ 1
#>
#> Data Summary:
#> Number of Effect Sizes: 13
#> Number of Fixed Effects: 1
#> Number of Random Effects: 1
#>
#> Random Components:
#> term var SD
#> intercept 0.28 0.5292
#>
#> Fixed Effects Estimates:
#> attribute estimate SE z_test p_value lower upper
#> (Intercept) -0.7112 0.1719 -4.137 3.513e-05 -1.048 -0.3743
#>
#> Model Fit Statistics:
#> logLik Dev AIC BIC AICc
#> -12.67 31.33 29.33 30.46 37.9
#>
#> Q Error: 152.233 (12), p < 0.0001
#>
#> I2 (General): 92.03957
#>
#> Residual Diagnostics:
#> n n_finite_raw mean_raw sd_raw rmse mae q_pearson mean_abs_studentized
#> 13 13 -0.02945 0.6938 0.6672 0.5958 13.16 0.9527
#> max_abs_studentized prop_abs_studentized_gt2 prop_abs_studentized_gt3
#> 1.434 0 0
#>
#> Normality (whitened residuals): test n_tested statistic p_value
#> shapiro_wilk_whitened 13 0.886 0.08592
#>
#> Heteroscedasticity trend (|raw residual| ~ fitted): n corr_abs_raw_fitted slope p_value
#> 13 NA NA NA
pb_uni <- publication_bias(
fit_uni,
method = "beta_selection",
p_value_tail = "two.sided",
maxit = 100,
integration_points = 31
)
summary(pb_uni)
#> Publication Bias Model (MARS)
#> Method: beta_selection
#> Structure: univariate
#> Effects: 13 | Studies: 13
#> P-value tail: two.sided
#> Convergence: 0 - Converged
#> Objective: 11.71
#>
#> Selection Parameters:
#> parameter estimate
#> alpha 1.0857
#> beta 0.5901
#> tau2 0.3713
#>
#> Selection Weight Profile (p-value grid):
#> model p_value weight
#> beta_selection 0.0100 0.6767
#> beta_selection 0.1189 0.8776
#> beta_selection 0.2278 0.9794
#> beta_selection 0.3367 1.0778
#> beta_selection 0.4456 1.1882
#> beta_selection 0.5544 1.3242
#> beta_selection 0.6633 1.5084
#> beta_selection 0.7722 1.7936
#> beta_selection 0.8811 2.3679
#> beta_selection 0.9900 6.5975
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> (Intercept) -0.9144 0.32 -2.8579 0.0043
#>
#> Adjusted vs Unadjusted Summary:
#> term unadjusted_estimate unadjusted_se unadjusted_ci_lower
#> (Intercept) -0.7112 0.1719 -1.0481
#> unadjusted_ci_upper estimate std_error adjusted_ci_lower adjusted_ci_upper
#> -0.3743 -0.9144 0.32 -1.5416 -0.2873
#> estimate_diff percent_change direction_changed ci_overlap
#> -0.2032 -28.5761 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1S, naive-t)
#> Intercept: -7.8519 SE: 3.2091 t: -2.4468 df: 12 p: 0.0308
#>
#> SE note:
#> SE rule: reported SEs use the larger of the Hessian-based and OPG-based values when both are finite, otherwise the finite one.
pb_uni_bb <- publication_bias(
fit_uni,
method = "beta_binomial",
n_bins = 8,
p_value_tail = "two.sided",
maxit = 100,
integration_points = 31
)
summary(pb_uni_bb)
#> Publication Bias Model (MARS)
#> Method: beta_binomial
#> Structure: univariate
#> Effects: 13 | Studies: 13
#> P-value tail: two.sided
#> Convergence: 0 - Converged
#> Objective: -29.08
#>
#> Selection Parameters:
#> parameter estimate
#> alpha 2753.5116
#> beta 0.2679
#> tau2 2.4061
#> n_bins 8.0000
#>
#> Selection Weight Profile (bin probabilities):
#> bin p_lower p_upper weight
#> 1 0.000 0.125 0.0000
#> 2 0.125 0.250 0.0000
#> 3 0.250 0.375 0.0000
#> 4 0.375 0.500 0.0000
#> 5 0.500 0.625 0.0000
#> 6 0.625 0.750 0.0000
#> 7 0.750 0.875 0.0007
#> 8 0.875 1.000 0.9993
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> (Intercept) -0.4604 1.1726 -0.3926 0.6946
#>
#> Adjusted vs Unadjusted Summary:
#> term unadjusted_estimate unadjusted_se unadjusted_ci_lower
#> (Intercept) -0.7112 0.1719 -1.0481
#> unadjusted_ci_upper estimate std_error adjusted_ci_lower adjusted_ci_upper
#> -0.3743 -0.4604 1.1726 -2.7587 1.8379
#> estimate_diff percent_change direction_changed ci_overlap
#> 0.2508 35.2681 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1S, naive-t)
#> Intercept: -7.8519 SE: 3.2091 t: -2.4468 df: 12 p: 0.0308
#>
#> SE note:
#> SE rule: reported SEs use the larger of the Hessian-based and OPG-based values when both are finite, otherwise the finite one. beta_binomial uses discretized p-value bins; coefficient magnitudes can differ materially from beta_selection.
pb_uni_logit <- publication_bias(
fit_uni,
method = "logistic_selection",
p_value_tail = "two.sided",
maxit = 100,
integration_points = 31
)
summary(pb_uni_logit)
#> Publication Bias Model (MARS)
#> Method: logistic_selection
#> Structure: univariate
#> Effects: 13 | Studies: 13
#> P-value tail: two.sided
#> Convergence: 0 - Converged
#> Objective: 12.67
#>
#> Selection Parameters:
#> parameter estimate
#> logit_intercept 7.6548
#> logit_slope 0.3333
#> tau2 0.2801
#>
#> Selection Weight Profile (p-value grid):
#> model p_value weight
#> logistic_selection 0.0100 0.9999
#> logistic_selection 0.1189 0.9998
#> logistic_selection 0.2278 0.9997
#> logistic_selection 0.3367 0.9997
#> logistic_selection 0.4456 0.9996
#> logistic_selection 0.5544 0.9996
#> logistic_selection 0.6633 0.9996
#> logistic_selection 0.7722 0.9996
#> logistic_selection 0.8811 0.9995
#> logistic_selection 0.9900 0.9995
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> (Intercept) -0.7111 0.5224 -1.3613 0.1734
#>
#> Adjusted vs Unadjusted Summary:
#> term unadjusted_estimate unadjusted_se unadjusted_ci_lower
#> (Intercept) -0.7112 0.1719 -1.0481
#> unadjusted_ci_upper estimate std_error adjusted_ci_lower adjusted_ci_upper
#> -0.3743 -0.7111 0.5224 -1.7351 0.3128
#> estimate_diff percent_change direction_changed ci_overlap
#> 1e-04 0.0079 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1S, naive-t)
#> Intercept: -7.8519 SE: 3.2091 t: -2.4468 df: 12 p: 0.0308
#>
#> SE note:
#> SE rule: reported SEs use the larger of the Hessian-based and OPG-based values when both are finite, otherwise the finite one.Selection and adjusted fixed effects:
pb_uni$selection
#> alpha beta tau2
#> 1.0856923 0.5901224 0.3713397
pb_uni$covariance_mode
#> [1] "publication_bias_tau2"
pb_uni$coefficients
#> term estimate std_error z_value p_value
#> 1 (Intercept) -0.914432 0.3199691 -2.857876 0.004264869
pb_uni_bb$selection_weights
#> bin p_lower p_upper weight
#> 1 1 0.000 0.125 2.931151e-22
#> 2 2 0.125 0.250 9.013700e-19
#> 3 3 0.250 0.375 1.413949e-15
#> 4 4 0.375 0.500 1.521505e-12
#> 5 5 0.500 0.625 1.283419e-09
#> 6 6 0.625 0.750 9.363074e-07
#> 7 7 0.750 0.875 6.790423e-04
#> 8 8 0.875 1.000 9.993200e-01
pb_uni_logit$selection
#> logit_intercept logit_slope tau2
#> 7.6548399 0.3333277 0.2800653
pb_uni_logit$coefficients
#> term estimate std_error z_value p_value
#> 1 (Intercept) -0.7111428 0.5224129 -1.361266 0.1734297For the univariate case, the publication-bias likelihood still
estimates a scalar tau2, so the selection parameter vector
includes tau2 along with the selection-function
parameters.
ref_uni <- extract_full_data_coefficients(fit_uni)
cmp_uni <- compare_pub_bias_3(pb_uni, pb_uni_bb, pb_uni_logit, ref = ref_uni)
cmp_uni
#> term estimate_beta_selection se_beta_selection estimate_beta_binomial
#> 1 (Intercept) -0.914432 0.3199691 -0.460373
#> se_beta_binomial estimate_logistic_selection se_logistic_selection
#> 1 1.172612 -0.7111428 0.5224129
#> estimate_all_data se_all_data
#> 1 -0.7111991 0.1718968
plot_pub_bias_compare_3(cmp_uni, main = "Univariate: publication-bias methods with all-data reference")
Funnel plot from the publication-bias fit:
funnel_plot(pb_uni)
Jackknife version of the complete-data fit and each publication-bias model:
set.seed(2026)
# Keep the vignette jackknife examples intentionally small/lightweight.
bcg_jack <- bcg[seq_len(min(6, nrow(bcg))), , drop = FALSE]
jack_complete <- mars_alt_estimation(
type = "jackknife",
data = bcg_jack,
formula = yi ~ 1,
studyID = "trial",
effectID = NULL,
sample_size = NULL,
variance = "vi",
varcov_type = "univariate",
structure = "univariate",
estimation_method = "MLE",
jackknife_target = "complete_data"
)
jack_beta_selection <- mars_alt_estimation(
type = "jackknife",
data = bcg_jack,
formula = yi ~ 1,
studyID = "trial",
effectID = NULL,
sample_size = NULL,
variance = "vi",
varcov_type = "univariate",
structure = "univariate",
estimation_method = "MLE",
jackknife_target = "beta_selection",
publication_bias_args = list(
p_value_tail = "two.sided",
maxit = 100,
integration_points = 31
)
)
jack_beta_binomial <- mars_alt_estimation(
type = "jackknife",
data = bcg_jack,
formula = yi ~ 1,
studyID = "trial",
effectID = NULL,
sample_size = NULL,
variance = "vi",
varcov_type = "univariate",
structure = "univariate",
estimation_method = "MLE",
jackknife_target = "beta_binomial",
publication_bias_args = list(
p_value_tail = "two.sided",
n_bins = 8,
maxit = 100,
integration_points = 31
)
)
jack_logistic <- mars_alt_estimation(
type = "jackknife",
data = bcg_jack,
formula = yi ~ 1,
studyID = "trial",
effectID = NULL,
sample_size = NULL,
variance = "vi",
varcov_type = "univariate",
structure = "univariate",
estimation_method = "MLE",
jackknife_target = "logistic_selection",
publication_bias_args = list(
p_value_tail = "two.sided",
maxit = 100,
integration_points = 31
)
)
jack_complete_stats <- compute_alt_stats(jack_complete, type = "jackknife")
jack_beta_selection_stats <- compute_alt_stats(jack_beta_selection, type = "jackknife")
jack_beta_binomial_stats <- compute_alt_stats(jack_beta_binomial, type = "jackknife")
jack_logistic_stats <- compute_alt_stats(jack_logistic, type = "jackknife")
jackknife_compare <- rbind(
complete_data = extract_alt_interval(jack_complete_stats),
beta_selection = extract_alt_interval(jack_beta_selection_stats),
beta_binomial = extract_alt_interval(jack_beta_binomial_stats),
logistic_selection = extract_alt_interval(jack_logistic_stats)
)
jackknife_compareThe univariate jackknife illustration uses only the first six trials so the vignette demonstrates the workflow without refitting all leave-one-out models from the full dataset.
jack_est <- jackknife_compare[, "point_estimate"]
jack_lo <- jackknife_compare[, "ci_lower"]
jack_hi <- jackknife_compare[, "ci_upper"]
xpos <- seq_len(nrow(jackknife_compare))
plot(
xpos, jack_est,
ylim = range(c(jack_lo, jack_hi)),
xaxt = "n", xlab = "", ylab = "Jackknife Estimate",
pch = 19, col = "#1F78B4",
main = "Univariate Jackknife: complete data and publication-bias targets"
)
segments(xpos, jack_lo, xpos, jack_hi, lwd = 2, col = "#1F78B4")
abline(h = 0, lty = 2, col = "gray60")
axis(1, at = xpos, labels = rownames(jackknife_compare), las = 2)
title(xlab = "Jackknife Target", line = 6)Multilevel Example
fit_ml <- mars(
formula = effect ~ 1 + (1 | district/study),
data = school,
studyID = "district",
variance = "var",
varcov_type = "multilevel",
structure = "multilevel",
estimation_method = "MLE"
)
summary(fit_ml)
#> Results generated with MARS:v 0.5.3
#> Friday, July 17, 2026
#>
#> Model Type:
#> multilevel
#>
#> Estimation Method:
#> Maximum Likelihood
#>
#> Model Formula:
#> effect ~ 1 + (1 | district/study)
#>
#> Data Summary:
#> Number of Effect Sizes: 56
#> Number of Fixed Effects: 1
#> Number of Random Effects: 2
#>
#> Random Data Summary:
#> Number unique district: 11
#> Number unique study: 56
#>
#> Random Components:
#> term var SD
#> district 0.05774 0.2403
#> study 0.03286 0.1813
#>
#> Fixed Effects Estimates:
#> attribute estimate SE z_test p_value lower upper
#> (Intercept) 0.1845 0.08048 2.292 0.02191 0.02671 0.3422
#>
#> Model Fit Statistics:
#> logLik Dev AIC BIC AICc
#> -8.395 24.79 22.79 28.87 24.5
#>
#> Q Error: 578.864 (55), p < 0.0001
#>
#> I2 (General):
#> names values
#> district 92.38
#> study 87.34
#>
#> I2 (Jackson): 98.69621
#> I2 (Between): 86.3727
#>
#> Residual Diagnostics:
#> n n_finite_raw mean_raw sd_raw rmse mae q_pearson mean_abs_studentized
#> 56 56 -0.06446 0.3258 0.3293 0.2632 96.8 1.065
#> max_abs_studentized prop_abs_studentized_gt2 prop_abs_studentized_gt3
#> 4.196 0.1071 0.03571
#>
#> Normality (whitened residuals): test n_tested statistic p_value
#> shapiro_wilk_whitened 56 0.9555 0.03766
#>
#> Heteroscedasticity trend (|raw residual| ~ fitted): n corr_abs_raw_fitted slope p_value
#> 56 NA NA NA
pb_ml <- publication_bias(
fit_ml,
method = "beta_selection",
p_value_tail = "two.sided",
maxit = 100,
integration_points = 31
)
summary(pb_ml)
#> Publication Bias Model (MARS)
#> Method: beta_selection
#> Structure: multilevel
#> Effects: 56 | Studies: 11
#> P-value tail: two.sided
#> Convergence: 0 - Converged
#> Objective: -142.5
#>
#> Selection Parameters:
#> parameter estimate
#> alpha 1757.9268
#> beta 0.0007
#>
#> Selection Weight Profile (p-value grid):
#> model p_value weight
#> beta_selection 0.0100 0
#> beta_selection 0.1189 0
#> beta_selection 0.2278 0
#> beta_selection 0.3367 0
#> beta_selection 0.4456 0
#> beta_selection 0.5544 0
#> beta_selection 0.6633 0
#> beta_selection 0.7722 0
#> beta_selection 0.8811 0
#> beta_selection 0.9900 0
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> (Intercept) 0.9458 0.0361 26.2164 0
#>
#> Adjusted vs Unadjusted Summary:
#> term unadjusted_estimate unadjusted_se unadjusted_ci_lower
#> (Intercept) 0.1845 0.0805 0.0267
#> unadjusted_ci_upper estimate std_error adjusted_ci_lower adjusted_ci_upper
#> 0.3422 0.9458 0.0361 0.8751 1.0165
#> estimate_diff percent_change direction_changed ci_overlap
#> 0.7613 412.7472 FALSE FALSE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1S, naive-t)
#> Intercept: 3.69 SE: 1.8503 t: 1.9943 df: 10 p: 0.0741
#>
#> SE note:
#> SE rule: reported SEs use the larger of the Hessian-based and OPG-based values when both are finite, otherwise the finite one. Covariance structure: used the fitted MARS within/between-study covariance blocks directly for the publication-bias likelihood.
pb_ml_bb <- publication_bias(
fit_ml,
method = "beta_binomial",
n_bins = 8,
p_value_tail = "two.sided",
maxit = 100,
integration_points = 31
)
summary(pb_ml_bb)
#> Publication Bias Model (MARS)
#> Method: beta_binomial
#> Structure: multilevel
#> Effects: 56 | Studies: 11
#> P-value tail: two.sided
#> Convergence: 0 - Converged
#> Objective: -400.3
#>
#> Selection Parameters:
#> parameter estimate
#> alpha 158.6091
#> beta 0.8091
#> n_bins 8.0000
#>
#> Selection Weight Profile (bin probabilities):
#> bin p_lower p_upper weight
#> 1 0.000 0.125 0.0000
#> 2 0.125 0.250 0.0000
#> 3 0.250 0.375 0.0000
#> 4 0.375 0.500 0.0000
#> 5 0.500 0.625 0.0000
#> 6 0.625 0.750 0.0011
#> 7 0.750 0.875 0.0332
#> 8 0.875 1.000 0.9656
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> (Intercept) 1.2761 0.034 37.5387 0
#>
#> Adjusted vs Unadjusted Summary:
#> term unadjusted_estimate unadjusted_se unadjusted_ci_lower
#> (Intercept) 0.1845 0.0805 0.0267
#> unadjusted_ci_upper estimate std_error adjusted_ci_lower adjusted_ci_upper
#> 0.3422 1.2761 0.034 1.2094 1.3427
#> estimate_diff percent_change direction_changed ci_overlap
#> 1.0916 591.7976 FALSE FALSE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1S, naive-t)
#> Intercept: 3.69 SE: 1.8503 t: 1.9943 df: 10 p: 0.0741
#>
#> SE note:
#> SE rule: reported SEs use the larger of the Hessian-based and OPG-based values when both are finite, otherwise the finite one. Covariance structure: used the fitted MARS within/between-study covariance blocks directly for the publication-bias likelihood. beta_binomial uses discretized p-value bins; coefficient magnitudes can differ materially from beta_selection.
pb_ml_logit <- publication_bias(
fit_ml,
method = "logistic_selection",
p_value_tail = "two.sided",
maxit = 100,
integration_points = 31
)
summary(pb_ml_logit)
#> Publication Bias Model (MARS)
#> Method: logistic_selection
#> Structure: multilevel
#> Effects: 56 | Studies: 11
#> P-value tail: two.sided
#> Convergence: 0 - Converged
#> Objective: 8.396
#>
#> Selection Parameters:
#> parameter estimate
#> logit_intercept 7.8501
#> logit_slope 11.7805
#>
#> Selection Weight Profile (p-value grid):
#> model p_value weight
#> logistic_selection 0.0100 1.0000
#> logistic_selection 0.1189 1.0000
#> logistic_selection 0.2278 1.0000
#> logistic_selection 0.3367 1.0000
#> logistic_selection 0.4456 1.0000
#> logistic_selection 0.5544 1.0000
#> logistic_selection 0.6633 1.0000
#> logistic_selection 0.7722 1.0000
#> logistic_selection 0.8811 0.9999
#> logistic_selection 0.9900 0.9997
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> (Intercept) 0.1844 0.087 2.1183 0.0341
#>
#> Adjusted vs Unadjusted Summary:
#> term unadjusted_estimate unadjusted_se unadjusted_ci_lower
#> (Intercept) 0.1845 0.0805 0.0267
#> unadjusted_ci_upper estimate std_error adjusted_ci_lower adjusted_ci_upper
#> 0.3422 0.1844 0.087 0.0138 0.355
#> estimate_diff percent_change direction_changed ci_overlap
#> -1e-04 -0.0299 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1S, naive-t)
#> Intercept: 3.69 SE: 1.8503 t: 1.9943 df: 10 p: 0.0741
#>
#> SE note:
#> SE rule: reported SEs use the larger of the Hessian-based and OPG-based values when both are finite, otherwise the finite one. Covariance structure: used the fitted MARS within/between-study covariance blocks directly for the publication-bias likelihood.Multilevel coefficients with standard errors:
pb_ml$covariance_mode
#> [1] "fixed_mars_structure"
pb_ml$covariance_structure
#> $structure
#> [1] "multilevel"
#>
#> $Tau
#> [,1] [,2]
#> [1,] 0.05773805 0.00000000
#> [2,] 0.00000000 0.03286473
pb_ml$coefficients
#> term estimate std_error z_value p_value
#> 1 (Intercept) 0.9457898 0.03607629 26.21638 1.728599e-151
pb_ml_bb$coefficients
#> term estimate std_error z_value p_value
#> 1 (Intercept) 1.276058 0.03399314 37.53869 2.154797e-308
pb_ml_logit$coefficients
#> term estimate std_error z_value p_value
#> 1 (Intercept) 0.1844002 0.08704963 2.118334 0.03414679For multilevel fits, publication_bias() now uses the
fitted MARS block covariance structure directly. In that case the output
stores covariance_mode = "fixed_mars_structure" and the
publication-bias selection parameters no longer include a separate
scalar tau2.
ref_ml <- extract_full_data_coefficients(fit_ml)
cmp_ml <- compare_pub_bias_3(pb_ml, pb_ml_bb, pb_ml_logit, ref = ref_ml)
cmp_ml
#> term estimate_beta_selection se_beta_selection estimate_beta_binomial
#> 1 (Intercept) 0.9457898 0.03607629 1.276058
#> se_beta_binomial estimate_logistic_selection se_logistic_selection
#> 1 0.03399314 0.1844002 0.08704963
#> estimate_all_data se_all_data
#> 1 0.1844554 0.0804815
plot_pub_bias_compare_3(cmp_ml, main = "Multilevel: publication-bias methods with all-data reference")
Jackknife version of the multilevel complete-data fit and each publication-bias model:
set.seed(2026)
school_jack <- subset(school, district %in% sort(unique(district))[1:4])
jack_ml_complete <- mars_alt_estimation(
type = "jackknife",
data = school_jack,
formula = effect ~ 1 + (1 | district/study),
studyID = "district",
effectID = NULL,
sample_size = NULL,
variance = "var",
varcov_type = "multilevel",
structure = "multilevel",
estimation_method = "MLE",
jackknife_target = "complete_data"
)
jack_ml_beta_selection <- mars_alt_estimation(
type = "jackknife",
data = school_jack,
formula = effect ~ 1 + (1 | district/study),
studyID = "district",
effectID = NULL,
sample_size = NULL,
variance = "var",
varcov_type = "multilevel",
structure = "multilevel",
estimation_method = "MLE",
jackknife_target = "beta_selection",
publication_bias_args = list(
p_value_tail = "two.sided",
maxit = 80,
integration_points = 31
)
)
jack_ml_beta_binomial <- mars_alt_estimation(
type = "jackknife",
data = school_jack,
formula = effect ~ 1 + (1 | district/study),
studyID = "district",
effectID = NULL,
sample_size = NULL,
variance = "var",
varcov_type = "multilevel",
structure = "multilevel",
estimation_method = "MLE",
jackknife_target = "beta_binomial",
publication_bias_args = list(
p_value_tail = "two.sided",
n_bins = 8,
maxit = 80,
integration_points = 31
)
)
jack_ml_logistic <- mars_alt_estimation(
type = "jackknife",
data = school_jack,
formula = effect ~ 1 + (1 | district/study),
studyID = "district",
effectID = NULL,
sample_size = NULL,
variance = "var",
varcov_type = "multilevel",
structure = "multilevel",
estimation_method = "MLE",
jackknife_target = "logistic_selection",
publication_bias_args = list(
p_value_tail = "two.sided",
maxit = 80,
integration_points = 31
)
)
jack_ml_complete_stats <- compute_alt_stats(jack_ml_complete, type = "jackknife")
jack_ml_beta_selection_stats <- compute_alt_stats(jack_ml_beta_selection, type = "jackknife")
jack_ml_beta_binomial_stats <- compute_alt_stats(jack_ml_beta_binomial, type = "jackknife")
jack_ml_logistic_stats <- compute_alt_stats(jack_ml_logistic, type = "jackknife")
jackknife_compare_ml <- rbind(
complete_data = extract_alt_interval(jack_ml_complete_stats),
beta_selection = extract_alt_interval(jack_ml_beta_selection_stats),
beta_binomial = extract_alt_interval(jack_ml_beta_binomial_stats),
logistic_selection = extract_alt_interval(jack_ml_logistic_stats)
)
jackknife_compare_mlThe multilevel jackknife illustration uses only the first four top-level clusters so the vignette stays quick to run while still showing the leave-one- cluster-out workflow.
jack_ml_est <- jackknife_compare_ml[, "point_estimate"]
jack_ml_lo <- jackknife_compare_ml[, "ci_lower"]
jack_ml_hi <- jackknife_compare_ml[, "ci_upper"]
xpos_ml <- seq_len(nrow(jackknife_compare_ml))
plot(
xpos_ml, jack_ml_est,
ylim = range(c(jack_ml_lo, jack_ml_hi)),
xaxt = "n", xlab = "", ylab = "Jackknife Estimate",
pch = 19, col = "#1F78B4",
main = "Multilevel Jackknife: complete data and publication-bias targets"
)
segments(xpos_ml, jack_ml_lo, xpos_ml, jack_ml_hi, lwd = 2, col = "#1F78B4")
abline(h = 0, lty = 2, col = "gray60")
axis(1, at = xpos_ml, labels = rownames(jackknife_compare_ml), las = 2)
title(xlab = "Jackknife Target", line = 6)Multivariate Example
becker09 <- na.omit(becker09)
fit_mv <- mars(
data = becker09,
studyID = "ID",
effectID = "numID",
sample_size = "N",
effectsize_type = "cor",
varcov_type = "weighted",
variable_names = c(
"Cognitive_Performance",
"Somatic_Performance",
"Selfconfidence_Performance",
"Somatic_Cognitive",
"Selfconfidence_Cognitive",
"Selfconfidence_Somatic"
)
)
summary(fit_mv)
#> Results generated with MARS:v 0.5.3
#> Friday, July 17, 2026
#>
#> Model Type:
#> multivariate
#>
#> Estimation Method:
#> Restricted Maximum Likelihood
#>
#> Model Formula:
#> NULL
#>
#> Data Summary:
#> Number of Effect Sizes: 48
#> Number of Fixed Effects: 6
#> Number of Random Effects: 6
#>
#> Random Components:
#> ri_1 ri_2 ri_3 ri_4 ri_5 ri_6
#> ri_1 0.13205 0.07631 -0.03146 0.003062 -0.018152 0.0022289
#> ri_2 0.99877 0.04420 -0.01699 0.001533 -0.009862 0.0008388
#> ri_3 -0.42101 -0.39300 0.04230 -0.007313 0.009632 -0.0122545
#> ri_4 0.23205 0.20075 -0.97922 0.001319 -0.001499 0.0022769
#> ri_5 -0.62516 -0.58706 0.58616 -0.516641 0.006384 -0.0024638
#> ri_6 0.09654 0.06279 -0.93778 0.986829 -0.485301 0.0040370
#>
#> Fixed Effects Estimates:
#> attribute estimate SE z_test p_value lower upper
#> 1 -0.09772 0.13777 -0.7093 4.781e-01 -0.3678 0.172310
#> 2 -0.17557 0.08620 -2.0367 4.168e-02 -0.3445 -0.006617
#> 3 0.31868 0.08407 3.7909 1.501e-04 0.1539 0.483451
#> 4 0.52720 0.03335 15.8069 2.786e-56 0.4618 0.592565
#> 5 -0.41756 0.04590 -9.0968 9.307e-20 -0.5075 -0.327592
#> 6 -0.40072 0.04182 -9.5817 9.543e-22 -0.4827 -0.318752
#>
#> Model Fit Statistics:
#> logLik Dev AIC BIC AICc
#> 18.63 -37.25 16.75 63.67 29.48
#>
#> Q Error: 230.642 (42), p < 0.0001
#>
#> I2 (General):
#> names values
#> ri_1 94.53
#> ri_2 85.26
#> ri_3 84.70
#> ri_4 14.71
#> ri_5 45.51
#> ri_6 34.56
#>
#> I2 (Jackson):
#> names values
#> ri_1 90.98
#> ri_2 77.93
#> ri_3 81.33
#> ri_4 16.14
#> ri_5 43.25
#> ri_6 31.43
#>
#> I2 (Between): 83.3943
#>
#> Residual Diagnostics:
#> n n_finite_raw mean_raw sd_raw rmse mae q_pearson mean_abs_studentized
#> 48 48 -0.02155 0.2127 0.2116 0.1646 133.5 1.419
#> max_abs_studentized prop_abs_studentized_gt2 prop_abs_studentized_gt3
#> 5.006 0.2708 0.1042
#>
#> Normality (whitened residuals): test n_tested statistic p_value
#> shapiro_wilk_whitened 48 0.9789 0.5353
#>
#> Heteroscedasticity trend (|raw residual| ~ fitted): n corr_abs_raw_fitted slope p_value
#> 48 -0.002888 -0.001089 0.9845
pb_mv <- publication_bias(
fit_mv,
method = "beta_selection",
p_value_tail = "two.sided",
maxit = 100,
integration_points = 31
)
summary(pb_mv)
#> Publication Bias Model (MARS)
#> Method: beta_selection
#> Structure: UN
#> Effects: 48 | Studies: 8
#> P-value tail: two.sided
#> Convergence: 0 - Converged
#> Objective: -24.91
#>
#> Selection Parameters:
#> parameter estimate
#> alpha 1.0157
#> beta 1.0412
#>
#> Selection Weight Profile (p-value grid):
#> model p_value weight
#> beta_selection 0.0100 0.9301
#> beta_selection 0.1189 0.9622
#> beta_selection 0.2278 0.9668
#> beta_selection 0.3367 0.9666
#> beta_selection 0.4456 0.9637
#> beta_selection 0.5544 0.9584
#> beta_selection 0.6633 0.9500
#> beta_selection 0.7722 0.9371
#> beta_selection 0.8811 0.9142
#> beta_selection 0.9900 0.8270
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> factor(data[[effectID]])1 -0.1147 2.0331 -0.0564 0.9550
#> factor(data[[effectID]])2 -0.1876 1.2242 -0.1533 0.8782
#> factor(data[[effectID]])3 0.3278 0.6447 0.5085 0.6111
#> factor(data[[effectID]])4 0.5284 0.8176 0.6463 0.5181
#> factor(data[[effectID]])5 -0.4196 0.7466 -0.5621 0.5741
#> factor(data[[effectID]])6 -0.4062 0.3007 -1.3507 0.1768
#>
#> Adjusted vs Unadjusted Summary:
#> term unadjusted_estimate unadjusted_se
#> factor(data[[effectID]])1 -0.0977 0.1378
#> factor(data[[effectID]])2 -0.1756 0.0862
#> factor(data[[effectID]])3 0.3187 0.0841
#> factor(data[[effectID]])4 0.5272 0.0334
#> factor(data[[effectID]])5 -0.4176 0.0459
#> factor(data[[effectID]])6 -0.4007 0.0418
#> unadjusted_ci_lower unadjusted_ci_upper estimate std_error adjusted_ci_lower
#> -0.3678 0.1723 -0.1147 2.0331 -4.0997
#> -0.3445 -0.0066 -0.1876 1.2242 -2.5870
#> 0.1539 0.4835 0.3278 0.6447 -0.9358
#> 0.4618 0.5926 0.5284 0.8176 -1.0740
#> -0.5075 -0.3276 -0.4196 0.7466 -1.8830
#> -0.4827 -0.3188 -0.4062 0.3007 -0.9957
#> adjusted_ci_upper estimate_diff percent_change direction_changed ci_overlap
#> 3.8702 -0.0170 -17.4097 FALSE TRUE
#> 2.2117 -0.0121 -6.8720 FALSE TRUE
#> 1.5915 0.0092 2.8723 FALSE TRUE
#> 2.1309 0.0012 0.2356 FALSE TRUE
#> 1.0437 -0.0021 -0.5008 FALSE TRUE
#> 0.1832 -0.0055 -1.3739 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1S, naive-t)
#> Intercept: -3.8099 SE: 1.2322 t: -3.092 df: 7 p: 0.0175
#>
#> SE note:
#> SE rule: reported SEs use the larger of the Hessian-based and OPG-based values when both are finite, otherwise the finite one. Covariance structure: used the fitted MARS within/between-study covariance blocks directly for the publication-bias likelihood.
pb_mv_bb <- publication_bias(
fit_mv,
method = "beta_binomial",
n_bins = 8,
p_value_tail = "two.sided",
maxit = 100,
integration_points = 31
)
summary(pb_mv_bb)
#> Publication Bias Model (MARS)
#> Method: beta_binomial
#> Structure: UN
#> Effects: 48 | Studies: 8
#> P-value tail: two.sided
#> Convergence: 0 - Converged
#> Objective: -25.85
#>
#> Selection Parameters:
#> parameter estimate
#> alpha 1.6265
#> beta 1.8036
#> n_bins 8.0000
#>
#> Selection Weight Profile (bin probabilities):
#> bin p_lower p_upper weight
#> 1 0.000 0.125 0.0908
#> 2 0.125 0.250 0.1325
#> 3 0.250 0.375 0.1535
#> 4 0.375 0.500 0.1599
#> 5 0.500 0.625 0.1540
#> 6 0.625 0.750 0.1366
#> 7 0.750 0.875 0.1077
#> 8 0.875 1.000 0.0650
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> factor(data[[effectID]])1 -0.0982 0.3344 -0.2937 0.7690
#> factor(data[[effectID]])2 -0.1849 0.1652 -1.1192 0.2630
#> factor(data[[effectID]])3 0.3379 0.0236 14.3232 0.0000
#> factor(data[[effectID]])4 0.5261 0.2708 1.9426 0.0521
#> factor(data[[effectID]])5 -0.4174 0.1167 -3.5777 0.0003
#> factor(data[[effectID]])6 -0.4125 0.1421 -2.9035 0.0037
#>
#> Adjusted vs Unadjusted Summary:
#> term unadjusted_estimate unadjusted_se
#> factor(data[[effectID]])1 -0.0977 0.1378
#> factor(data[[effectID]])2 -0.1756 0.0862
#> factor(data[[effectID]])3 0.3187 0.0841
#> factor(data[[effectID]])4 0.5272 0.0334
#> factor(data[[effectID]])5 -0.4176 0.0459
#> factor(data[[effectID]])6 -0.4007 0.0418
#> unadjusted_ci_lower unadjusted_ci_upper estimate std_error adjusted_ci_lower
#> -0.3678 0.1723 -0.0982 0.3344 -0.7536
#> -0.3445 -0.0066 -0.1849 0.1652 -0.5086
#> 0.1539 0.4835 0.3379 0.0236 0.2917
#> 0.4618 0.5926 0.5261 0.2708 -0.0047
#> -0.5075 -0.3276 -0.4174 0.1167 -0.6461
#> -0.4827 -0.3188 -0.4125 0.1421 -0.6910
#> adjusted_ci_upper estimate_diff percent_change direction_changed ci_overlap
#> 0.5572 -0.0005 -0.4901 FALSE TRUE
#> 0.1389 -0.0093 -5.3047 FALSE TRUE
#> 0.3842 0.0192 6.0364 FALSE TRUE
#> 1.0568 -0.0011 -0.2147 FALSE TRUE
#> -0.1887 0.0001 0.0329 FALSE TRUE
#> -0.1341 -0.0118 -2.9453 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1S, naive-t)
#> Intercept: -3.8099 SE: 1.2322 t: -3.092 df: 7 p: 0.0175
#>
#> SE note:
#> SE rule: reported SEs use the larger of the Hessian-based and OPG-based values when both are finite, otherwise the finite one. Covariance structure: used the fitted MARS within/between-study covariance blocks directly for the publication-bias likelihood. beta_binomial uses discretized p-value bins; coefficient magnitudes can differ materially from beta_selection.
pb_mv_logit <- publication_bias(
fit_mv,
method = "logistic_selection",
p_value_tail = "two.sided",
maxit = 100,
integration_points = 31
)
summary(pb_mv_logit)
#> Publication Bias Model (MARS)
#> Method: logistic_selection
#> Structure: UN
#> Effects: 48 | Studies: 8
#> P-value tail: two.sided
#> Convergence: 0 - Converged
#> Objective: -25.6
#>
#> Selection Parameters:
#> parameter estimate
#> logit_intercept -22.1202
#> logit_slope 379.3210
#>
#> Selection Weight Profile (p-value grid):
#> model p_value weight
#> logistic_selection 0.0100 1
#> logistic_selection 0.1189 1
#> logistic_selection 0.2278 1
#> logistic_selection 0.3367 1
#> logistic_selection 0.4456 1
#> logistic_selection 0.5544 1
#> logistic_selection 0.6633 1
#> logistic_selection 0.7722 1
#> logistic_selection 0.8811 1
#> logistic_selection 0.9900 0
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> factor(data[[effectID]])1 -0.1116 0.1449 -0.7704 0.4410
#> factor(data[[effectID]])2 -0.1909 0.0899 -2.1241 0.0337
#> factor(data[[effectID]])3 0.3173 0.0798 3.9763 0.0001
#> factor(data[[effectID]])4 0.5279 0.0903 5.8436 0.0000
#> factor(data[[effectID]])5 -0.4172 0.0656 -6.3553 0.0000
#> factor(data[[effectID]])6 -0.4010 0.0585 -6.8536 0.0000
#>
#> Adjusted vs Unadjusted Summary:
#> term unadjusted_estimate unadjusted_se
#> factor(data[[effectID]])1 -0.0977 0.1378
#> factor(data[[effectID]])2 -0.1756 0.0862
#> factor(data[[effectID]])3 0.3187 0.0841
#> factor(data[[effectID]])4 0.5272 0.0334
#> factor(data[[effectID]])5 -0.4176 0.0459
#> factor(data[[effectID]])6 -0.4007 0.0418
#> unadjusted_ci_lower unadjusted_ci_upper estimate std_error adjusted_ci_lower
#> -0.3678 0.1723 -0.1116 0.1449 -0.3957
#> -0.3445 -0.0066 -0.1909 0.0899 -0.3671
#> 0.1539 0.4835 0.3173 0.0798 0.1609
#> 0.4618 0.5926 0.5279 0.0903 0.3508
#> -0.5075 -0.3276 -0.4172 0.0656 -0.5458
#> -0.4827 -0.3188 -0.4010 0.0585 -0.5157
#> adjusted_ci_upper estimate_diff percent_change direction_changed ci_overlap
#> 0.1724 -0.0139 -14.2464 FALSE TRUE
#> -0.0148 -0.0154 -8.7532 FALSE TRUE
#> 0.4737 -0.0014 -0.4281 FALSE TRUE
#> 0.7049 0.0007 0.1268 FALSE TRUE
#> -0.2885 0.0004 0.0970 FALSE TRUE
#> -0.2864 -0.0003 -0.0819 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1S, naive-t)
#> Intercept: -3.8099 SE: 1.2322 t: -3.092 df: 7 p: 0.0175
#>
#> SE note:
#> SE rule: reported SEs use the larger of the Hessian-based and OPG-based values when both are finite, otherwise the finite one. Covariance structure: used the fitted MARS within/between-study covariance blocks directly for the publication-bias likelihood.Multivariate coefficients with standard errors:
pb_mv$covariance_mode
#> [1] "fixed_mars_structure"
pb_mv$covariance_structure
#> $structure
#> [1] "UN"
#>
#> $Tau
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.132049322 0.0763072125 -0.031464964 0.003062058 -0.018151847
#> [2,] 0.076307212 0.0442044023 -0.016994035 0.001532691 -0.009862248
#> [3,] -0.031464964 -0.0169940346 0.042299545 -0.007313398 0.009632500
#> [4,] 0.003062058 0.0015326912 -0.007313398 0.001318676 -0.001499049
#> [5,] -0.018151847 -0.0098622484 0.009632500 -0.001499049 0.006384347
#> [6,] 0.002228934 0.0008387928 -0.012254466 0.002276871 -0.002463753
#> [,6]
#> [1,] 0.0022289341
#> [2,] 0.0008387928
#> [3,] -0.0122544658
#> [4,] 0.0022768707
#> [5,] -0.0024637529
#> [6,] 0.0040369638
pb_mv$coefficients
#> term estimate std_error z_value p_value
#> 1 factor(data[[effectID]])1 -0.1147351 2.0331230 -0.05643295 0.9549969
#> 2 factor(data[[effectID]])2 -0.1876318 1.2241533 -0.15327475 0.8781816
#> 3 factor(data[[effectID]])3 0.3278376 0.6447066 0.50850669 0.6110981
#> 4 factor(data[[effectID]])4 0.5284380 0.8175725 0.64635000 0.5180527
#> 5 factor(data[[effectID]])5 -0.4196498 0.7465865 -0.56209131 0.5740538
#> 6 factor(data[[effectID]])6 -0.4062260 0.3007483 -1.35071747 0.1767860
pb_mv_bb$coefficients
#> term estimate std_error z_value p_value
#> 1 factor(data[[effectID]])1 -0.09820103 0.33439228 -0.2936701 7.690100e-01
#> 2 factor(data[[effectID]])2 -0.18488000 0.16518513 -1.1192291 2.630424e-01
#> 3 factor(data[[effectID]])3 0.33792092 0.02359257 14.3231928 1.567529e-46
#> 4 factor(data[[effectID]])4 0.52606395 0.27080513 1.9425923 5.206544e-02
#> 5 factor(data[[effectID]])5 -0.41742137 0.11667222 -3.5777272 3.465949e-04
#> 6 factor(data[[effectID]])6 -0.41252285 0.14207659 -2.9035244 3.689882e-03
pb_mv_logit$coefficients
#> term estimate std_error z_value p_value
#> 1 factor(data[[effectID]])1 -0.1116439 0.14491129 -0.7704293 4.410453e-01
#> 2 factor(data[[effectID]])2 -0.1909345 0.08988834 -2.1241292 3.365934e-02
#> 3 factor(data[[effectID]])3 0.3173197 0.07980366 3.9762543 7.000927e-05
#> 4 factor(data[[effectID]])4 0.5278644 0.09033241 5.8435772 5.109160e-09
#> 5 factor(data[[effectID]])5 -0.4171534 0.06563874 -6.3552921 2.080310e-10
#> 6 factor(data[[effectID]])6 -0.4010487 0.05851647 -6.8536040 7.201224e-12The same idea applies to multivariate fits: the publication-bias
likelihood is evaluated using the fitted MARS within-study covariance
blocks and the fitted between-study Tau matrix, rather than
replacing them with a single scalar heterogeneity term.
ref_mv <- extract_full_data_coefficients(fit_mv)
cmp_mv <- compare_pub_bias_3(pb_mv, pb_mv_bb, pb_mv_logit, ref = ref_mv)
cmp_mv
#> term estimate_beta_selection se_beta_selection
#> 1 factor(data[[effectID]])1 -0.1147351 2.0331230
#> 2 factor(data[[effectID]])2 -0.1876318 1.2241533
#> 3 factor(data[[effectID]])3 0.3278376 0.6447066
#> 4 factor(data[[effectID]])4 0.5284380 0.8175725
#> 5 factor(data[[effectID]])5 -0.4196498 0.7465865
#> 6 factor(data[[effectID]])6 -0.4062260 0.3007483
#> estimate_beta_binomial se_beta_binomial estimate_logistic_selection
#> 1 -0.09820103 0.33439228 -0.1116439
#> 2 -0.18488000 0.16518513 -0.1909345
#> 3 0.33792092 0.02359257 0.3173197
#> 4 0.52606395 0.27080513 0.5278644
#> 5 -0.41742137 0.11667222 -0.4171534
#> 6 -0.41252285 0.14207659 -0.4010487
#> se_logistic_selection estimate_all_data se_all_data
#> 1 0.14491129 -0.09772205 0.13777405
#> 2 0.08988834 -0.17556678 0.08620031
#> 3 0.07980366 0.31868398 0.08406615
#> 4 0.09033241 0.52719585 0.03335216
#> 5 0.06563874 -0.41755862 0.04590195
#> 6 0.05851647 -0.40072034 0.04182127
plot_pub_bias_compare_3(cmp_mv, main = "Multivariate: publication-bias methods with all-data reference")
Jackknife version of the multivariate complete-data fit and each publication-bias model:
set.seed(2026)
becker09_jack <- becker09[seq_len(min(6, nrow(becker09))), , drop = FALSE]
jack_mv_complete <- mars_alt_estimation(
type = "jackknife",
data = becker09_jack,
studyID = "ID",
effectID = "numID",
sample_size = "N",
effectsize_type = "cor",
varcov_type = "weighted",
variable_names = c(
"Cognitive_Performance",
"Somatic_Performance",
"Selfconfidence_Performance",
"Somatic_Cognitive",
"Selfconfidence_Cognitive",
"Selfconfidence_Somatic"
),
structure = "UN",
estimation_method = "MLE",
jackknife_target = "complete_data"
)
jack_mv_beta_selection <- mars_alt_estimation(
type = "jackknife",
data = becker09_jack,
studyID = "ID",
effectID = "numID",
sample_size = "N",
effectsize_type = "cor",
varcov_type = "weighted",
variable_names = c(
"Cognitive_Performance",
"Somatic_Performance",
"Selfconfidence_Performance",
"Somatic_Cognitive",
"Selfconfidence_Cognitive",
"Selfconfidence_Somatic"
),
structure = "UN",
estimation_method = "MLE",
jackknife_target = "beta_selection",
publication_bias_args = list(
p_value_tail = "two.sided",
maxit = 80,
integration_points = 31
)
)
jack_mv_beta_binomial <- mars_alt_estimation(
type = "jackknife",
data = becker09_jack,
studyID = "ID",
effectID = "numID",
sample_size = "N",
effectsize_type = "cor",
varcov_type = "weighted",
variable_names = c(
"Cognitive_Performance",
"Somatic_Performance",
"Selfconfidence_Performance",
"Somatic_Cognitive",
"Selfconfidence_Cognitive",
"Selfconfidence_Somatic"
),
structure = "UN",
estimation_method = "MLE",
jackknife_target = "beta_binomial",
publication_bias_args = list(
p_value_tail = "two.sided",
n_bins = 8,
maxit = 80,
integration_points = 31
)
)
jack_mv_logistic <- mars_alt_estimation(
type = "jackknife",
data = becker09_jack,
studyID = "ID",
effectID = "numID",
sample_size = "N",
effectsize_type = "cor",
varcov_type = "weighted",
variable_names = c(
"Cognitive_Performance",
"Somatic_Performance",
"Selfconfidence_Performance",
"Somatic_Cognitive",
"Selfconfidence_Cognitive",
"Selfconfidence_Somatic"
),
structure = "UN",
estimation_method = "MLE",
jackknife_target = "logistic_selection",
publication_bias_args = list(
p_value_tail = "two.sided",
maxit = 80,
integration_points = 31
)
)
jack_mv_complete_stats <- compute_alt_stats(jack_mv_complete, type = "jackknife")
jack_mv_beta_selection_stats <- compute_alt_stats(jack_mv_beta_selection, type = "jackknife")
jack_mv_beta_binomial_stats <- compute_alt_stats(jack_mv_beta_binomial, type = "jackknife")
jack_mv_logistic_stats <- compute_alt_stats(jack_mv_logistic, type = "jackknife")
mv_term <- pb_mv$coefficients$term[1]
jackknife_compare_mv <- rbind(
complete_data = extract_alt_interval(jack_mv_complete_stats, term = mv_term),
beta_selection = extract_alt_interval(jack_mv_beta_selection_stats, term = mv_term),
beta_binomial = extract_alt_interval(jack_mv_beta_binomial_stats, term = mv_term),
logistic_selection = extract_alt_interval(jack_mv_logistic_stats, term = mv_term)
)
jackknife_compare_mvThe multivariate jackknife illustration also uses a small subset so the chunk demonstrates the interface without dominating vignette runtime.
jack_mv_est <- jackknife_compare_mv[, "point_estimate"]
jack_mv_lo <- jackknife_compare_mv[, "ci_lower"]
jack_mv_hi <- jackknife_compare_mv[, "ci_upper"]
xpos_mv <- seq_len(nrow(jackknife_compare_mv))
plot(
xpos_mv, jack_mv_est,
ylim = range(c(jack_mv_lo, jack_mv_hi)),
xaxt = "n", xlab = "", ylab = "Jackknife Estimate",
pch = 19, col = "#1F78B4",
main = paste("Multivariate Jackknife for", mv_term)
)
segments(xpos_mv, jack_mv_lo, xpos_mv, jack_mv_hi, lwd = 2, col = "#1F78B4")
abline(h = 0, lty = 2, col = "gray60")
axis(1, at = xpos_mv, labels = rownames(jackknife_compare_mv), las = 2)
title(xlab = "Jackknife Target", line = 6)Method comparison note: - beta_selection and
beta_binomial use different selection-function
parameterizations. - Coefficients and SEs can differ materially,
especially when beta_binomial uses coarse binning or when
Hessian diagnostics indicate weak identifiability. - The
all_data reference intervals come from the original MARS
fit, while publication-bias intervals come from the publication-bias
likelihood with the conservative Hessian-versus-OPG SE rule. - In
multilevel and multivariate settings, publication-bias fits now retain
the fitted MARS covariance structure through
covariance_mode = "fixed_mars_structure".
The same publication-bias interface can be used across univariate,
multilevel, and multivariate mars objects, and both methods
support funnel plotting.