User-Specified Heterogeneity in MARS

mars() usually estimates between-study heterogeneity from the data. In some workflows, the analyst may want to fix heterogeneity at a known or user-selected value instead. Common examples include sensitivity analyses, comparisons against an external estimate, or network meta-analysis models that use a prespecified common heterogeneity value.

The tau2 argument fixes the between-study variance or covariance and estimates the fixed effects conditional on that value. When tau2 = NULL, the default, MARS estimates heterogeneity.

library(mars)

Univariate Meta-analysis

For univariate models, tau2 is one non-negative value.

fit_uni <- mars(
  data = teacher_expectancy,
  studyID = "study",
  effectID = NULL,
  sample_size = NULL,
  formula = yi ~ 1,
  variance = "vi",
  varcov_type = "univariate",
  structure = "univariate",
  estimation_method = "MLE",
  tau2 = 0.04
)

summary(fit_uni)
#> Results generated with MARS:v 0.5.2 
#> Friday, May 15, 2026
#> 
#> Model Type: 
#> univariate
#> 
#> Estimation Method: 
#> Maximum Likelihood
#> 
#> Model Formula: 
#> yi ~ 1
#> 
#> Data Summary: 
#> Number of Effect Sizes: 19
#> Number of Fixed Effects: 1
#> Number of Random Effects: 1
#> 
#> Random Components: 
#> Between-study variance/covariance fixed by user-supplied tau2.
#>       term  var  SD
#>  intercept 0.04 0.2
#> 
#> Fixed Effects Estimates: 
#>    attribute estimate      SE z_test p_value    lower  upper
#>  (Intercept)  0.09817 0.06296  1.559   0.119 -0.02523 0.2216
#> 
#> Model Fit Statistics: 
#>  logLik   Dev   AIC   BIC  AICc
#>  -3.622 13.24 9.244 10.19 13.86
#> 
#> Q Error: 35.83 (18), p = 0.0074
#> 
#> I2 (General): 61.82526
#> 
#> Residual Diagnostics: 
#>   n n_finite_raw mean_raw sd_raw   rmse    mae q_pearson mean_abs_studentized
#>  19           19  0.06551 0.3589 0.3554 0.2499     20.06               0.8197
#>  max_abs_studentized prop_abs_studentized_gt2 prop_abs_studentized_gt3
#>                2.585                   0.1053                        0
#> 
#> Normality (whitened residuals):                   test n_tested statistic p_value
#>  shapiro_wilk_whitened       19    0.9302  0.1747
#> 
#> Heteroscedasticity trend (|raw residual| ~ fitted):   n corr_abs_raw_fitted slope p_value
#>  19                  NA    NA      NA
fit_uni$est_values$Tau
#>      [,1]
#> [1,] 0.04

The summary notes that the random component was fixed by the user-supplied tau2 value.

Multilevel Meta-analysis

For multilevel models, supply one value per random-effect component. Named vectors are matched to the random-effect component names.

fit_ml <- mars(
  data = school,
  formula = effect ~ 1 + (1 | district/study),
  studyID = "district",
  effectID = NULL,
  sample_size = NULL,
  variance = "var",
  varcov_type = "multilevel",
  structure = "multilevel",
  estimation_method = "MLE",
  tau2 = c(district = 0.03, study = 0.01)
)

summary(fit_ml)
#> Results generated with MARS:v 0.5.2 
#> Friday, May 15, 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: 
#> Between-study variance/covariance fixed by user-supplied tau2.
#>      term  var     SD
#>  district 0.03 0.1732
#>     study 0.01 0.1000
#> 
#> Fixed Effects Estimates: 
#>    attribute estimate      SE z_test  p_value   lower  upper
#>  (Intercept)   0.1926 0.05859  3.287 0.001012 0.07776 0.3074
#> 
#> Model Fit Statistics: 
#>  logLik   Dev   AIC   BIC  AICc
#>  -15.07 38.14 32.14 34.16 33.85
#> 
#> Q Error: 578.864 (55), p < 0.0001
#> 
#> I2 (General): 
#>     names values
#>  district  86.29
#>     study  67.73
#> 
#> I2 (Jackson): 97.5396
#> I2 (Between): 73.67202
#> 
#> Residual Diagnostics: 
#>   n n_finite_raw mean_raw sd_raw  rmse   mae q_pearson mean_abs_studentized
#>  56           56 -0.07259 0.3258 0.331 0.266     202.7                1.544
#>  max_abs_studentized prop_abs_studentized_gt2 prop_abs_studentized_gt3
#>                6.034                   0.3214                  0.08929
#> 
#> Normality (whitened residuals):                   test n_tested statistic p_value
#>  shapiro_wilk_whitened       56     0.955 0.03593
#> 
#> Heteroscedasticity trend (|raw residual| ~ fitted):   n corr_abs_raw_fitted slope p_value
#>  56                  NA    NA      NA
fit_ml$est_values$Tau
#>          district study
#> district     0.03  0.00
#> study        0.00  0.01

Multivariate Meta-analysis

For multivariate models, the expected tau2 shape follows the selected between-study covariance structure:

DIAG2: one common variance value
DIAG1: one variance value per outcome
UN: a full positive semidefinite covariance matrix

The next examples use the Becker data and are shown as templates because the model can take longer on some systems.

vars <- c(
  "Cognitive_Performance",
  "Somatic_Performance",
  "Selfconfidence_Performance",
  "Somatic_Cognitive",
  "Selfconfidence_Cognitive",
  "Selfconfidence_Somatic"
)

fit_mv_common <- mars(
  data = becker09,
  studyID = "ID",
  effectID = "numID",
  sample_size = "N",
  effectsize_type = "cor",
  varcov_type = "weighted",
  variable_names = vars,
  structure = "DIAG2",
  estimation_method = "MLE",
  tau2 = 0.02
)

Tau_user <- diag(0.03, length(vars))
Tau_user[lower.tri(Tau_user)] <- 0.005
Tau_user[upper.tri(Tau_user)] <- 0.005

fit_mv_un <- mars(
  data = becker09,
  studyID = "ID",
  effectID = "numID",
  sample_size = "N",
  effectsize_type = "cor",
  varcov_type = "weighted",
  variable_names = vars,
  structure = "UN",
  estimation_method = "MLE",
  tau2 = Tau_user
)

fit_mv_un$est_values$Tau

Network Meta-analysis

network_meta() passes tau2 to the underlying MARS fit. For multivariate network models with heterogeneity = "common", supply one value. With heterogeneity = "flexible", supply one value per comparison. To fix a full multivariate heterogeneity covariance in network mode, pass a positive semidefinite matrix. For multilevel network models, supply one value per random-effect component.

nma_dat <- data.frame(
  study = c("s1", "s1", "s2", "s3", "s4", "s4", "s5", "s6"),
  trt1 = c("A", "A", "A", "A", "B", "A", "B", "A"),
  trt2 = c("B", "C", "B", "C", "C", "B", "C", "C"),
  yi = c(0.20, 0.42, 0.12, 0.48, 0.26, 0.15, 0.31, 0.44),
  vi = c(0.04, 0.05, 0.03, 0.06, 0.05, 0.04, 0.05, 0.04),
  stringsAsFactors = FALSE
)

nma_fixed <- network_meta(
  data = nma_dat,
  study_id = "study",
  treatment_1 = "trt1",
  treatment_2 = "trt2",
  effect = "yi",
  variance = "vi",
  reference = "A",
  model_type = "multivariate",
  heterogeneity = "common",
  estimation_method = "MLE",
  tau2 = 0.02
)

nma_fixed$fit$est_values$Tau
#>      [,1] [,2] [,3]
#> [1,] 0.02 0.00 0.00
#> [2,] 0.00 0.02 0.00
#> [3,] 0.00 0.00 0.02

Tau_nma <- diag(0.02, 3)
Tau_nma[lower.tri(Tau_nma)] <- 0.003
Tau_nma[upper.tri(Tau_nma)] <- 0.003

nma_fixed_cov <- network_meta(
  data = nma_dat,
  study_id = "study",
  treatment_1 = "trt1",
  treatment_2 = "trt2",
  effect = "yi",
  variance = "vi",
  reference = "A",
  model_type = "multivariate",
  heterogeneity = "flexible",
  estimation_method = "MLE",
  tau2 = Tau_nma
)

Practical Notes

Fixed tau2 values are treated as user inputs, not estimated parameters. They are therefore not assigned heterogeneity standard errors, and model comparisons should be interpreted as conditional on the supplied heterogeneity values.

tau2 cannot be combined with scale_formula, because scale-location models estimate heterogeneity as a function of predictors. It also cannot be combined with lasso = TRUE in the current implementation.

mars authors

2026-05-15

Univariate Meta-analysis

Multilevel Meta-analysis

Multivariate Meta-analysis

Network Meta-analysis

Practical Notes