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)
  
}

1 rxx = 1, rmm = 1, ryy = .70

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 
##   0

2 rxx = 1, rmm = .70, ryy = 1

Unreliability 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.02881537

3 rxx = .70, rmm = 1, ryy = 1

Unreliability 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