These analyses demonstrate with simple examples the correctness of the bias formulas given in Savalai (2019). Based on true score effects for a simple mediation model and reliability estimates, a SEM is fitted to the attenuated implied variance- covariance matrix.The difference in the estimated effects to the data-generating true effects give the same results as the analytic bias formulas.
# Clear workspace
rm(list = ls())
# Load packages
suppressPackageStartupMessages({
library("lavaan")
})
#
# Analytic bias for regressions weights
#
# Note: Assumes variances of 1 for the true scores.
#
# @param rxx A reliability estimate for X
# @param rmm A reliability estimate for M
# @param MONX A standardized regression weight of M on X
# @param YONM A standardized regression weight of Y on M
# @param YONX A standardized regression weight of Y on X
# @return A numeric vector
# @source https://doi.org/10.1037/met0000181
#
<- function(rxx = 1, rmm = 1, ryy = 1,
bias_a
MONX, YONM, YONX) {
# eq A6
<- rxx * MONX
a
# eq A7a
<- rmm * (YONM + MONX * YONX - rxx * (MONX * YONM + YONX) * MONX) /
b 1 - rxx * rmm * MONX^2)
(
# eq A8
<- rxx * ((YONX + MONX * YONM) - MONX * rmm * (YONM + MONX * YONX)) /
c 1 - rxx * rmm * MONX^2)
(
# biases
<- c(bias_a = MONX - a, bias_b = YONM - b, bias_c = YONX - c)
out
return(out)
}
#
# Generate attenuated variance-covariance matrix from
# standardized regression weights of a mediation model
#
# @param a A standardized regression weight of M on X
# @param b A standardized regression weight of Y on X
# @param c A standardized regression weight of Y on M
# @param rxx Reliability of X
# @param rmm Reliability of M
# @param ryy Reliability of Y
# @return A 3x3 variance-covariance matrix
#
<- function(a, b, c,
get_r rxx = 1, rmm = 1, ryy = 1) {
<- diag(1 / c(rxx, rmm, ryy))
m lower.tri(m)] <- m[upper.tri(m)] <- c(a, a * b + c, b + a * c)
m[colnames(m) <- rownames(m) <- c("X", "M", "Y")
return(m)
}
Unreliability in Y does not lead to bias.
# Attenuated effects
<- coef(sem(model = "M ~ X; Y ~ M + X;",
(b sample.cov = get_r(a = .3, b = .3, c = .3,
rxx = 1, rmm = 1, ryy = .70),
sample.nobs = 1000)))
## M~X Y~M Y~X M~~M Y~~Y
## 0.300 0.300 0.300 0.909 1.193
# Analytic bias
bias_a(ryy = .70, MONX = .30, YONM = .30, YONX = .30)
## bias_a bias_b bias_c
## 0 0 0
# True effect minus attenuated effect corresponds to analytic bias
30 - b["M~X"] .
## M~X
## 0
30 - b["Y~M"] .
## Y~M
## -5.551115e-17
30 - b["Y~X"] .
## Y~X
## 0
Unreliability in M does not bias a, but leads to biases in b and c.
# Attenuated effects
<- coef(sem(model = "M ~ X; Y ~ M + X;",
(b sample.cov = get_r(a = .3, b = .3, c = .3,
rxx = 1, rmm = .70, ryy = 1),
sample.nobs = 1000)))
## M~X Y~M Y~X M~~M Y~~Y
## 0.300 0.204 0.329 1.337 0.791
# Analytic bias
bias_a(rmm = .70, MONX = .30, YONM = .30, YONX = .30)
## bias_a bias_b bias_c
## 0.00000000 0.09605123 -0.02881537
# True effect minus attenuated effect corresponds to analytic bias
30 - b["M~X"] .
## M~X
## 0
30 - b["Y~M"] .
## Y~M
## 0.09605123
30 - b["Y~X"] .
## Y~X
## -0.02881537
Unreliability in X leads to biases in all effects.
# Attenuated effects
<- coef(sem(model = "M ~ X; Y ~ M + X;",
(b sample.cov = get_r(a = .3, b = .3, c = .3,
rxx = .70, rmm = 1, ryy = 1),
sample.nobs = 1000)))
## M~X Y~M Y~X M~~M Y~~Y
## 0.210 0.329 0.204 0.936 0.791
# Analytic bias
bias_a(rxx = .70, MONX = .30, YONM = .30, YONX = .30)
## bias_a bias_b bias_c
## 0.09000000 -0.02881537 0.09605123
# True effect minus attenuated effect corresponds to analytic bias
30 - b["M~X"] .
## M~X
## 0.09
30 - b["Y~M"] .
## Y~M
## -0.02881537
30 - b["Y~X"] .
## Y~X
## 0.09605123