Skip to content

Instantly share code, notes, and snippets.

@tkasci
Created December 30, 2018 16:01
Show Gist options
  • Save tkasci/e2a485942f2e95ff7f533856f4184dc3 to your computer and use it in GitHub Desktop.
Save tkasci/e2a485942f2e95ff7f533856f4184dc3 to your computer and use it in GitHub Desktop.
Compare two correlation matrices in R
# From this thread on ResearchGate: https://www.researchgate.net/post/Significance_of_the_difference_between_two_correlation_matrices
# Slightly modified by me
# A test if two correlation matrices are equal
# Correlation matries obtained from:
#
# Ryan, J. J., Paolo, A. M., Miller, D. A., & Morris, J. (1997).
# Exploratory factor analysis of the Wechsler Adult Intelligence Scale-Revised in a sample of brain-damaged women.
# Archives of Clinical Neuropsychology, 12(7), 683-689.
#
# and
#
# Allen, D. N., Huegel, S. G., Seaton, B. E., Goldstein, G., Gurklis, J. A., & Van Kammen, D. P. (1998).
# Confirmatory factor analysis of the WAIS-R in patients with schizophrenia.
# Schizophrenia research, 34(1), 87-94.
# load R structural equation modeling package lavaan
library( lavaan )
# Varible names
# Replace with yours (separated by a space, no tabs )
varnames <- c("Inf", "DSp", "Voc", "Ari", "Com", "Sim", "PC", "PA", "BD", "OA", "Dsy")
# The number of variables (check)
( nvars <- length( varnames ) )
# The number of observations in each study
# Replace with yours
nobs1 <- 152
nobs2 <- 196
# The correlation matrices
# replace with yours
cormat1 <- matrix(c(1.00,.44,.80,.57,.71,.68,.42,.37,.32,.37,.35,
.44,1.00,.49,.68,.49,.55,.43,.36,.45,.39,.53,
.80,.49,1.00,.56,.76,.70,.40,.38,.31,.38,.36,
.57,.68,.56,1.00,.58,.61,.59,.54,.57,.58,.60,
.71,.49,.76,.58,1.00,.67,.47,.45,.42,.46,.41,
.68,.55,.70,.61,.67,1.00,.63,.51,.53,.49,.52,
.42,.43,.40,.59,.47,.63,1.00,.71,.73,.71,.65,
.37,.36,.38,.54,.45,.51,.71,1.00,.65,.66,.61,
.32,.45,.31,.57,.42,.53,.73,.65,1.00,.77,.70,
.37,.39,.38,.58,.46,.49,.71,.66,.77,1.00,.65,
.35,.53,.36,.60,.41,.52,.65,.61,.70,.65,1.00), nvars, nvars, dimnames = list( varnames, varnames ) )
cormat2 <- matrix(c(1.00,.34,.72,.45,.52,.44,.42,.24,.37,.20,.32,
.34,1.00,.38,.43,.29,.31,.17,.19,.29,.13,.26,
.72,.38,1.00,.39,.70,.59,.35,.20,.27,.17,.23,
.45,.43,.39,1.00,.40,.23,.27,.34,.43,.16,.26,
.52,.29,.70,.40,1.00,.54,.43,.32,.27,.12,.19,
.44,.31,.59,.23,.54,1.00,.45,.31,.33,.33,.16,
.42,.17,.35,.27,.43,.45,1.00,.43,.55,.49,.21,
.24,.19,.20,.34,.32,.31,.43,1.00,.57,.49,.25,
.37,.29,.27,.43,.27,.33,.55,.57,1.00,.62,.36,
.20,.13,.17,.16,.12,.33,.49,.49,.62,1.00,.30,
.32,.26,.23,.26,.19,.16,.21,.25,.36,.30,1.00), nvars, nvars, dimnames = list( varnames, varnames ) )
# To make things easy, the variable names are changed into x1 to xn
newvarnames <- paste( 'x',1:nvars, sep='' )
# Add the new variable names to the matrices
dimnames( cormat1 ) <- dimnames( cormat2 ) <- list( newvarnames, newvarnames )
# Here's the model
model <- paste( 'X', 1:nvars, ' =~ ', 'x', 1:nvars, ';', sep='' )
# Fit the model
fit <- sem( model, sample.cov = list( cormat1, cormat2 ), sample.nobs = c( nobs1, nobs2 ), group.equal = 'lv.covariances' )
# This is the outcome of the (Chi^2) test
lavTestLRT( fit )[ 2, 5:7 ]
fitMeasures(fit, c("tli","nnfi","rmsea","srmr"))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment