2
$\begingroup$

If we have two variables, $X$ and $M$, which follow a multivariate normal distribution, have a mean of 0, are correlated ($r_{XM}$), and have known reliabilities ($\rho_X$ and $\rho_M$), then the reliability of the interaction term $XM$ ($XM = X * M$) is:

$\rho_{XM} = \frac{\rho_X*\rho_M + r_{XM}^2 }{1 +r_{XM}^2} $

See link1 (Eq. 10), link2 (Eq. 19). Reliability can be thought of as a latent property of a measurement that is estimated. This tutorial is a good resource. It is generally defined as the fraction of an observed score variance that was not error:

Reliability of $X$ = $\rho_X = 1 - \frac{\sigma^2_\epsilon}{V_x}$

Where $\sigma^2_\epsilon$ is variance of the measurement error and $V_X$ is the observed variance. In practice, how reliability is estimated will vary depending on what kind of measurement we are interested. See this extensive tutorial that reviews some of the most common methods. Note that how reliability was estimated is not germane to my question, only that it exists and is known.

Now, my question is: If we add a third variable $Z$, following all the assumptions above, which is correlated with both $X$ and $M$, with a known reliability ($\rho_Z$), what is the reliability of the interaction term $XMZ$ ($XMZ = X*M*Z$)?

Edits:

As an example, we can simulate some data and compute the reliability of the terms and interaction terms (note that in this simulation I'm assuming that $cor(X,MZ) = cor(M,XZ) = cor(Z,XM) = 0$, which may not actually be the case in truth, this is why I'd like to know the formula for the reliability.

add.noise = function(x,rel){
  x2 = base::scale(x = (x+stats::rnorm(length(x), 0, sqrt((1-rel)/rel))),
                   center = TRUE, scale = TRUE)
  return(x2)
}

estimate_three_way_rel=function(r.X.M,r.X.Z,r.M.Z,rel.X,rel.M,rel.Z,N){
  
  # Simulate two identical data sets
  
  cor.mat = matrix(data = c(1,r.X.M,r.X.Z,
                            r.X.M,1,r.M.Z,
                            r.X.Z,r.M.Z,1),nrow = 3,byrow = TRUE)
  
  d1 = MASS::mvrnorm(n=N, mu=rep(0,3), Sigma=cor.mat, empirical=TRUE)
  d2 = d1
  
  # add in random noise proportional to the reliability
  
  d1[,1] = add.noise(x = d1[,1],rel = rel.X)
  d1[,2] = add.noise(x = d1[,2],rel = rel.M)
  d1[,3] = add.noise(x = d1[,3],rel = rel.Z)
  
  d2[,1] = add.noise(x = d2[,1],rel = rel.X)
  d2[,2] = add.noise(x = d2[,2],rel = rel.M)
  d2[,3] = add.noise(x = d2[,3],rel = rel.Z)
  
  # compute reliabilities
  
  rel.x = cor(d1[,1],d2[,1])
  rel.m = cor(d1[,2],d2[,2])
  rel.z = cor(d1[,3],d2[,3])
  rel.xm = cor(c(d1[,1]*d1[,2]),c(d2[,1]*d2[,2]))
  rel.xz = cor(c(d1[,1]*d1[,3]),c(d2[,1]*d2[,3]))
  rel.mz = cor(c(d1[,2]*d1[,3]),c(d2[,2]*d2[,3]))
  rel.xmz = cor(c(d1[,1]*d1[,2]*d1[,3]),c(d2[,1]*d2[,2]*d2[,3]))
  
  # return reliabilities
  
  return(c(rel.x,rel.m,rel.z,rel.xm,rel.xz,rel.mz,rel.xmz))
  
}

Now compute the reliabilities for a simulated data set:

    demo.int = estimate_three_way_rel(
                                r.X.M=.3,
                                r.X.Z=.5,
                                r.M.Z=.1,
                                
                                rel.X=.8,
                                rel.M=.5,
                                rel.Z=.2,
                                
                                N=100000)
> demo.int[c(1:3)]
[1] 0.7989486 0.4997133 0.2002855

The reliability of the $X$, $M$, and $Z$ terms match with that was specified.

The formula for two-way reliability yields:

adjust.cor = function(cor1,rel1,rel2){
  newcor = cor1*sqrt(rel1*rel2)
  return(newcor)
}

two.way.rel = function(cor1,rel1,rel2){
  newrel=((rel1*rel2)+cor1^2)/(1+cor1^2)
  return(newrel)
}

interaction.rel = function(r.X.M,r.X.Z,r.M.Z,rel.X,rel.M,rel.Z){
  
  r.X.M.new=adjust.cor(cor1 = r.X.M, rel1 = rel.X, rel2 = rel.M)
  r.X.Z.new=adjust.cor(cor1 = r.X.Z, rel1 = rel.X, rel2 = rel.Z)
  r.M.Z.new=adjust.cor(cor1 = r.M.Z, rel1 = rel.M, rel2 = rel.Z)
  
  rel.XM = two.way.rel(cor1 = r.X.M.new, rel1 = rel.X, rel2 = rel.M)
  rel.XZ = two.way.rel(cor1 = r.X.Z.new, rel1 = rel.X, rel2 = rel.Z)
  rel.MZ = two.way.rel(cor1 = r.M.Z.new, rel1 = rel.M, rel2 = rel.Z)
  
 return(c(rel.XM,rel.XZ,rel.MZ)) 
  
}

interaction.rel(  r.X.M=.3,
              r.X.Z=.5,
              r.M.Z=.1,
              
              rel.X=.8,
              rel.M=.5,
              rel.Z=.2)

[1] 0.4208494 0.1923077 0.1008991

This matches the reliability of the two-way interaction terms from the simulated data set:

demo.int[c(4:6)]
[1] 0.41461804 0.18411966 0.09994802

Now, the reliability of the three-way interaction term in this simulation is:

demo.int[7]
[1] 0.1125685

I'd like to know this value without having to run a simulation.

Improve this question
$\endgroup$
11
  • 1
    $\begingroup$ What is the definition of "reliability"? Is it just an alias for "standard deviation"? $\endgroup$
    – Zhanxiong
    Commented Feb 18, 2023 at 18:23
  • 1
    $\begingroup$ No, it is measurement consistency (see Wikipedia for example: en.m.wikipedia.org/wiki/Reliability_(statistics)). You might say it is the proportion of the variance and is not attributable to measurement, error or noise. $\endgroup$
    – David B
    Commented Feb 18, 2023 at 18:34
  • 1
    $\begingroup$ Have you verified that your equation for $\rho_{XM}$ matches equation 19 of your link 2? I gave it a few minutes of thought and concluded that they did not match. I think the math contained in link 2 is exactly what needs to be done to give the reliability of a three variable product. $\endgroup$
    – R Carnell
    Commented Feb 28, 2023 at 17:54
  • 1
    $\begingroup$ Thanks @RCarnell. I'm assuming the mean=0 for both X and M, once you account for that then it simplifies to the equation I gave. I agree that probably I need to expand/follow the math, unfortunately I'm not able to follow all of it. One major piece is that I'm not sure what C^2(X,Y) means (eq.s 14-18 in link 2), and so I can't follow why it simplifies to cor(x,y)^2 (eq. 19). $\endgroup$
    – David B
    Commented Feb 28, 2023 at 18:24
  • 1
    $\begingroup$ Hi @RCarnell. I did find one paper that solves for $V(XMY)$ - ccsenet.org/journal/index.php/ijsp/article/view/62832 (bottom of pg. 3). Is that helpful enough to avoid numerical integration? $\endgroup$
    – David B
    Commented Mar 1, 2023 at 15:32

1 Answer 1

Reset to default
3
+50
$\begingroup$

Using the result for the variance of three way interaction of correlated variables here:

(shortened version)

$$V(XMZ) = \sigma^2_{XMZ} = E(X^2M^2Z^2) - [E(XMZ)]^2$$

(reference notation)

$$ = \sigma_{11}(\sigma_{22}\sigma_{33} + 2\sigma_{23}^2) + 2\sigma_{12}(\sigma_{12}\sigma_{33} + 2\sigma_{13}\sigma_{23}) + 2\sigma_{13}(2\sigma_{12}\sigma_{23} + \sigma_{13}\sigma_{22})$$

(notation here)

$$ = \sigma_X^2(\sigma_M^2\sigma_Z^2 + 2\sigma_{MZ}^2) + 2\sigma_{XM}(\sigma_{XM}\sigma_Z^2 + 2\sigma_{XZ}\sigma_{MZ}) + 2\sigma_{XZ}(2\sigma_{XM}\sigma_{MZ} + \sigma_{XZ}\sigma_M^2)$$

Using the notation of this reference:

$$\rho_X = \frac{\sigma_{T,X}^2}{\sigma_{X}^2}$$

$$\rho_{XMZ} = \frac{\sigma_{T,XMZ}^2}{\sigma_{XMZ}^2}$$

$$ = \frac{\sigma_{T,X}^2(\sigma_{T,M}^2\sigma_{T,Z}^2 + 2\sigma_{MZ}^2) + 2\sigma_{XM}(\sigma_{XM}\sigma_{T,Z}^2 + 2\sigma_{XZ}\sigma_{MZ}) + 2\sigma_{XZ}(2\sigma_{XM}\sigma_{MZ} + \sigma_{XZ}\sigma_{T,M}^2)}{\sigma_X^2(\sigma_M^2\sigma_Z^2 + 2\sigma_{MZ}^2) + 2\sigma_{XM}(\sigma_{XM}\sigma_Z^2 + 2\sigma_{XZ}\sigma_{MZ}) + 2\sigma_{XZ}(2\sigma_{XM}\sigma_{MZ} + \sigma_{XZ}\sigma_M^2)}$$

Simulation

Starting with my own simulation...

set.seed(1403949)
N <- 100000
S <- matrix(c(1, .3, .7,
              .3, 3, .4,
              .7, .4, 5), nrow = 3, byrow = TRUE)
error <- 2
Xt <- MASS::mvrnorm(N, mu = c(0,0,0), Sigma = S)
X <- Xt + rnorm(3*N, 0, sqrt(error))

## Check Sample Variance vs Theoretical variance

var(Xt)
#>           [,1]      [,2]      [,3]
#> [1,] 0.9979354 0.2911005 0.7033502
#> [2,] 0.2911005 2.9848208 0.4028383
#> [3,] 0.7033502 0.4028383 5.0007352
S
#>      [,1] [,2] [,3]
#> [1,]  1.0  0.3  0.7
#> [2,]  0.3  3.0  0.4
#> [3,]  0.7  0.4  5.0

var(X)
#>           [,1]      [,2]      [,3]
#> [1,] 2.9939103 0.2894055 0.6979344
#> [2,] 0.2894055 5.0110680 0.4188717
#> [3,] 0.6979344 0.4188717 7.0286428
S + diag(3)*error
#>      [,1] [,2] [,3]
#> [1,]  3.0  0.3  0.7
#> [2,]  0.3  5.0  0.4
#> [3,]  0.7  0.4  7.0

for (i in 1:2) 
{
  for (j in (i+1):3) 
  {
    print(paste(i, j, var(X[,i]*X[,j])))
    print((S[i,i]+error)*(S[j,j]+error) + S[i,j]^2)
  }
}
#> [1] "1 2 15.1145101901148"
#> [1] 15.09
#> [1] "1 3 21.6049535467542"
#> [1] 21.49
#> [1] "2 3 35.0236533929537"
#> [1] 35.16

var(apply(Xt, 1, prod))
#> [1] 19.64158
S[1,1]*(S[2,2]*S[3,3] + 2*S[2,3]^2) + 
  2*S[1,2]*(S[1,2]*S[3,3] + 2*S[1,3]*S[2,3]) + 
  2*S[1,3]*(2*S[1,2]*S[2,3] + S[1,3]*S[2,2])
#> [1] 19.832

var(apply(X, 1, prod))
#> [1] 112.8418
(S[1,1]+error)*((S[2,2]+error)*(S[3,3]+error) + 2*S[2,3]^2) + 
  2*S[1,2]*(S[1,2]*(S[3,3] + error) + 2*S[1,3]*S[2,3]) + 
  2*S[1,3]*(2*S[1,2]*S[2,3] + S[1,3]*(S[2,2] + error))
#> [1] 112.792

## Check sample reliability vs theoretical reliability

1 - error / (S[1,1] + error)
#> [1] 0.3333333
var(Xt[,1]) / var(X[,1])
#> [1] 0.3333217
cor(Xt[,1], X[,1])^2
#> [1] 0.3335976

(S[1,1]*S[2,2] + S[1,2]^2) / ((S[1,1] + error)*(S[2,2] + error) + S[1,2]^2)
#> [1] 0.2047714
var(Xt[,1]*Xt[,2]) / var(X[,1]*X[,2])
#> [1] 0.2038291
cor(Xt[,1]*Xt[,2], X[,1]*X[,2])^2
#> [1] 0.2067311

var(apply(Xt, 1, prod)) / var(apply(X, 1, prod))
#> [1] 0.1740631
cor(apply(Xt, 1, prod), apply(X, 1, prod))^2
#> [1] 0.1795337
(S[1,1]*(S[2,2]*S[3,3] + 2*S[2,3]^2) + 
    2*S[1,2]*(S[1,2]*S[3,3] + 2*S[1,3]*S[2,3]) + 
    2*S[1,3]*(2*S[1,2]*S[2,3] + S[1,3]*S[2,2])) / 
  ((S[1,1]+error)*((S[2,2]+error)*(S[3,3]+error) + 2*S[2,3]^2) + 
     2*S[1,2]*(S[1,2]*(S[3,3] + error) + 2*S[1,3]*S[2,3]) + 
     2*S[1,3]*(2*S[1,2]*S[2,3] + S[1,3]*(S[2,2] + error)))
#> [1] 0.1758281

Continuing with the OP example:

set.seed(29389445)

add.noise = function(x,rel){
  x2 = base::scale(x = (x+stats::rnorm(length(x), 0, sqrt((1-rel)/rel))),
                   center = T,scale = T)
  return(x2)
}

estimate_three_way_rel=function(r.X.M,r.X.Z,r.M.Z,rel.X,rel.M,rel.Z,N){
  
  # Simulate two identical data sets
  
  cor.mat = matrix(data = c(1,r.X.M,r.X.Z,
                            r.X.M,1,r.M.Z,
                            r.X.Z,r.M.Z,1),nrow = 3,byrow = TRUE)
  
  d1 = MASS::mvrnorm(n=N, mu=rep(0,3), Sigma=cor.mat, empirical=TRUE)
  d2 = d1
  
  # add in random noise proportional to the reliability
  
  d1[,1] = add.noise(x = d1[,1],rel = rel.X)
  d1[,2] = add.noise(x = d1[,2],rel = rel.M)
  d1[,3] = add.noise(x = d1[,3],rel = rel.Z)
  
  d2[,1] = add.noise(x = d2[,1],rel = rel.X)
  d2[,2] = add.noise(x = d2[,2],rel = rel.M)
  d2[,3] = add.noise(x = d2[,3],rel = rel.Z)
  
  # compute reliabilities
  
  rel.x = cor(d1[,1],d2[,1])
  rel.m = cor(d1[,2],d2[,2])
  rel.z = cor(d1[,3],d2[,3])
  rel.xm = cor(c(d1[,1]*d1[,2]),c(d2[,1]*d2[,2]))
  rel.xz = cor(c(d1[,1]*d1[,3]),c(d2[,1]*d2[,3]))
  rel.mz = cor(c(d1[,2]*d1[,3]),c(d2[,2]*d2[,3]))
  rel.xmz = cor(c(d1[,1]*d1[,2]*d1[,3]),c(d2[,1]*d2[,2]*d2[,3]))
  
  # return reliabilities
  
  return(c(rel.x,rel.m,rel.z,rel.xm,rel.xz,rel.mz,rel.xmz))
  
}

demo.int = estimate_three_way_rel(
  r.X.M=.3,
  r.X.Z=.5,
  r.M.Z=.1,
  
  rel.X=.8,
  rel.M=.5,
  rel.Z=.2,
  
  N=100000)


three.way.rel <- function(r.X.M, r.X.Z, r.M.Z, rel.X, rel.M, rel.Z)
{
  S = matrix(data = c(1,r.X.M,r.X.Z,
                      r.X.M,1,r.M.Z,
                      r.X.Z,r.M.Z,1),nrow = 3,byrow = TRUE)
  
  (S[1,1]*(S[2,2]*S[3,3] + 2*S[2,3]^2) + 
      2*S[1,2]*(S[1,2]*S[3,3] + 2*S[1,3]*S[2,3]) + 
      2*S[1,3]*(2*S[1,2]*S[2,3] + S[1,3]*S[2,2])) /
    ((S[1,1]/rel.X)*((S[2,2]/rel.M)*(S[3,3]/rel.Z) + 2*S[2,3]^2) + 
       2*S[1,2]*(S[1,2]*(S[3,3]/rel.Z) + 2*S[1,3]*S[2,3]) + 
       2*S[1,3]*(2*S[1,2]*S[2,3] + S[1,3]*(S[2,2]/rel.M)))
}
three.way.rel(r.X.M=.3,
              r.X.Z=.5,
              r.M.Z=.1,
              
              rel.X=.8,
              rel.M=.5,
              rel.Z=.2)
#> [1] 0.1251289

demo.int[7]
#> [1] 0.1249188
Improve this answer
$\endgroup$
4
  • $\begingroup$ Amazing! Thank you! It's surprising to me that for the reliability of a 2-way interaction term I first had to adjust the individual correlations for reliability (via the adjust.cor() function), but your formula works without prior adjustment. Is the adjustment implicit in the math in some way that I'm not seeing? $\endgroup$
    – David B
    Commented Mar 2, 2023 at 19:11
  • 1
    $\begingroup$ The reason it was simpler in my formulation is that I stayed in the variance space instead of switching to the correlation and reliability space. The $S$ matrix gives the $\sigma_T$ values and the reliability then allows you to switch to the overall $\sigma^2 = \sigma_T^2 / \rho$ $\endgroup$
    – R Carnell
    Commented Mar 2, 2023 at 20:25
  • $\begingroup$ I see, thanks! Am I correct then that I should treat S as a variance-covariance matrix? $\endgroup$
    – David B
    Commented Mar 2, 2023 at 20:38
  • 1
    $\begingroup$ Yes. The S is a variance-covariance matrix. You had called it a correlation matrix in your code, and that was fine because of how you simulated data, but I used the variance-covariance matrix in my derivation. $\endgroup$
    – R Carnell
    Commented Mar 2, 2023 at 20:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.