MASMR-Workflow.Rmdmasmr() implements a two-stage workflow:
library(mars)
teacher_expectancy
#> study author year weeks setting tester n1i n2i yi vi
#> 1 1 Rosenthal et al. 1974 2 group aware 77 339 0.03 0.0156
#> 2 2 Conn et al. 1968 21 group aware 60 198 0.12 0.0216
#> 3 3 Jose & Cody 1971 19 group aware 72 72 -0.14 0.0279
#> 4 4 Pellegrini & Hicks 1972 0 group aware 11 22 1.18 0.1391
#> 5 5 Pellegrini & Hicks 1972 0 group blind 11 22 0.26 0.1362
#> 6 6 Evans & Rosenthal 1969 3 group aware 129 348 -0.06 0.0106
#> 7 7 Fielder et al. 1971 17 group blind 110 636 -0.02 0.0106
#> 8 8 Claiborn 1969 24 group aware 26 99 -0.32 0.0484
#> 9 9 Kester 1969 0 group aware 75 74 0.27 0.0269
#> 10 10 Maxwell 1970 1 indiv blind 32 32 0.80 0.0630
#> 11 11 Carter 1970 0 group blind 22 22 0.54 0.0912
#> 12 12 Flowers 1966 0 group blind 43 38 0.18 0.0497
#> 13 13 Keshock 1970 1 indiv blind 24 24 -0.02 0.0835
#> 14 14 Henrikson 1970 2 indiv blind 19 32 0.23 0.0841
#> 15 15 Fine 1972 17 group aware 80 79 -0.18 0.0253
#> 16 16 Grieger 1970 5 group blind 72 72 -0.06 0.0279
#> 17 17 Rosenthal & Jacobson 1968 1 group aware 65 255 0.30 0.0193
#> 18 18 Fleming & Anttonen 1971 2 group blind 233 224 0.07 0.0088
#> 19 19 Ginsburg 1970 7 group aware 65 67 -0.07 0.0303
fit_masmr <- masmr(
data = teacher_expectancy,
methodology_mods = c("factor(tester)"),
substantive_mods = c("year", "weeks", "factor(setting)"),
formula = yi ~ 1 + year + weeks + factor(setting) + factor(tester),
studyID = "study",
variance = "vi",
varcov_type = "univariate",
structure = "univariate",
estimation_method = "REML"
)
summary(fit_masmr)
#> Methodology-Adjusted Substantive Meta-Regression (MASMR)
#> Workflow: formula
#> Methodology moderators: factor(tester)
#> Substantive moderators: year, weeks, factor(setting)
#>
#> Stage 1 (Methodology Adjustment)
#> Results generated with MARS:v 0.5.1.1
#> Thursday, April 2, 2026
#>
#> Model Type:
#> univariate
#>
#> Estimation Method:
#> Restricted Maximum Likelihood
#>
#> Model Formula:
#> yi ~ 1 + factor(tester)
#> <environment: 0x55b241848e58>
#>
#> Data Summary:
#> Number of Effect Sizes: 19
#> Number of Fixed Effects: 2
#> Number of Random Effects: 1
#>
#> Random Components:
#> term var SD
#> intercept 0.02476 0.1574
#>
#> Fixed Effects Estimates:
#> attribute estimate SE z_test p_value lower upper
#> (Intercept) 0.04635 0.07187 0.6448 0.5190 -0.09452 0.1872
#> blind 0.10248 0.11209 0.9142 0.3606 -0.11722 0.3222
#>
#> Model Fit Statistics:
#> logLik Dev AIC BIC AICc
#> -3.838 7.675 13.68 16.17 23.86
#>
#> Q Error: 35.157 (17), p = 0.0059
#>
#> I2 (General): 49.037
#>
#> Residual Diagnostics:
#> n n_finite_raw mean_raw sd_raw rmse mae q_pearson mean_abs_studentized
#> 19 19 0.06879 0.3547 0.352 0.2418 22.86 0.907
#> max_abs_studentized prop_abs_studentized_gt2 prop_abs_studentized_gt3
#> 2.846 0.1053 0
#>
#> Normality (whitened residuals): test n_tested statistic p_value
#> shapiro_wilk_whitened 19 0.9141 0.08801
#>
#> Heteroscedasticity trend (|raw residual| ~ fitted): n corr_abs_raw_fitted slope p_value
#> 19 -0.1174 -0.5868 0.6322
#>
#> Stage 2 (Substantive Model on Adjusted Effects)
#> Results generated with MARS:v 0.5.1.1
#> Thursday, April 2, 2026
#>
#> Model Type:
#> univariate
#>
#> Estimation Method:
#> Restricted Maximum Likelihood
#>
#> Model Formula:
#> yi_masmr ~ 1 + year + weeks + factor(setting)
#> <environment: 0x55b24184aa68>
#>
#> Data Summary:
#> Number of Effect Sizes: 19
#> Number of Fixed Effects: 4
#> Number of Random Effects: 1
#>
#> Random Components:
#> term var SD
#> intercept 0.01345 0.116
#>
#> Fixed Effects Estimates:
#> attribute estimate SE z_test p_value lower upper
#> (Intercept) 37.12140 53.047587 0.6998 0.48407 -66.84996 141.092756
#> year -0.01872 0.026926 -0.6952 0.48693 -0.07149 0.034055
#> weeks -0.01349 0.006082 -2.2184 0.02653 -0.02541 -0.001572
#> indiv 0.16226 0.184128 0.8812 0.37819 -0.19862 0.523146
#>
#> Model Fit Statistics:
#> logLik Dev AIC BIC AICc
#> -1.5 3.001 13 16.54 50.71
#>
#> Q Error: 25.48 (15), p = 0.0439
#>
#> I2 (General): 35.03433
#>
#> Residual Diagnostics:
#> n n_finite_raw mean_raw sd_raw rmse mae q_pearson mean_abs_studentized
#> 19 19 0.05031 0.3186 0.3142 0.2032 19.35 0.8655
#> max_abs_studentized prop_abs_studentized_gt2 prop_abs_studentized_gt3
#> 2.861 0.05263 0
#>
#> Normality (whitened residuals): test n_tested statistic p_value
#> shapiro_wilk_whitened 19 0.9312 0.1821
#>
#> Heteroscedasticity trend (|raw residual| ~ fitted): n corr_abs_raw_fitted slope p_value
#> 19 0.3177 0.5347 0.1851The stage-2 model (fit_masmr$stage2) represents
relationships with substantive moderators after methodology-associated
variation is removed in stage 1.
The table below aligns coefficient estimates from both stages for quick inspection.
extract_coef_table <- function(fit) {
b <- if(!is.null(fit$beta_r)) fit$beta_r else fit$beta_fe
data.frame(
term = colnames(fit$design_matrix),
estimate = as.numeric(b),
stringsAsFactors = FALSE
)
}
coef_stage1 <- extract_coef_table(fit_masmr$stage1)
coef_stage2 <- extract_coef_table(fit_masmr$stage2)
coef_compare <- merge(
coef_stage1, coef_stage2,
by = "term", all = TRUE,
suffixes = c("_stage1_methodology", "_stage2_substantive")
)
coef_compare$present_in <- ifelse(
is.na(coef_compare$estimate_stage1_methodology), "stage2_only",
ifelse(is.na(coef_compare$estimate_stage2_substantive), "stage1_only", "both")
)
coef_compare
#> term estimate_stage1_methodology estimate_stage2_substantive
#> 1 (Intercept) 0.04634783 37.1213966
#> 2 factor(setting)indiv NA 0.1622616
#> 3 factor(tester)blind 0.10247924 NA
#> 4 weeks NA -0.0134931
#> 5 year NA -0.0187189
#> present_in
#> 1 both
#> 2 stage2_only
#> 3 stage1_only
#> 4 stage2_only
#> 5 stage2_onlyNA values are expected here: MASMR intentionally uses
different predictor blocks across stages. A term appears as
NA in the stage where it is not included by design.
If you only want directly comparable rows, filter to shared terms:
subset(coef_compare, present_in == "both")
#> term estimate_stage1_methodology estimate_stage2_substantive
#> 1 (Intercept) 0.04634783 37.1214
#> present_in
#> 1 bothFor multilevel models, keep the random-effects structure in
formula and call masmr() the same way. Random
terms are preserved across both stages.
fit_masmr_ml <- masmr(
data = your_data,
methodology_mods = c("risk_bias", "blinded"),
substantive_mods = c("dose", "duration"),
formula = yi ~ 1 + risk_bias + blinded + dose + duration + (1 | study/effect),
studyID = "study",
variance = "vi",
varcov_type = "univariate",
structure = "multilevel",
estimation_method = "REML"
)When using effectsize_type = "smd" with multivariate
covariates:
fit_masmr_mv <- masmr(
data = your_data,
methodology_mods = c("risk_bias", "blinded"),
substantive_mods = c("dose", "duration"),
studyID = "study",
effectID = "effect",
sample_size = "n",
effectsize_type = "smd",
effectsize_name = "yi",
variance = "vi",
varcov_type = "univariate",
structure = "UN",
estimation_method = "REML"
)In this mode, moderator blocks are passed through
multivariate_covs internally for stage 1 and stage 2.