Fit a Multilevel Model Specification

Centers requested predictors, builds an lme4 formula, and fits the specified mixed-effects or multilevel model.

Usage

mlm_fit(spec, REML = TRUE, optimizer = "bobyqa", maxfun = 20000)

# S3 method for class 'mlm_result'
print(x, ...)

Arguments

spec: An "mlm_spec" object created with mlm_spec.
REML: Logical; whether to use restricted maximum likelihood for Gaussian models. Use REML = FALSE when comparing nested models that differ in fixed effects.
optimizer: Optimizer passed to lme4, such as "bobyqa", "Nelder_Mead", or "nloptwrap".
maxfun: Maximum number of function evaluations.
x: An "mlm_result" object returned by mlm_fit().
...: Ignored.

Value

An object of class "mlm_result": a list with the fitted lme4 model (fit), the model formula (formula), the data including centered predictors (data), R code that recreates the centered predictors (centering_code), the model specification (spec), and the estimation settings (settings: REML, optimizer, and maxfun). Printing the object shows a compact summary with the fixed-effects table.

Details

For Gaussian models, mlm_fit() calls lme4::lmer(). For supported generalized mixed models, it calls lme4::glmer() or lme4::glmer.nb(). The returned fit element is the underlying lme4 model object, so users can pass it to standard lme4 methods.

The function performs requested grand-mean and cluster-mean centering before fitting. Centering uses the analysis sample: means are computed from the rows that are complete on every model variable, which are the rows lme4 keeps. A grand-mean-centered predictor therefore has a mean of zero, and a cluster-mean-centered predictor has cluster means of zero, among the observations that enter the model. Rows with missing values stay in the returned data but are not used for estimation.

mlm_fit() stops with an explanatory message when a specification cannot be fit as requested, for example when a categorical predictor is centered, a random slope is requested for a variable that is not a fixed predictor, or no random effects are specified. Users should inspect convergence messages, singular fits, residual patterns, and whether the random-effects structure is supported by the data.

Examples

dat <- example_hsb(n_schools = 4, min_students = 5, max_students = 6)
spec <- mlm_spec(
  outcome = "mathscore",
  fixed = list(ses = list(center = "CWC")),
  grouping = list(schoolid = "schoolid"),
  random = list(schoolid = list(intercept = TRUE, slopes = character(),
    correlation = TRUE)),
  data = dat
)
result <- mlm_fit(spec)
result
#> <mlm_result>
#>   Formula:      mathscore ~ ses_CWC + (1 | schoolid)
#>   Distribution: gaussian (identity link), estimated by REML
#>   Observations: 22 of 22 rows
#>   Groups:       schoolid = 4
#>   AIC:          163.63
#>   Singular fit: no
#> 
#>    Predictor     b   SE Statistic      p         95% CI
#>  (Intercept) 72.78 3.46     21.06 < .001 [66.01, 79.55]
#>      ses_CWC  3.22 3.77      0.85  0.393 [-4.17, 10.62]
#> 
#> Use mlm_apa_tables(), mlm_latex_equations(), or mlm_publication_bundle() for reporting output.
summary(result$fit)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: mathscore ~ ses_CWC + (1 | schoolid)
#>    Data: spec$data
#> Control: control
#> 
#> REML criterion at convergence: 155.6
#> 
#> Scaled residuals: 
#>      Min       1Q   Median       3Q      Max 
#> -1.33488 -0.59912 -0.08057  0.34960  2.33011 
#> 
#> Random effects:
#>  Groups   Name        Variance Std.Dev.
#>  schoolid (Intercept) 30.61    5.532   
#>  Residual             93.78    9.684   
#> Number of obs: 22, groups:  schoolid, 4
#> 
#> Fixed effects:
#>             Estimate Std. Error t value
#> (Intercept)   72.780      3.455  21.065
#> ses_CWC        3.224      3.773   0.854
#> 
#> Correlation of Fixed Effects:
#>         (Intr)
#> ses_CWC 0.000