cid <int> | y <int> | x1 <dbl> | x2 <int> |
|---|---|---|---|
| 1 | 10 | 0 | 2 |
| 1 | 5 | 0 | 0 |
| 1 | 5 | 0 | 1 |
| 1 | 7 | 0 | 2 |
| 2 | 5 | 1 | 2 |
| 2 | 5 | 1 | 2 |
| 2 | 5 | 1 | 2 |
| 2 | 8 | 1 | 1 |
| 3 | 8 | 1 | 0 |
| 3 | 5 | 1 | 0 |
Unlocking the Multilevel Equation
Multicollinearity, Effect Size, Design Effects, and Power Analysis
Mark Lai
University of Southern California
2024-03-22
Overview
- The mixed model equation
- Multicollinearity and variance partitioning (with Venn diagrams)
- Standardized effect size
- Design effect/variance inflation
- Power analysis
Mixed Model Equation
Multilevel Data
Mixed Model Equations
y=Xγ+Zu+e
[105575558]=[102100101102112112112111][γ0γ1γ2]+[1210111211121212][u01u02u11u12]+[e11e12e13e14e21e22e23e24]
[105575558]=[102100101102112112112111][γ0γ1γ2]+[1210111211121212][u01u02u11u12]+[e11e12e13e14e21e22e23e24]
For cluster j, yj=Xjγ+Zjuj+ej
Mixed Model Equations (cont’d)
y=Xγ+Zu+e
- Z = [Z0⋯Zq−1], and each component is block-diagonal
- e is a vector of errors;1 V(e)=σ2I
- u: vector of length qJ of random effects; V(uj)=T=σ2D∀j
D=[τ20/σ2τ01/σ2τ21/σ2⋮⋯⋱τ0q−1/σ2τ1q−1/σ2⋯τ2q−1/σ2] - u and e are mutually independent, and both are indepedent of X and Z
Multicollinearity
Variance Partitioning
In regression, we look at
- Xc⊤Xc/N, covariance among columns of Xc (assume mean centered predictors)
If X1 and X2 are uncorrelated, xc1⊤xc2=0, and variance accounted for by the two predictors are separate
If X1 and X2 are correlated, xc1⊤xc2≠0, then adding X2 . . .
- changes the coefficient of X1, and
- increases the standard error of the coefficient
Rearranging the MLM Equation
[105575558]=[102100101102112112112111][γ0γ1γ2]+[1210111211121212][u01u02u11u12]+[e11e12e13e14e21e22e23e24]
[105575558]=[1021210010101111021211211112121121211112][γ0γ1γ2u01u02u11u12]+[e11e12e13e14e21e22e23e24]
In MLM, we also look at
Z⊤Z, cross-products among columns of Z, and
Xc⊤Z, cross-products between columns of Xc and those of Z
Example 1: Random Intercepts Only
- Residual maker matrix: MZ0 = I−Z0(Z⊤0Z0)−1Z⊤0
- MZ0Z0 = 0
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0.67 -0.33 -0.33 0.00 0.00 0.00
[2,] -0.33 0.67 -0.33 0.00 0.00 0.00
[3,] -0.33 -0.33 0.67 0.00 0.00 0.00
[4,] 0.00 0.00 0.00 0.67 -0.33 -0.33
[5,] 0.00 0.00 0.00 -0.33 0.67 -0.33
[6,] 0.00 0.00 0.00 -0.33 -0.33 0.67Example 2a: Level-2 Predictor
[1] -1 -1 -1 1 1 1
Level-2 predictor can only account for level-2 variance
Example 2b: Pure Level-1 Predictor
[1] -2 0 2 -3 2 1
Pure level-1 predictor can only account for level-1 variance
Example 2c: General Level-1 Predictor
[1] -1 1 2 -3 0 1
Level-1 predictor can account for both level-2 and level-1 variance
Example 3: Random slopes
- Without cluster-mean centering, Z1=diag(x11,x12,…,x1J)
[,1] [,2]
[1,] -1 0
[2,] 1 0
[3,] 2 0
[4,] 0 -3
[5,] 0 0
[6,] 0 1- For general level-1 predictor
- Random slope variance (τ21), if omitted, will be redistributed to the fixed effect, random intercept variance (τ20), and σ2
- For purely level-1 predictor, Z⊤1Z0 = 0
- Random slope variance (τ21), if omitted, will be redistributed to only fixed effect and σ2
Revisiting Lai & Kwok (2015)
Using simulations, we found that, when clustering is ignored, SEs of the following are underestimated
- Lv-2 predictor
- General Lv-1 predictor (i.e., with unequal cluster means)
- Lv-1 predictor with random slopes
Because they all correlate with Z!
Additional Questions
- Can we use this framework to quantify model misspecification?
- How does omitting components of X affect estimates of γ and τ?
- How does omitting components of Z affect estimates of γ and τ?
- How do we generalize this beyond two-level models?
- E.g., Crossed random effects (e.g., Lai, 2019)1
Effect Size
Total Variance of y
Conditional on X
V(Zu+e)=σ2V=σ2[Z(D⊗IJ)Z⊤+I],
Unconditional (Xc = grand-mean centered X)
Xcγγ⊤Xc⊤⏟fixed+σ2Z(D⊗IJ)Z⊤⏟random lv-2+σ2I⏟random lv-1.
Average Variance of Each Observation
- Random slopes imply nonconstant variance across observations1
- V(xijuj) depends on xij
Average Variance of Each Observation (cont’d)
Tr[V(y)]/N=γ⊤Xc⊤Xcγ/N+σ2Tr[Z(D⊗IJ)Z⊤]/N+σ2
Using Summary Statistics
Tr[V(y)]/N=γ⊤ΣXγ+σ2[Tr(DKZ)+1],
Let Xr be the design matrix of variables with random effects1
- ΣX = Xc⊤Xc/N (covariance matrix of X)
- NKZ = Xr⊤Xr (uncentered crossproduct matrix)
- KZ=ˉXr⊤ˉXr+ΣrX
Intraclass Correlations
Unconditional or Conditional on the fixed effects
- i.e., proportion of variance at level 2 out of all variance unaccounted for by the fixed effects
ICCD=Tr(DKZ)Tr(DKZ)+1,
Multilevel R2
Proportion of variance by fixed effects (see Johnson, 2014; Rights & Sterba, 2019)1
R2GLMM=γ⊤ΣXγγ⊤ΣXγ+σ2[Tr(DKZ)+1].
Code Example
Linear mixed model fit by REML ['lmerMod']
Formula: mAch ~ ses + meanses + (ses | school)
Data: Hsb82
REML criterion at convergence: 46561.4
Scaled residuals:
Min 1Q Median 3Q Max
-3.1671 -0.7269 0.0163 0.7547 2.9646
Random effects:
Groups Name Variance Std.Dev. Corr
school (Intercept) 2.695 1.6418
ses 0.453 0.6731 -0.21
Residual 36.796 6.0659
Number of obs: 7185, groups: school, 160
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.6740 0.1506 84.176
ses 2.1903 0.1218 17.976
meanses 3.7812 0.3826 9.882 [,1]
[1,] 8.190918Code Example (cont’d)
library(r2mlm)
library(MuMIn)
# R^2 using formula, r2mlm::r2mlm(), and MuMIn::r.squaredGLMM()
list(formula = vf / (vf + vr + sigma(m1)^2),
r2mlm = r2mlm(m1),
MuMIn = MuMIn::r.squaredGLMM(m1))$formula
[,1]
[1,] 0.170797
$r2mlm
$r2mlm$Decompositions
total
fixed 0.170816580
slope variation 0.005737776
mean variation 0.056202223
sigma2 0.767243421
$r2mlm$R2s
total
f 0.170816580
v 0.005737776
m 0.056202223
fv 0.176554356
fvm 0.232756579
$MuMIn
R2m R2c
[1,] 0.1708167 0.232756
Confidence interval of R2
- Delta method
- Bootstrapping (with
lme4::bootMer()or R package bootmlm)1 - Bayesian estimation
Additional Questions
- What is the bias of sample R2 estimator?
- Can we compute adjusted R2?
Standardized Mean Difference
SMD = Mean difference / SDp1
For cluster-randomized trials
Yij=Tjγ+Z0ju0j+eij,V(u0j)=τ20
Tj = 0 (control) and 1 (treatment)
δ=γ√σ2[Tr(DKZ)+1]=γ√τ20+σ2
More general form:
δ=adjusted/unadjusted treatment effect√average conditional variance of y
δ=E[E(Y|T=1)−E(Y|T=0)|C]√Tr[V(Xγ|T)+V(Zu|T)]/N+σ2
where C is the subset of covariates for adjusted differences
Linking R2 to δ
f2GLMM=R2GLMM1−R2GLMM
- For two-group designs with one predictor, 2f = d
Additional Questions
- What ranges of effect sizes do published MLM studies have?
Standardized Coefficients
γs=γSDx√Tr[V(y)]/N
Linear mixed model fit by REML ['lmerMod']
Formula: mAch ~ ses + meanses + (ses | school)
Data: Hsb82
REML criterion at convergence: 46561.4
Scaled residuals:
Min 1Q Median 3Q Max
-3.1671 -0.7269 0.0163 0.7547 2.9646
Random effects:
Groups Name Variance Std.Dev. Corr
school (Intercept) 2.695 1.6418
ses 0.453 0.6731 -0.21
Residual 36.796 6.0659
Number of obs: 7185, groups: school, 160
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.6740 0.1506 84.176
ses 2.1903 0.1218 17.976
meanses 3.7812 0.3826 9.882
Correlation of Fixed Effects:
(Intr) ses
ses -0.084
meanses -0.013 -0.256For cluster means and cluster-mean centered variables, it is still natural to use the total SD
- SDSES = 0.78
- ΣX = [0.610.170.170.17]
- KZ=[1000.61], σ2D = [2.7−0.23−0.230.45]
- Tr[V(y)]/N = γ⊤ΣXγ + σ2[Tr(DKZ)+1] = 8.19 + 2.97 + 36.8 = 47.96
- γsses = 2.19 / √47.96 = 0.25; γsmeanses = 3.78 / √47.96 = 0.43
Design Effect
Design Effect/Variance Inflation
Design effect: Expected impact of design on sampling variance of an estimator1
For the sample mean (ˆμ):
- Simple random sample: VSRS(ˆμSRS) = ˜σ2/N
- Assuming constant total variance: ˜σ2=τ20+σ2
- With clustering:2 VMLM(ˆγ0) = τ20/J+σ2/N
Deff(ˆμ)=VSRS(ˆγ0)VMLM(ˆμSRS)=1+(n−1)ICC
Deff(ˆμ) Does Not Inform About Random Slopes
- Even when Deff(ˆμ) is close to 1, random slope variance τ21 can still be large
- Random slopes can be independent from random intercepts
Deff for ˆγ
ˆγ=(X⊤V−1X)−1X⊤V−1y
VMLM(ˆγGLS)=σ2(X⊤V−1X)−1VSRS(ˆγOLS)=˜σ2(X⊤X)−1
where ˜σ2=σ2[Tr(DKZ)+1]
Unpacking V(ˆγGLS)
σ2(X⊤V−1X)−1=σ2{X⊤[Z(D⊗IJ)Z⊤+I]−1X}−1=σ2{X⊤[I−Z(D⊗I−1J+Z⊤Z)Z⊤]X}−1=σ2{X⊤X−X⊤Z(D⊗I−1J+Z⊤Z)Z⊤X}−1=σ2{[(X⊤X)−1+(X⊤X)−1X⊤Z(D⊗IJ)Z⊤X(X⊤X)−1]−1}−1=σ2[(X⊤X)−1+(X⊤X)−1X⊤Z(D⊗IJ)Z⊤X(X⊤X)−1]
We can substitute X with Xc. As Z is block diagonal, Xc⊤Z = [Xc⊤1Z1|Xc⊤2Z2|⋯|Xc⊤JZJ], and
NG=Xc⊤Z(D⊗IJ)Z⊤Xc=J∑j=1Xc⊤jZjDZ⊤jXcj
is the cross-product matrix of Xc⊤Zu, and Gkl seems equal to
{˜nNΣXkZDΣ⊤XlZwhen both Xl and Xk are purely level-1˜nN[ΣXkZDΣ⊤XlZ+KXkZDK⊤XlZ]otherwise
where ˜n=∑Jj=1n2j/N (reduced to n if cluster sizes are equal across clusters). So
σ2(Xc⊤V−1Xc)−=σ2[Σ−1X+˜nΣ−1XGΣ−1X],
where G = (ΣXZDΣ⊤XZ+FKXZDK⊤XZF), F is a diagonal filtering matrix in the form diag[f0,f1,…,fp−1], where fk = 1 if Xk has a level-2 component and 0 otherwise.
The OLS standard error and the MLM standard error have different rates of converagence to 0
For coefficient γk
Deff(ˆγk)={Σ−1X}kk+n[Σ−1XGΣ−1X]kk{Σ−1X}kk[Tr(DKZ)+1]=(1−ICCD)+(1−ICCD)n[Σ−1XGΣ−1X]kkΣ−1Xkk
- Generally speaking, deff is larger with a larger cluster size, n
- Deff is fixed for constant n (as in Deff[ˆμ])
Example: Growth curve modeling with Tx × slope interaction
yi=[1Ti001Ti1Ti1Ti22Ti1Ti33Ti1Ti44Ti][γ0γ1γ2γ3]+[1011121314][u0u1]+[ei1ei2ei3ei4ei5]
γ=[1,0,−0.1,−0.3]⊤,D=[0.50.20.1],σ2=1.5
icc_d r2_GLMM
0.655 0.056 Example cont’d
# Design effect
sigmaxz <- sigmax[, c(1, 3)]
Kxz <- Kx[, c(1, 3)]
F_mat <- diag(c(0, 1, 0, 1)) # filtering
G_mat <- (sigmaxz %*% D_mat %*% t(sigmaxz)) +
F_mat %*% Kxz %*% D_mat %*% t(Kxz) %*% F_mat
sigmax_inv <- MASS::ginv(sigmax)
(1 - icc) * (1 + num_obs * diag(sigmax_inv %*% G_mat %*% sigmax_inv) /
diag(sigmax_inv))[1] NaN 0.6321839 0.6896552 0.6896552Power Analysis
For H0: γk=0
- Wald statistic: z=ˆγk/√ˆV(ˆγk)
- In large sample, z∼N(z1,1), z1 = γk/√V(ˆγk)
- With significance level α, power is approximately Φ(z1−α/2−z1)+Φ(zα/2−z1)
- Growth model example (Tx × time interaction)
- Line: analytic formula; Points: 1,000 simulation samples
Summary
- Multicollinearity can happen between random-effect predictors, and between random- and fixed-effect predictors
- With random coefficients, one can do variance partitioning using the average variance of an observation
- ICCD quantifies the proportion of variance due to random effects
- R2 quantifies the proportion of variance due to fixed effects
Summary (cont’d)
- Standardized coefficients and standardized effect size can be defined
- Design effect quantifies variance inflation of an estimator due to cluster sampling
- Can be applied to fixed-effect estimators
- Closed-form expression for power can be obtained using summary statistics
Thank You!
If you have any suggestions or feedback, please email me at hokchiol@usc.edu
References
Aguinis, H., & Culpepper, S. A. (2015). An expanded decision-making procedure for examining cross-level interaction effects With multilevel modeling. Organizational Research Methods, 18(2), 155–176. https://doi.org/10.1177/1094428114563618
Hedges, L. V. (2009). Effect sizes in nested designs. In H. Cooper, L. V. Hedges, & J. C. Valentine, The handbook of research synthesis and meta-analysis (2nd ed., pp. 337–355). Russell Sage Foundation.
Johnson, P. C. D. (2014). Extension of Nakagawa & Schielzeth’s R 2 GLMM to random slopes models. Methods in Ecology and Evolution, 5(9), 944–946. https://doi.org/10.1111/2041-210X.12225
Lai, M. H. C. (2019). Correcting fixed effect standard errors when a crossed random effect was ignored for balanced and unbalanced designs. Journal of Educational and Behavioral Statistics, 44(4), 448–472. https://doi.org/10.3102/1076998619843168
Lai, M. H. C., & Kwok, O. (2015). Examining the rule of thumb of not using multilevel modeling: The “design effect smaller than two” rule. The Journal of Experimental Education, 83(3), 423–438. https://doi.org/10.1080/00220973.2014.907229
Luo, W., & Kwok, O. (2009). The impacts of ignoring a crossed factor in analyzing cross-classified data. Multivariate Behavioral Research, 44(2), 182–212. https://doi.org/10.1080/00273170902794214
Murayama, K., Usami, S., & Sakaki, M. (2022). Summary-statistics-based power analysis: A new and practical method to determine sample size for mixed-effects modeling. Psychological Methods. https://doi.org/10.1037/met0000330
Rights, J. D., & Sterba, S. K. (2019). Quantifying explained variance in multilevel models: An integrative framework for defining R-squared measures. Psychological Methods, 24(3), 309–338. https://doi.org/10.1037/met0000184
Rights, J. D., & Sterba, S. K. (2023). R-squared measures for multilevel models with three or more levels. Multivariate Behavioral Research, 58(2), 340–367. https://doi.org/10.1080/00273171.2021.1985948
Snijders, T. A. B., & Bosker, R. J. (1993). Standard errors and sample sizes for two-level research. Journal of Educational Statistics, 18(3), 237–259. https://doi.org/10.3102/10769986018003237
Supplemental Slides
Multicollinearity with Crossed Random Effects
y=Xγ+ZAτA+ZBτB+ϵ
- (MX0ZA)⊤MX0ZB = collinearity between the two crossed levels1
- When there are only random intercepts at Levels A and B,
- For balanced designs (i.e., fully crossed), ZA and ZB are orthogonal (i.e., MX0Z⊤AMX0ZB = 0)
- If Level B is omitted, its variance goes to Level 1
- For unbalanced designs (i.e., partially crossed)
- If Level B is omitted, some of its variance goes to Level A
- For balanced designs (i.e., fully crossed), ZA and ZB are orthogonal (i.e., MX0Z⊤AMX0ZB = 0)
3 x 10 sparse Matrix of class "dgCMatrix"
a 2 2 2 2 2 2 2 2 2 2
b 2 2 2 2 2 2 2 2 2 2
c 2 2 2 2 2 2 2 2 2 23 x 10 Matrix of class "dgeMatrix"
A B C D E F G H I J
a 0 0 0 0 0 0 0 0 0 0
b 0 0 0 0 0 0 0 0 0 0
c 0 0 0 0 0 0 0 0 0 0(Unconditional) ICC and R2
From a random-intercepts model without fixed-effect predictors,
ICC = maximum R2 for a level-2 predictor
1 - ICC = maximum R2 for a purely level-1 predictor
(Conditional) ICC and R2
From a conditional model with fixed-effect predictors,
ICC = Proportion of variance accounted for by level-2 random effects
R2 = Proportion of variance accounted for by all fixed-effect predictors
SMD With Covariates
With covariates X2, Yij=Tjγ1+x⊤2ijγ2+Z0ju0j+eij
- →T⊤X2=0 and no interactions
δ=γ1√γ⊤2ΣX2γ2+τ20+σ2
With random slopes, Yij=Tjγ1+x⊤2ijγ2+z⊤ijuj+eij
- and also →T⊤Z=0
δ=γ1√γ⊤2ΣX2γ2+σ2[Tr(DKZ)+1]
Code example (data)
Code example (cont’d)
library(lme4)
# Fit MLM
m1 <- lmer(y ~ treat + x2 + (x2 | clus_id), data = sim_data)
# Variance by x2
Sigma_x <- with(sim_data, { x2c <- x2 - mean(x2); mean(x2c^2) })
va_x2 <- Sigma_x * fixef(m1)[["x2"]]^2
# Variance by random intercepts and random slopes
# Using parameterization in lme4
Kz <- crossprod(cbind(1, sim_data$x2)) / nrow(sim_data)
va_z <- sum(VarCorr(m1)$clus_id * Kz) # Tr(Tau, Kz)
# Show all components
c(gamma_1 = fixef(m1)[["treat"]], va_x2 = va_x2,
va_z = va_z, va_e = sigma(m1)^2) gamma_1 va_x2 va_z va_e
0.257714138 0.001712517 0.537828861 0.989420368 [1] 0.2084203Additional Notes on SMD
- SE and CI: Delta method, bootstrap, Bayesian (also approximate variance as in Hedges, 2009)
- Using unadjusted (marginal) mean difference when →T⊤X2≠0
- Unadjusted mean difference:
E(Y|T=1,u) - E(Y|T=1,u) = γ1+Σ−1TTΣTX2γ21 - Unadjusted within-arm variance by fixed effects:
V(ˆY|T,u)=γ⊤ΣXγ−γ21ΣTT−2γ1ΣTX2γ⊤2−Σ−1TTγ⊤2ΣX2TΣTX2 - More complicated when T also correlates with Z
- Unadjusted mean difference:
Choice of Denominator in SMD
- Square root of total variance should be default1
- Longitudinal: between-person variance
- i.e., variance at time 0 (or time t?)
- Multisite trials: Exclude treatment effect heterogeneity (i.e., random slope of T)?
- i.e., total variance without the intervention
Power Analysis
Wait, This Is Not New . . .
What can still be done?
- Implementing formula-based approach in other software
- Express V(ˆγGLS) in terms of effect sizes
- Extend to models with other error covariance structure and more levels
- Handling uncertainty in nuisance parameters (e.g., ICC)1
Handling Nuisance Parameters in Power Analysis
- There are many nuisance parameters (e.g., ICC, correlation among predictors) that affect power
- Usually we just choose some convenient values
- A more disciplined approach: incorporate our uncertainty of those parameters as Bayesian priors
- See the HCB shiny application by Winnie Tse
Unlocking the Multilevel Equation Multicollinearity, Effect Size, Design Effects, and Power Analysis Mark Lai University of Southern California 2024-03-22