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
#
bias_a <- function(rxx = 1, rmm = 1, ryy = 1,
MONX, YONM, YONX) {
# eq A6
a <- rxx * MONX
# eq A7a
b <- rmm * (YONM + MONX * YONX - rxx * (MONX * YONM + YONX) * MONX) /
(1 - rxx * rmm * MONX^2)
# eq A8
c <- rxx * ((YONX + MONX * YONM) - MONX * rmm * (YONM + MONX * YONX)) /
(1 - rxx * rmm * MONX^2)
# biases
out <- c(bias_a = MONX - a, bias_b = YONM - b, bias_c = YONX - c)
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
#
get_r <- function(a, b, c,
rxx = 1, rmm = 1, ryy = 1) {
m <- diag(1 / c(rxx, rmm, ryy))
m[lower.tri(m)] <- m[upper.tri(m)] <- c(a, a * b + c, b + a * c)
colnames(m) <- rownames(m) <- c("X", "M", "Y")
return(m)
}Unreliability in Y does not lead to bias.
# Attenuated effects
(b <- coef(sem(model = "M ~ X; Y ~ M + X;",
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
## 0Unreliability in M does not bias a, but leads to biases in b and c.
# Attenuated effects
(b <- coef(sem(model = "M ~ X; Y ~ M + X;",
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.02881537Unreliability in X leads to biases in all effects.
# Attenuated effects
(b <- coef(sem(model = "M ~ X; Y ~ M + X;",
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