MARS Model Examples

This vignette shows three common ways to use mars():

a univariate meta-analysis with one effect size per study,
a multivariate meta-analysis with multiple dependent effect sizes, and
a multilevel meta-analysis with nested effect sizes.

The examples use datasets included with the package and keep the fitted models small enough to run quickly during package installation.

library(mars)

Univariate Meta-Analysis

Use structure = "univariate" when each row contributes a single effect size with its own sampling variance. The bcg data include log risk-ratio estimates (yi) and their sampling variances (vi) from tuberculosis vaccine trials.

fit_uni <- mars(
  data = bcg,
  formula = yi ~ 1,
  studyID = "trial",
  variance = "vi",
  varcov_type = "univariate",
  structure = "univariate",
  estimation_method = "MLE"
)

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: 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

Moderators can be added directly to the formula. Here, absolute latitude is used as a study-level predictor.

fit_uni_mod <- mars(
  data = bcg,
  formula = yi ~ 1 + ablat,
  studyID = "trial",
  variance = "vi",
  varcov_type = "univariate",
  structure = "univariate",
  estimation_method = "MLE"
)

summary(fit_uni_mod)
#> Results generated with MARS:v 0.5.2 
#> Friday, May 15, 2026
#> 
#> Model Type: 
#> univariate
#> 
#> Estimation Method: 
#> Maximum Likelihood
#> 
#> Model Formula: 
#> yi ~ 1 + ablat
#> 
#> Data Summary: 
#> Number of Effect Sizes: 13
#> Number of Fixed Effects: 2
#> Number of Random Effects: 1
#> 
#> Random Components: 
#>       term     var     SD
#>  intercept 0.03435 0.1853
#> 
#> Fixed Effects Estimates: 
#>    attribute estimate       SE z_test   p_value    lower    upper
#>  (Intercept)  0.28211 0.187184  1.507 1.318e-01 -0.08477  0.64898
#>        ablat -0.02951 0.005488 -5.377 7.560e-08 -0.04027 -0.01875
#> 
#> Model Fit Statistics: 
#>  logLik   Dev   AIC   BIC  AICc
#>  -7.686 25.37 21.37 23.07 43.77
#> 
#> Q Error: 30.733 (11), p = 0.0012
#> 
#> I2 (General): 51.71087
#> 
#> Residual Diagnostics: 
#>   n n_finite_raw mean_raw sd_raw rmse    mae q_pearson mean_abs_studentized
#>  13           13 -0.03533 0.5817 0.56 0.3935     18.46                1.021
#>  max_abs_studentized prop_abs_studentized_gt2 prop_abs_studentized_gt3
#>                2.704                   0.1538                        0
#> 
#> Normality (whitened residuals):                   test n_tested statistic p_value
#>  shapiro_wilk_whitened       13    0.9727   0.924
#> 
#> Heteroscedasticity trend (|raw residual| ~ fitted):   n corr_abs_raw_fitted   slope p_value
#>  13              0.1005 0.09779  0.7439

Multivariate Meta-Analysis

Use effectsize_type = "cor" for correlations stored in a wide study-level dataset. In this workflow, variable_names identifies the correlation columns, and mars() builds the working multivariate dataset internally.

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.2 
#> Friday, May 15, 2026
#> 
#> Model Type: 
#> multivariate
#> 
#> Estimation Method: 
#> Restricted Maximum Likelihood
#> 
#> Model Formula: 
#> NULL
#> 
#> Data Summary: 
#> Number of Effect Sizes: 54
#> 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.1221  0.07342 -0.02923  0.002719 -0.015945  0.003653
#> ri_2  0.8992  0.05460 -0.02882  0.003697 -0.012822  0.001786
#> ri_3 -0.3645 -0.53747  0.05266 -0.009075  0.012947 -0.013151
#> ri_4  0.1934  0.39320 -0.98293  0.001619 -0.002136  0.002353
#> ri_5 -0.5668 -0.68157  0.70076 -0.659537  0.006482 -0.002520
#> ri_6  0.1487  0.10871 -0.81499  0.831698 -0.445197  0.004945
#> 
#> Fixed Effects Estimates: 
#>  attribute estimate      SE  z_test   p_value    lower   upper
#>          1 -0.02286 0.12024 -0.1901 8.492e-01 -0.25853  0.2128
#>          2 -0.06819 0.08607 -0.7923 4.282e-01 -0.23688  0.1005
#>          3  0.23486 0.08614  2.7264 6.403e-03  0.06602  0.4037
#>          4  0.54220 0.03456 15.6895 1.784e-55  0.47447  0.6099
#>          5 -0.44472 0.04717 -9.4278 4.187e-21 -0.53717 -0.3523
#>          6 -0.38870 0.04407 -8.8199 1.145e-18 -0.47507 -0.3023
#> 
#> Model Fit Statistics: 
#>  logLik    Dev   AIC   BIC  AICc
#>   15.39 -30.79 23.21 73.73 33.98
#> 
#> Q Error: 217.704 (48), p < 0.0001
#> 
#> I2 (General): 
#>  names values
#>   ri_1  92.91
#>   ri_2  85.43
#>   ri_3  84.98
#>   ri_4  14.81
#>   ri_5  41.05
#>   ri_6  34.69
#> 
#> I2 (Jackson): 
#>  names values
#>   ri_1  89.19
#>   ri_2  79.35
#>   ri_3  82.78
#>   ri_4  17.29
#>   ri_5  39.51
#>   ri_6  33.65
#> 
#> I2 (Between): 81.26981
#> 
#> Residual Diagnostics: 
#>   n n_finite_raw mean_raw sd_raw   rmse   mae q_pearson mean_abs_studentized
#>  54           54 -0.01122 0.2431 0.2411 0.184     109.9                1.217
#>  max_abs_studentized prop_abs_studentized_gt2 prop_abs_studentized_gt3
#>                  4.1                   0.1852                  0.05556
#> 
#> Normality (whitened residuals):                   test n_tested statistic p_value
#>  shapiro_wilk_whitened       54    0.9797  0.4891
#> 
#> Heteroscedasticity trend (|raw residual| ~ fitted):   n corr_abs_raw_fitted     slope p_value
#>  54           0.0002964 0.0001388  0.9983

The result estimates the average correlation vector and a between-study covariance structure for the dependent correlations.

Multilevel Meta-Analysis

Use structure = "multilevel" when effect sizes are nested. Random effects are specified in the model formula with the same general style as mixed-model formulas. In the school data, effects are nested within studies, and studies are nested within districts.

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

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: 
#>      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

Fixed moderators can be combined with the random-effects structure in the same formula.

fit_ml_mod <- mars(
  data = school,
  formula = effect ~ 1 + year + (1 | district/study),
  studyID = "district",
  variance = "var",
  varcov_type = "multilevel",
  structure = "multilevel",
  estimation_method = "MLE"
)

summary(fit_ml_mod)
#> Results generated with MARS:v 0.5.2 
#> Friday, May 15, 2026
#> 
#> Model Type: 
#> multilevel
#> 
#> Estimation Method: 
#> Maximum Likelihood
#> 
#> Model Formula: 
#> effect ~ 1 + year + (1 | district/study)
#> 
#> Data Summary: 
#> Number of Effect Sizes: 56
#> Number of Fixed Effects: 2
#> Number of Random Effects: 2
#> 
#> Random Data Summary: 
#> Number unique district: 11 
#> Number unique study: 56
#> 
#> Random Components: 
#>      term     var     SD
#>  district 0.05767 0.2402
#>     study 0.03287 0.1813
#> 
#> Fixed Effects Estimates: 
#>    attribute  estimate        SE  z_test p_value     lower    upper
#>  (Intercept) -9.956118 17.092194 -0.5825  0.5602 -43.45620 23.54397
#>         year  0.005094  0.008585  0.5933  0.5530  -0.01173  0.02192
#> 
#> Model Fit Statistics: 
#>  logLik   Dev   AIC   BIC  AICc
#>  -8.389 28.78 24.78 32.88 27.84
#> 
#> Q Error: 550.26 (54), p < 0.0001
#> 
#> I2 (General): 
#>     names values
#>  district  92.34
#>     study  87.30
#> 
#> I2 (Jackson): 
#>     names values
#>  district  98.01
#>     study  98.01
#> 
#> I2 (Between): 82.55685
#> 
#> Residual Diagnostics: 
#>   n n_finite_raw mean_raw sd_raw   rmse    mae q_pearson mean_abs_studentized
#>  56           56 -0.05808 0.3193 0.3217 0.2568     94.61                1.055
#>  max_abs_studentized prop_abs_studentized_gt2 prop_abs_studentized_gt3
#>                4.094                   0.1071                  0.03571
#> 
#> Normality (whitened residuals):                   test n_tested statistic p_value
#>  shapiro_wilk_whitened       56    0.9577 0.04771
#> 
#> Heteroscedasticity trend (|raw residual| ~ fitted):   n corr_abs_raw_fitted  slope p_value
#>  56              0.1208 0.4562  0.3752

Choosing a Structure

The same entry point, mars(), handles each model family. The main choice is how the effect sizes are organized:

structure = "univariate" for one effect size per independent study unit,
effectsize_type = "cor" or "smd" for dependent multivariate effect sizes,
structure = "multilevel" for nested effect sizes described with random effects in the formula.

After fitting, summary() provides the main estimates and standard errors for each model.

Brandon LeBeau

2026-05-15

Univariate Meta-Analysis

Multivariate Meta-Analysis

Multilevel Meta-Analysis

Choosing a Structure