Publication Bias Models with MARS
mars authors
2026-04-02
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.1.1
#> Thursday, April 2, 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"
)
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.24
#>
#> Selection Parameters:
#> parameter estimate
#> alpha 1.0330
#> beta 0.2139
#> tau2 0.5664
#>
#> Selection Weight Profile (p-value grid):
#> model p_value weight
#> beta_selection 0.0100 0.1870
#> beta_selection 0.1189 0.2224
#> beta_selection 0.2278 0.2520
#> beta_selection 0.3367 0.2877
#> beta_selection 0.4456 0.3343
#> beta_selection 0.5544 0.3999
#> beta_selection 0.6633 0.5014
#> beta_selection 0.7722 0.6852
#> beta_selection 0.8811 1.1473
#> beta_selection 0.9900 8.0641
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> (Intercept) -1.0912 0.2881 -3.7873 2e-04
#>
#> 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 -1.0912 0.2881 -1.6559 -0.5265
#> estimate_diff percent_change direction_changed ci_overlap
#> -0.38 -53.4245 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1)
#> 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"
)
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: 12.4
#>
#> Selection Parameters:
#> parameter estimate
#> alpha 0.5407
#> beta 0.2341
#> tau2 0.3591
#> n_bins 8.0000
#>
#> Selection Weight Profile (bin probabilities):
#> bin p_lower p_upper weight
#> 1 0.000 0.125 0.1072
#> 2 0.125 0.250 0.0651
#> 3 0.250 0.375 0.0575
#> 4 0.375 0.500 0.0575
#> 5 0.500 0.625 0.0629
#> 6 0.625 0.750 0.0767
#> 7 0.750 0.875 0.1148
#> 8 0.875 1.000 0.4583
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> (Intercept) -0.6824 0.1927 -3.5419 4e-04
#>
#> 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.6824 0.1927 -1.0601 -0.3048
#> estimate_diff percent_change direction_changed ci_overlap
#> 0.0288 4.0426 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1)
#> 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"
)
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: 13.92
#>
#> Selection Parameters:
#> parameter estimate
#> logit_intercept 9.0929
#> logit_slope 0.4472
#> tau2 0.4256
#>
#> 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 0.9999
#> logistic_selection 0.3367 0.9999
#> logistic_selection 0.4456 0.9999
#> logistic_selection 0.5544 0.9999
#> logistic_selection 0.6633 0.9999
#> logistic_selection 0.7722 0.9999
#> logistic_selection 0.8811 0.9999
#> logistic_selection 0.9900 0.9999
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> (Intercept) -0.7724 0.3054 -2.5287 0.0114
#>
#> 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.7724 0.3054 -1.371 -0.1737
#> estimate_diff percent_change direction_changed ci_overlap
#> -0.0612 -8.5986 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1)
#> 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.0330169 0.2138633 0.5663630
pb_uni$covariance_mode
#> [1] "publication_bias_tau2"
pb_uni$coefficients
#> term estimate std_error z_value p_value
#> 1 (Intercept) -1.091154 0.2881115 -3.787261 0.0001523172
pb_uni_bb$selection_weights
#> bin p_lower p_upper weight
#> 1 1 0.000 0.125 0.10719722
#> 2 2 0.125 0.250 0.06507810
#> 3 3 0.250 0.375 0.05746743
#> 4 4 0.375 0.500 0.05747199
#> 5 5 0.500 0.625 0.06291976
#> 6 6 0.625 0.750 0.07672810
#> 7 7 0.750 0.875 0.11482686
#> 8 8 0.875 1.000 0.45831054
pb_uni_logit$selection
#> logit_intercept logit_slope tau2
#> 9.0929409 0.4471942 0.4255670
pb_uni_logit$coefficients
#> term estimate std_error z_value p_value
#> 1 (Intercept) -0.7723525 0.3054325 -2.528718 0.01144801For 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) -1.091154 0.2881115 -0.682448
#> se_beta_binomial estimate_logistic_selection se_logistic_selection
#> 1 0.1926782 -0.7723525 0.3054325
#> 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
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_compare
#> point_estimate ci_lower ci_upper
#> complete_data -0.7071685 -0.7867671 -0.6292833
#> beta_selection -1.1713117 -4.1793912 -0.7214840
#> beta_binomial -0.6908430 -1.7810121 0.6252479
#> logistic_selection -0.7525089 -0.8203594 -0.6654832
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.1.1
#> Thursday, April 2, 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"
)
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: -111
#>
#> Selection Parameters:
#> parameter estimate
#> alpha 1.0047
#> beta 0.0034
#>
#> Selection Weight Profile (p-value grid):
#> model p_value weight
#> beta_selection 0.0100 0.0034
#> beta_selection 0.1189 0.0039
#> beta_selection 0.2278 0.0044
#> beta_selection 0.3367 0.0051
#> beta_selection 0.4456 0.0061
#> beta_selection 0.5544 0.0077
#> beta_selection 0.6633 0.0101
#> beta_selection 0.7722 0.0150
#> beta_selection 0.8811 0.0286
#> beta_selection 0.9900 0.3375
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> (Intercept) 0.3259 9e-04 359.3879 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.3259 9e-04 0.3241 0.3277
#> estimate_diff percent_change direction_changed ci_overlap
#> 0.1415 76.6906 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1)
#> 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"
)
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: -302.1
#>
#> Selection Parameters:
#> parameter estimate
#> alpha 122.3738
#> beta 0.3659
#> 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.0006
#> 7 0.750 0.875 0.0195
#> 8 0.875 1.000 0.9798
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> (Intercept) -1.3124 0.0201 -65.3197 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.3124 0.0201 -1.3518 -1.273
#> estimate_diff percent_change direction_changed ci_overlap
#> -1.4968 -811.496 TRUE FALSE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1)
#> 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"
)
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.395
#>
#> Selection Parameters:
#> parameter estimate
#> logit_intercept 11.7447
#> logit_slope 12.3280
#>
#> 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 1
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> (Intercept) 0.1845 0.0871 2.1186 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.1845 0.0871 0.0138 0.3551
#> estimate_diff percent_change direction_changed ci_overlap
#> 0 -2e-04 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1)
#> 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.3259153 0.000906862 359.3879 0
pb_ml_bb$coefficients
#> term estimate std_error z_value p_value
#> 1 (Intercept) -1.312393 0.02009184 -65.31968 0
pb_ml_logit$coefficients
#> term estimate std_error z_value p_value
#> 1 (Intercept) 0.1844551 0.0870634 2.11863 0.03412178For 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.3259153 0.000906862 -1.312393
#> se_beta_binomial estimate_logistic_selection se_logistic_selection
#> 1 0.02009184 0.1844551 0.0870634
#> 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_ml
#> point_estimate ci_lower ci_upper
#> complete_data 0.23973526 0.14570705 0.3043585
#> beta_selection 0.04208093 0.01522397 0.1019976
#> beta_binomial 0.20373586 0.08254856 0.3298935
#> logistic_selection 0.16972049 0.08100785 0.1962584The 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.1.1
#> Thursday, April 2, 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.13204 0.07630 -0.03144 0.003058 -0.018150 0.0022236
#> ri_2 0.99878 0.04420 -0.01698 0.001532 -0.009861 0.0008355
#> ri_3 -0.42074 -0.39284 0.04229 -0.007313 0.009630 -0.0122546
#> ri_4 0.23176 0.20062 -0.97927 0.001319 -0.001498 0.0022776
#> ri_5 -0.62514 -0.58701 0.58603 -0.516444 0.006384 -0.0024641
#> ri_6 0.09632 0.06256 -0.93797 0.987249 -0.485426 0.0040360
#>
#> Fixed Effects Estimates:
#> attribute estimate SE z_test p_value lower upper
#> 1 -0.09773 0.13777 -0.7093 4.781e-01 -0.3677 0.172299
#> 2 -0.17556 0.08620 -2.0367 4.168e-02 -0.3445 -0.006613
#> 3 0.31868 0.08406 3.7911 1.500e-04 0.1539 0.483438
#> 4 0.52720 0.03335 15.8070 2.785e-56 0.4618 0.592564
#> 5 -0.41756 0.04590 -9.0968 9.305e-20 -0.5075 -0.327591
#> 6 -0.40071 0.04182 -9.5821 9.510e-22 -0.4827 -0.318750
#>
#> 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.69
#> 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.42
#>
#> I2 (Between): 83.393
#>
#> 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.5351
#>
#> Heteroscedasticity trend (|raw residual| ~ fitted): n corr_abs_raw_fitted slope p_value
#> 48 -0.002893 -0.001091 0.9844
pb_mv <- publication_bias(
fit_mv,
method = "beta_selection",
p_value_tail = "two.sided"
)
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: -17.08
#>
#> Selection Parameters:
#> parameter estimate
#> alpha 0.9761
#> beta 0.8709
#>
#> Selection Weight Profile (p-value grid):
#> model p_value weight
#> beta_selection 0.0100 0.9522
#> beta_selection 0.1189 0.9111
#> beta_selection 0.2278 0.9125
#> beta_selection 0.3367 0.9219
#> beta_selection 0.4456 0.9372
#> beta_selection 0.5544 0.9590
#> beta_selection 0.6633 0.9901
#> beta_selection 0.7722 1.0375
#> beta_selection 0.8811 1.1248
#> beta_selection 0.9900 1.5442
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> factor(data[[effectID]])1 -0.0573 0.1721 -0.3330 0.7391
#> factor(data[[effectID]])2 -0.2066 0.1568 -1.3172 0.1878
#> factor(data[[effectID]])3 0.2559 0.1346 1.9009 0.0573
#> factor(data[[effectID]])4 0.5083 0.2685 1.8933 0.0583
#> factor(data[[effectID]])5 -0.4056 0.2949 -1.3752 0.1691
#> factor(data[[effectID]])6 -0.3925 0.0697 -5.6313 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.0573 0.1721 -0.3946
#> -0.3445 -0.0066 -0.2066 0.1568 -0.5139
#> 0.1539 0.4834 0.2559 0.1346 -0.0080
#> 0.4618 0.5926 0.5083 0.2685 -0.0179
#> -0.5075 -0.3276 -0.4056 0.2949 -0.9837
#> -0.4827 -0.3187 -0.3925 0.0697 -0.5291
#> adjusted_ci_upper estimate_diff percent_change direction_changed ci_overlap
#> 0.2800 0.0404 41.3524 FALSE TRUE
#> 0.1008 -0.0310 -17.6571 FALSE TRUE
#> 0.5198 -0.0628 -19.6984 FALSE TRUE
#> 1.0345 -0.0189 -3.5841 FALSE TRUE
#> 0.1725 0.0120 2.8643 FALSE TRUE
#> -0.2559 0.0082 2.0547 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1)
#> 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"
)
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: -17.57
#>
#> Selection Parameters:
#> parameter estimate
#> alpha 1.2740
#> beta 0.9773
#> n_bins 8.0000
#>
#> Selection Weight Profile (bin probabilities):
#> bin p_lower p_upper weight
#> 1 0.000 0.125 0.0777
#> 2 0.125 0.250 0.0993
#> 3 0.250 0.375 0.1134
#> 4 0.375 0.500 0.1243
#> 5 0.500 0.625 0.1336
#> 6 0.625 0.750 0.1420
#> 7 0.750 0.875 0.1501
#> 8 0.875 1.000 0.1596
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> factor(data[[effectID]])1 -0.1744 0.2298 -0.7589 0.4479
#> factor(data[[effectID]])2 -0.2887 0.1571 -1.8379 0.0661
#> factor(data[[effectID]])3 0.3471 0.0878 3.9538 0.0001
#> factor(data[[effectID]])4 0.5116 0.1128 4.5341 0.0000
#> factor(data[[effectID]])5 -0.4142 0.1094 -3.7847 0.0002
#> factor(data[[effectID]])6 -0.4337 0.0530 -8.1894 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.1744 0.2298 -0.6249
#> -0.3445 -0.0066 -0.2887 0.1571 -0.5967
#> 0.1539 0.4834 0.3471 0.0878 0.1750
#> 0.4618 0.5926 0.5116 0.1128 0.2905
#> -0.5075 -0.3276 -0.4142 0.1094 -0.6286
#> -0.4827 -0.3187 -0.4337 0.0530 -0.5376
#> adjusted_ci_upper estimate_diff percent_change direction_changed ci_overlap
#> 0.2761 -0.0767 -78.4704 FALSE TRUE
#> 0.0192 -0.1132 -64.4697 FALSE TRUE
#> 0.5191 0.0284 8.9112 FALSE TRUE
#> 0.7328 -0.0155 -2.9493 FALSE TRUE
#> -0.1997 0.0034 0.8131 FALSE TRUE
#> -0.3299 -0.0330 -8.2440 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1)
#> 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"
)
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: -16.75
#>
#> Selection Parameters:
#> parameter estimate
#> logit_intercept 3.5331
#> logit_slope 23.1095
#>
#> 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 0.9999
#> logistic_selection 0.8811 0.9984
#> logistic_selection 0.9900 0.9774
#>
#> Bias-Adjusted Coefficients:
#> term estimate std_error z_value p_value
#> factor(data[[effectID]])1 -0.0991 0.2421 -0.4095 0.6822
#> factor(data[[effectID]])2 -0.2357 0.1671 -1.4110 0.1583
#> factor(data[[effectID]])3 0.2973 0.1220 2.4359 0.0149
#> factor(data[[effectID]])4 0.5148 0.0662 7.7805 0.0000
#> factor(data[[effectID]])5 -0.4169 0.0808 -5.1584 0.0000
#> factor(data[[effectID]])6 -0.4183 0.0615 -6.7981 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.0991 0.2421 -0.5736
#> -0.3445 -0.0066 -0.2357 0.1671 -0.5632
#> 0.1539 0.4834 0.2973 0.1220 0.0581
#> 0.4618 0.5926 0.5148 0.0662 0.3852
#> -0.5075 -0.3276 -0.4169 0.0808 -0.5753
#> -0.4827 -0.3187 -0.4183 0.0615 -0.5389
#> adjusted_ci_upper estimate_diff percent_change direction_changed ci_overlap
#> 0.3753 -0.0014 -1.4307 FALSE TRUE
#> 0.0917 -0.0602 -34.2677 FALSE TRUE
#> 0.5365 -0.0214 -6.7152 FALSE TRUE
#> 0.6445 -0.0123 -2.3418 FALSE TRUE
#> -0.2585 0.0006 0.1532 FALSE TRUE
#> -0.2977 -0.0176 -4.3802 FALSE TRUE
#>
#> Small-Study Effects Check (Egger-type):
#> Model: SND ~ precision
#> Method: Egger-type small-study-effects regression with study-cluster robust SE (CR1)
#> 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.132040025 0.0763014664 -0.031441429 0.003058166 -0.018150222
#> [2,] 0.076301466 0.0441998006 -0.016984791 0.001531606 -0.009860779
#> [3,] -0.031441429 -0.0169847911 0.042292885 -0.007313107 0.009629524
#> [4,] 0.003058166 0.0015316061 -0.007313107 0.001318666 -0.001498453
#> [5,] -0.018150222 -0.0098607790 0.009629524 -0.001498453 0.006384168
#> [6,] 0.002223569 0.0008355477 -0.012254561 0.002277566 -0.002464062
#> [,6]
#> [1,] 0.0022235691
#> [2,] 0.0008355477
#> [3,] -0.0122545613
#> [4,] 0.0022775661
#> [5,] -0.0024640625
#> [6,] 0.0040360241
pb_mv$coefficients
#> term estimate std_error z_value p_value
#> 1 factor(data[[effectID]])1 -0.05731339 0.17209590 -0.3330317 7.391104e-01
#> 2 factor(data[[effectID]])2 -0.20655358 0.15681039 -1.3172187 1.877653e-01
#> 3 factor(data[[effectID]])3 0.25590590 0.13462471 1.9008836 5.731726e-02
#> 4 factor(data[[effectID]])4 0.50830044 0.26847412 1.8932940 5.831878e-02
#> 5 factor(data[[effectID]])5 -0.40559668 0.29493108 -1.3752253 1.690616e-01
#> 6 factor(data[[effectID]])6 -0.39248053 0.06969686 -5.6312513 1.789068e-08
pb_mv_bb$coefficients
#> term estimate std_error z_value p_value
#> 1 factor(data[[effectID]])1 -0.1744103 0.22982856 -0.7588711 4.479297e-01
#> 2 factor(data[[effectID]])2 -0.2887357 0.15710154 -1.8378924 6.607826e-02
#> 3 factor(data[[effectID]])3 0.3470794 0.08778424 3.9537777 7.692687e-05
#> 4 factor(data[[effectID]])4 0.5116472 0.11284349 4.5341311 5.784109e-06
#> 5 factor(data[[effectID]])5 -0.4141613 0.10942945 -3.7847338 1.538732e-04
#> 6 factor(data[[effectID]])6 -0.4337486 0.05296469 -8.1893924 2.625483e-16
pb_mv_logit$coefficients
#> term estimate std_error z_value p_value
#> 1 factor(data[[effectID]])1 -0.09912321 0.24206557 -0.4094891 6.821808e-01
#> 2 factor(data[[effectID]])2 -0.23571435 0.16706042 -1.4109527 1.582586e-01
#> 3 factor(data[[effectID]])3 0.29728114 0.12204165 2.4358991 1.485483e-02
#> 4 factor(data[[effectID]])4 0.51484944 0.06617171 7.7805071 7.223443e-15
#> 5 factor(data[[effectID]])5 -0.41691705 0.08082318 -5.1583846 2.490896e-07
#> 6 factor(data[[effectID]])6 -0.41826593 0.06152684 -6.7981048 1.060043e-11The 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.05731339 0.17209590
#> 2 factor(data[[effectID]])2 -0.20655358 0.15681039
#> 3 factor(data[[effectID]])3 0.25590590 0.13462471
#> 4 factor(data[[effectID]])4 0.50830044 0.26847412
#> 5 factor(data[[effectID]])5 -0.40559668 0.29493108
#> 6 factor(data[[effectID]])6 -0.39248053 0.06969686
#> estimate_beta_binomial se_beta_binomial estimate_logistic_selection
#> 1 -0.1744103 0.22982856 -0.09912321
#> 2 -0.2887357 0.15710154 -0.23571435
#> 3 0.3470794 0.08778424 0.29728114
#> 4 0.5116472 0.11284349 0.51484944
#> 5 -0.4141613 0.10942945 -0.41691705
#> 6 -0.4337486 0.05296469 -0.41826593
#> se_logistic_selection estimate_all_data se_all_data
#> 1 0.24206557 -0.09772505 0.13776971
#> 2 0.16706042 -0.17555552 0.08619682
#> 3 0.12204165 0.31868112 0.08406129
#> 4 0.06617171 0.52719552 0.03335209
#> 5 0.08082318 -0.41755660 0.04590156
#> 6 0.06152684 -0.40071392 0.04181907
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_mv
#> point_estimate ci_lower ci_upper
#> complete_data -0.08983328 -0.20035300 0.01635273
#> beta_selection 0.02405734 -0.08801925 0.06878734
#> beta_binomial -0.14637973 -0.27612223 0.01877218
#> logistic_selection -0.07121728 -0.17278774 0.02747445The 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.