Methods of scaling applied to easier examples

A one-factor model.

Model for Reference-Group Method

Model for Marker-Variable Method

Model for Effects-Scaling Method

OneFactor.r

#### Methods of Scaling and Identification

## Some easier examples.
## Demonstrates three methods of scaling in a one-factor model:
## 1. Reference-Group Method - Constrain latent variable's variance and mean;
## 2. Marker-Variable Method - Constrain one loading and that indicator's
##    intercept;
## 3. Effects-Scaling Method - Constrain sums of loadings and intercepts.

## Compare results with lavaan's results

## One Factor:
#  "Positive Affect" factor for 7th grade

# Data in Appendix A of:
# Little, T., Slegers, D., & Card, N. (2006). A non-arbitrary method of
# identifying and scaling latent variables in SEM and MACS models.
# Structural Equation Modeling, 13(1), 59-72.

## Load package
library(OpenMx)

## Get data
# Vectors of correlations (row-by-row), standard deviations, and means, 
# and sample size.
vcor <- c(
   1.00000,
   0.75854,  1.00000,
   0.76214,  0.78705,  1.00000)

vmean <- c(3.13552, 2.99061, 3.06945)
vsd <- c(0.66770, 0.68506, 0.70672)
n <- 380

# Variable names
names <- c("pos1", "pos2", "pos3")

# Get full correlation matrix
mcor <- matrix( , 3, 3)                          # Empty matrix
mcor[upper.tri(mcor, diag = TRUE)] <- vcor       # Fill the upper triangle
mcor <- pmax(mcor, t(mcor), na.rm = TRUE)        # Fill the lower triangle

# Get co/variance matrix
mcov <- outer(vsd, vsd) * mcor

# Name the rows and columns
dimnames(mcov) <- list(names, names); mcov

names(vmean) <- names   # OpenMx requires the means be named

# Get data into OpenMx format
data <- mxData(observed = mcov, type = "cov", means = vmean, numObs = n)


### Method 1: Reference-Group Method
## Constrain latent variance to 1
## Constrain latent mean to 0

## Collect the bits and pieces needed by OpenMx
# Loadings
loadings <- mxPath(from = "POS", to = names, arrows = 1,
   free = TRUE, values = 0.5,
   labels = c("lambda1", "lambda2", "lambda3"))

# Latent variance - Constrain variance to 1
varLatent <- mxPath(from = "POS", arrows = 2,
   free = FALSE, values = 1,
   labels = "phi")

# Latent mean - Constrain mean to 0
means <- mxPath(from = "one", to = "POS", arrows = 1,
   free = FALSE, values = 0,
   labels = "kappa")

# Residual variances
varRes <- mxPath(from = names, arrows = 2,
   free = TRUE, values = 1,
   labels = c("theta1", "theta2", "theta3"))

# Intercepts
intercepts <- mxPath(from = "one", to = names, arrows = 1,
   free = TRUE, values = 1,
   labels = c("tau1", "tau2", "tau3"))

## Setup the model
model1 <- mxModel("One Factor Model", type = "RAM",
   manifestVars = names, latentVars = "POS",
   data, loadings, varLatent, means, varRes, intercepts)

## Run the model and get summary
fit1 <- mxRun(model1)
summary(fit1)

# These models are just-identified.
# Number of variables is 3;
# Therefore, number of pieces of information in 
# co/variance matrix: (3 X 4) / 2 = 6
# plus 3 means = 9 pieces of information.
# Number of parameters:
#   3 loadings
#   3 redisual variances
#   3 intercepts
#   Total of 9 parameters

# Therefore, degrees of freedom = 0,
# chi square = 0,
# and other fit indices are either 1 or 0.
# There is a small discrepancy - chi sq is not quite 0.
# Note: the log likelihoods for the model and the saturated model differ.
# OpenMx needs to be told to estimate reference 
# (saturated and independence) models.

summary(fit1, refModels = mxRefModels(fit1, run = TRUE))
coef(fit1)

########################
## Check with lavaan
library(lavaan)

m1 <- "
  # Loadings
  POS =~ NA*lambda1*pos1 + lambda2*pos2 + lambda3*pos3

  # Latent variance - Constrained to 1
  POS ~~ 1*phi*POS

  # Latent mean - Constrained to 0
  POS ~ 0*kappa*1

  # Residual variances
  pos1 ~~ theta1*pos1
  pos2 ~~ theta2*pos2
  pos3 ~~ theta3*pos3

  # Intercepts 
  pos1 ~ tau1*1
  pos2 ~ tau2*1
  pos3 ~ tau3*1
"

m1_fit <- sem(m1, sample.cov = mcov, sample.nobs = n, sample.mean = vmean)
OpenMx <- coef(fit1)
lavaan <- coef(m1_fit)

# Get coefs in same order
lavaan <- lavaan[match(names(OpenMx), names(lavaan))]
cbind(OpenMx, lavaan)
########################


### Method 2: Marker-Variable Method
## Constrain third loading to 1
## Constrain third intercept to 0

## Collect the bits and pieces needed by OpenMx
# Loadings - Constrain 3rd loading to 1
loadings <- mxPath(from = "POS", to = names, arrows = 1,
   free = c(TRUE, TRUE, FALSE), values = c(0.5, 0.5, 1),
   labels = c("lambda1", "lambda2", "lambda3"))

# Latent variance
varLatent <- mxPath(from = "POS", arrows = 2,
   free = TRUE, values = 1,
   labels = "phi")

# Latent mean
means <- mxPath(from = "one", to = "POS", arrows = 1,
   free = TRUE, values = 1,
   labels = "kappa")

# Residual variances
varRes <- mxPath(from = names, arrows = 2,
   free = TRUE, values = 1,
   labels = c("theta1", "theta2", "theta3"))

# Intercepts - Constrain 3rd intercept to 0
intercepts <- mxPath(from = "one", to = names, arrows = 1,
   free = c(TRUE, TRUE, FALSE), values = c(1, 1, 0),
   labels = c("tau1", "tau2", "tau3"))

## Setup the model
model2 <- mxModel("One Factor Model", type = "RAM",
   manifestVars = names, latentVars = "POS",
   data, loadings, varLatent, means, varRes, intercepts)

## Run the model and get summary
fit2 <- mxRun(model2)
summary(fit2, refModels = mxRefModels(fit2, run = TRUE))
coef(fit2)

########################
## Check with lavaan
library(lavaan)

m2 <- "
  # Loadings - Constrain 3rd loading to 1
  POS =~ NA*lambda1*pos1 + lambda2*pos2 + 1*lambda3*pos3

  # Latent variance
  POS ~~ phi*POS

  # Latent mean
  POS ~ kappa*1

  # Residual variances
  pos1 ~~ theta1*pos1
  pos2 ~~ theta2*pos2
  pos3 ~~ theta3*pos3

  # Intercepts - Constrain 3rd intercept to 0
  pos1 ~ tau1*1
  pos2 ~ tau2*1
  pos3 ~ 0*tau3*1
"

m2_fit <- sem(m2, sample.cov = mcov, sample.nobs = n, sample.mean = vmean)
OpenMx <- coef(fit2)
lavaan <- coef(m2_fit)

# Get coefs in same order
lavaan <- lavaan[match(names(OpenMx), names(lavaan))]      
cbind(OpenMx, lavaan)
########################


### Method 3: Effects-Scaling Method
## Constrain sum of loadings to equal number of loadings
## Constrain sum of intercepts to equal 0

## Collect the bits and pieces needed by OpenMx
# Loadings
loadings <- mxPath(from = "POS", to = names, arrows = 1,
   free = TRUE, values = 0.5,
   labels = c("lambda1", "lambda2", "lambda3"))

# Latent variance
varLatent <- mxPath(from = "POS", arrows = 2,
   free = TRUE, values = 1,
   labels = "phi")

# Latent mean
means <- mxPath(from = "one", to = "POS", arrows = 1,
   free = TRUE, values = 1,
   labels = "kappa")

# Residual variances
varRes <- mxPath(from = names, arrows = 2,
   free = TRUE, values = 1,
   labels = c("theta1", "theta2", "theta3"))

# Intercepts
intercepts <- mxPath(from = "one", to = names, arrows = 1,
   free = TRUE, values = 1,
   labels = c("tau1", "tau2", "tau3"))

# Constraints
conLoad <- mxConstraint(lambda1 + lambda2 + lambda3 == 3)
conInter <- mxConstraint(tau1 + tau2 + tau3 == 0)

## Setup the model
model3 <- mxModel("One Factor Model", type = "RAM",
   manifestVars = names, latentVars = "POS",
   data, loadings, means, varLatent, varRes, intercepts,
   conLoad, conInter)

## Run the model and get summary
fit3 <- mxRun(model3)
summary(fit3, refModels = mxRefModels(fit3, run = TRUE))
coef(fit3)

########################
## Check with lavaan
library(lavaan)

m3 <- "
  # Loadings
  POS =~ NA*lambda1*pos1 + lambda2*pos2 + lambda3*pos3

  # Latent variance
  POS ~~ phi*POS

  # Latent mean
  POS ~ kappa*1

  # Residual variances
  pos1 ~~ theta1*pos1
  pos2 ~~ theta2*pos2
  pos3 ~~ theta3*pos3

  # Intercepts 
  pos1 ~ tau1*1
  pos2 ~ tau2*1
  pos3 ~ tau3*1

  # Constraints
  lambda1 + lambda2 + lambda3 == 3
  tau1 + tau2 + tau3 == 0
"

m3_fit <- sem(m3, sample.cov = mcov, sample.nobs = n, sample.mean = vmean)
OpenMx <- coef(fit3)
lavaan <- coef(m3_fit)

# Get coefs in same order
lavaan <- lavaan[match(names(OpenMx), names(lavaan))]
cbind(OpenMx, lavaan)
########################

A two-factor model

Model for Reference-Group Method

Model for Marker-Variable Method

Model for Effects-Scaling Method

TwoFactor.r

#### Methods of Scaling and Identification

## Some easier examples.
## Demonstrates three methods of scaling in two-factor model:
## 1. Reference-Group Method - Constrain latent variables' variances and means;
## 2. Marker-Variable Method - Constrain one loading and that indicator's
##    intercept in both factors;
## 3. Effects-Scaling Method - Constrain sums of loadings and intercepts
##    for both factors.

## Compare results with lavaan's results

## Two Factors
#  "Positive Affect" and "Negative Affect" factors for 7th grade

# Data in Appendix A of:
# Little, T., Slegers, D., & Card, N. (2006). A non-arbitrary method of
# identifying and scaling latent variables in SEM and MACS models.
# Structural Equation Modeling, 13(1), 59-72.

## Load package
library(OpenMx)

## Get data
# Vectors of correlations (row-by-row), standard deviations, and means,
# and sample size.
vcor <- c(
   1.00000,
   0.75854,  1.00000,
   0.76214,  0.78705,  1.00000,
   0.02766,  0.00973, -0.05762,  1.00000,
  -0.06112, -0.06105, -0.14060,  0.78501,  1.00000,
  -0.02222, -0.05180, -0.10250,  0.81616,  0.81076,  1.00000)

vmean <- c(3.13552, 2.99061, 3.06945, 1.70069, 1.52705, 1.54483)
vsd <- c(0.66770, 0.68506, 0.70672, 0.71418, 0.66320, 0.65276)
n <- 380

# Variable names
names <- c("pos1", "pos2", "pos3", "neg1", "neg2", "neg3")

# Get full correlation matrix
mcor <- matrix( , 6, 6)                          # Empty matrix
mcor[upper.tri(mcor, diag = TRUE)] <- vcor       # Fill the upper triangle
mcor <- pmax(mcor, t(mcor), na.rm = TRUE)        # Fill the lower triangle

# Get co/variance matrix
mcov <- outer(vsd, vsd) * mcor

# Name the rows and columns
dimnames(mcov) <- list(names, names)
mcov

names(vmean) <- names   # OpenMx requires the means be named

# Get data into OpenMx format
data <- mxData(observed = mcov, type = "cov", means = vmean, numObs = n)


### Method 1: Reference-Group Method
## Constrain latent variances to 1
## Constrain latent means to 0

## Collect the bits and pieces needed by OpenMx
# Loadings
loadings1 <- mxPath(from = "POS", to = c("pos1", "pos2", "pos3"), arrows = 1,
   free = TRUE, values = 0.5,
   labels = c("lambda1", "lambda2", "lambda3"))

loadings2 <- mxPath(from = "NEG", to = c("neg1", "neg2", "neg3"), arrows = 1,
   free = TRUE, values = 0.5,
   labels = c("lambda4", "lambda5", "lambda6"))

# Latent variances and covariance - constrain variances to 1
varLatent <- mxPath(from = c("POS", "NEG"), arrows = 2, 
   connect = "unique.pairs",
   free = c(FALSE, TRUE, FALSE), values = 1,
   labels = c("phi11", "phi12", "phi22"))

# Latent means - constrain means to 0
means <- mxPath(from = "one", to = c("POS", "NEG"), arrows = 1,
   free = FALSE, values = 0,
   labels = c("kappa1", "kappa2"))

# Residual variances
varRes <- mxPath(from = names, arrows = 2,
   free = TRUE, values = 1, 
   labels = c("theta1", "theta2", "theta3", "theta4", "theta5", "theta6"))

# Intercepts
intercepts <- mxPath(from = "one", to = names, arrows = 1,
   free = TRUE, values = 1,
   labels = c("tau1", "tau2", "tau3", "tau4", "tau5", "tau6"))

## Setup the model
model1 <- mxModel("Two Factor Model", type = "RAM",
   manifestVars = names, latentVars = c("POS", "NEG"),
   data, loadings1, loadings2, varLatent, means, varRes, intercepts)

## Run the model and get summary
fit1 <- mxRun(model1)
summary(fit1)

# Number of variables is 6;
# Number of pieces of information in covariance matrix: (6 X 7) / 2 = 21
# plus 6 means = 27 pieces of information
# Number of parameters:
#   6 loadings (3 per factor)
#   1 covariance between factors
#   6 residual variances (3 per factor)
#   6 intercepts (3 per factor)
#   Total of 19 parameters
 
# Therefore, degrees of freedom = 8,
# and chi sq and fit indices are calculated.
# But make sure OpenMx estimates reference models (saturated and independence)
# upon which to base calculations for chi sq and fit indices.

summary(fit1, refModels = mxRefModels(fit1, run = TRUE))
coef(fit1)

########################
## Check with lavaan
library(lavaan)

m1 <- "
  # Loadings
  POS =~ NA*lambda1*pos1 + lambda2*pos2 + lambda3*pos3
  NEG =~ NA*lambda4*neg1 + lambda5*neg2 + lambda6*neg3

  # Latent variances and covariance - constrain variances to 1
  POS ~~ 1*phi11*POS
  NEG ~~ 1*phi22*NEG
  POS ~ phi12*NEG

  # Latent means - constrain means to 0
  POS ~ 0*kappa1*1
  NEG ~ 0*kappa2*1

  # Residual variances
  pos1 ~~ theta1*pos1
  pos2 ~~ theta2*pos2
  pos3 ~~ theta3*pos3
  neg1 ~~ theta4*neg1
  neg2 ~~ theta5*neg2
  neg3 ~~ theta6*neg3

  # Intercepts 
  pos1 ~ tau1*1
  pos2 ~ tau2*1
  pos3 ~ tau3*1
  neg1 ~ tau4*1
  neg2 ~ tau5*1
  neg3 ~ tau6*1
"

m1_fit <- sem(m1, sample.cov = mcov, sample.nobs = n, sample.mean = vmean)
summary(m1_fit)
OpenMx <- coef(fit1)
lavaan <- coef(m1_fit)

# Get coefs in same order
lavaan <- lavaan[match(names(OpenMx), names(lavaan))]
cbind(OpenMx, lavaan)
########################


### Method 2: Marker-Variable Method
## Constrain 3rd loading for POS and 1st loading for NEG to 1
## Constrain 3rd intercept for POS and 1st intercept for NEG to 0

## Collect the bits and pieces needed by OpenMx
# Loadings - Constrain 3rd loading for POS & 1st loading for NEG to 1
loadings1 <- mxPath(from = "POS", to = c("pos1", "pos2", "pos3"), arrows = 1,
   free = c(TRUE, TRUE, FALSE), values = c(0.5, 0.5, 1),
   labels = c("lambda1", "lambda2", "lambda3"))

loadings2 <- mxPath(from = "NEG", to = c("neg1", "neg2", "neg3"), arrows = 1,
   free = c(FALSE, TRUE, TRUE), values = c(1, 0.5, 0.5),
   labels = c("lambda4", "lambda5", "lambda6"))

# Latent variances and covariance
varLatent <- mxPath(from = c("POS", "NEG"), arrows = 2, 
   connect = "unique.pairs",
   free = TRUE, values = c(1, 0.5, 1),
   labels = c("phi11", "phi12", "phi22"))

# Latent means
means <- mxPath(from = "one", to = c("POS", "NEG"), arrows = 1,
   free = TRUE, values = 1,
   labels = c("kappa1", "kappa2"))

# Residual variances
varRes <- mxPath(from = names, arrows = 2,
   free = TRUE, values = 1,
   labels = c("theta1", "theta2", "theta3", "theta4", "theta5", "theta6"))

# Intercepts - constrain 3rd intercept for POS & 1st intercept for NEG to 0
intercepts <- mxPath(from = "one", to = names, arrows = 1,
   free = c(TRUE, TRUE, FALSE, FALSE, TRUE, TRUE),
   values = c(1, 1, 0, 0, 1, 1),
   labels = c("tau1", "tau2", "tau3", "tau4", "tau5", "tau6"))

## Setup the model
model2 <- mxModel("Two Factor Model", type = "RAM",
   manifestVars = names, latentVars = c("POS", "NEG"),
   data, loadings1, loadings2, varLatent, means, varRes, intercepts)

## Run the model and get summary
fit2 <- mxRun(model2)
summary(fit2, refModels = mxRefModels(fit2, run = TRUE))
coef(fit2)

########################
## Check with lavaan
library(lavaan)

m2 <- "
  # Loadings - Constrain 3rd loading for POS & 1st loading for NEG to 1
  POS =~ NA*lambda1*pos1 + lambda2*pos2 + 1*lambda3*pos3
  NEG =~  1*lambda4*neg1 + lambda5*neg2 + lambda6*neg3

  # Latent variances and covariance
  POS ~~ phi11*POS
  NEG ~~ phi22*NEG
  POS ~~ phi12*NEG

  # Latent means 
  POS ~ kappa1*1
  NEG ~ kappa2*1

  # Residual variances
  pos1 ~~ theta1*pos1
  pos2 ~~ theta2*pos2
  pos3 ~~ theta3*pos3
  neg1 ~~ theta4*neg1
  neg2 ~~ theta5*neg2
  neg3 ~~ theta6*neg3  

  # Intercepts - Constrain 3rd intercept for POS & 1st intercept for NEG to 0
  pos1 ~ tau1*1
  pos2 ~ tau2*1
  pos3 ~ 0*tau3*1
  neg1 ~ 0*tau4*1
  neg2 ~ tau5*1
  neg3 ~ tau6*1
"

m2_fit <- sem(m2, sample.cov = mcov, sample.nobs = n, sample.mean = vmean)
OpenMx <- coef(fit2)
lavaan <- coef(m2_fit)

# Get coefs in same order
lavaan <- lavaan[match(names(OpenMx), names(lavaan))]
cbind(OpenMx, lavaan)
########################


### Method 3: Effects-Scaling Method
## Constrain sum of loadings to equal number of loadings
## Constrain sum of intercepts to equal 0

## Collect the bits and pieces needed by OpenMx
# Loadings
loadings1 <- mxPath(from = "POS", to = c("pos1", "pos2", "pos3"), arrows = 1,
   free = TRUE, values = 0.5,
   labels = c("lambda1", "lambda2", "lambda3"))

loadings2 <- mxPath(from = "NEG", to = c("neg1", "neg2", "neg3"), arrows = 1,
   free = TRUE, values = 0.5,
   labels = c("lambda4", "lambda5", "lambda6"))

# Latent variances and covariance
varFac <- mxPath(from = c("POS", "NEG"), arrows = 2, 
   connect = "unique.pairs",
   free = TRUE, values = 1, labels = c("phi11", "phi12", "phi22"))

# Latent means
means <- mxPath(from = "one", to = c("POS", "NEG"), arrows = 1,
   free = TRUE, values = 1,
   labels = c("kappa1", "kappa2"))

# Residual variances
varRes <- mxPath(from = names, arrows = 2,
   free = TRUE, values = 1,
   labels = c("theta1", "theta2", "theta3", "theta4", "theta5", "theta6"))

# Intercepts
intercepts <- mxPath(from = "one", to = names, arrows = 1,
   free = TRUE, values = 1,
   labels = c("tau1", "tau2", "tau3", "tau4", "tau5", "tau6"))

# Constraints
conLoadPOS <- mxConstraint(lambda1 + lambda2 + lambda3 == 3)
conLoadNEG <- mxConstraint(lambda4 + lambda5 + lambda6 == 3)
conInterPOS <- mxConstraint(tau1 + tau2 + tau3 == 0)
conInterNEG <- mxConstraint(tau4 + tau5 + tau6 == 0)

## Setup the model
model3 <- mxModel("Two Factor Model", type = "RAM",
   manifestVars = names, latentVars = c("POS", "NEG"),
   data, loadings1, loadings2, varLatent, means, varRes, intercepts,
   conLoadPOS, conLoadNEG, conInterPOS, conInterNEG)

## Run the model and get summary
fit3 <- mxRun(model3)
summary(fit3, refModels = mxRefModels(fit3, run = TRUE))
coef(fit3)

########################
## Check with lavaan
library(lavaan)

m3 <- "
  # Loadings
  POS =~ NA*lambda1*pos1 + lambda2*pos2 + lambda3*pos3
  NEG =~ NA*lambda4*neg1 + lambda5*neg2 + lambda6*neg3

  # Latent variances and covariance
  POS ~~ phi11*POS
  NEG ~~ phi22*NEG
  POS ~~ phi12*NEG

  # Latent means
  POS ~ kappa1*1
  NEG ~ kappa2*1

  # Residual variances
  pos1 ~~ theta1*pos1
  pos2 ~~ theta2*pos2
  pos3 ~~ theta3*pos3
  neg1 ~~ theta4*neg1
  neg2 ~~ theta5*neg2
  neg3 ~~ theta6*neg3

  # Intercepts 
  pos1 ~ tau1*1
  pos2 ~ tau2*1
  pos3 ~ tau3*1
  neg1 ~ tau4*1
  neg2 ~ tau5*1
  neg3 ~ tau6*1

  # Constraints
  lambda1 + lambda2 + lambda3 == 3
  lambda4 + lambda5 + lambda6 == 3

  tau1 + tau2 + tau3 == 0
  tau4 + tau5 + tau6 == 0
"

m3_fit <- sem(m3, sample.cov = mcov, sample.nobs = n, sample.mean = vmean)
summary(m3_fit)
OpenMx <- coef(fit3)
lavaan <- coef(m3_fit)

# Get coefs in same order
lavaan <- lavaan[match(names(OpenMx), names(lavaan))]
cbind(OpenMx, lavaan)
########################

A one-factor two-group model

Model for Reference-Group Method

Model for Marker-Variable Method

Model for Effects-Scaling Method

OneFactorTwoGroup.r

#### Methods of Scaling and Identification

## Some easier examples.
## Demonstrates three methods of scaling in one-factor, two-group model:
## 1. Reference-Group Method - Constrain latent variable's variance and mean;
## 2. Marker-Variable Method - Constrain one loading and that indicator's
##    intercept;
## 3. Effects-Scaling Method - Constrain sums of loadings and intercepts.

## Following Little et al's lead, assume strong metric invariance:
## corresponding loadings and intercepts constrained to equality across groups

## Compare results from OpenMx with lavaan's results

## One Factor, Two Groups
#  "Positive Affect" factor for 7th and 8th grades

# Data in Appendix A of:
# Little, T., Slegers, D., & Card, N. (2006). A non-arbitrary method of
# identifying and scaling latent variables in SEM and MACS models.
# Structural Equation Modeling, 13(1), 59-72.

## Load package
library(OpenMx)

## Get data
# Vectors of correlations (row-by-row), standard deviations, and means,
# and sample size.
# 7th grade
vcor7 <- c(
   1.00000,
   0.75854,  1.00000,
   0.76214,  0.78705,  1.00000)

vmean7 <- c(3.13552, 2.99061, 3.06945)
vsd7 <- c(0.66770, 0.68506, 0.70672)
n7 <- 380

# 8th grade
vcor8 <- c(
   1.00000,
   0.81366,  1.00000,
   0.84980,  0.83523,  1.00000)

vmean8 <- c(3.07338, 2.84716, 2.97882)
vsd8 <- c(0.70299, 0.71780, 0.76208)
n8 <- 379

# Variable names
names <- c("pos1", "pos2", "pos3")

# Get full correlation matrix for each Grade
mcor7 <- matrix( , 3, 3)                           # Empty matrix
mcor7[upper.tri(mcor7, diag = TRUE)] <- vcor7      # Fill the upper triangle
mcor7 <- pmax(mcor7, t(mcor7), na.rm = TRUE)       # Fill the lower triangle

mcor8 <- matrix( , 3, 3)                           # Empty matrix
mcor8[upper.tri(mcor8, diag = TRUE)] <- vcor8      # Fill the upper triangle
mcor8 <- pmax(mcor8, t(mcor8), na.rm = TRUE)       # Fill the lower triangle

# Get co/variance matrices
mcov7 <- outer(vsd7, vsd7) * mcor7
mcov8 <- outer(vsd8, vsd8) * mcor8

# Name the rows and columns
dimnames(mcov7) <- list(names, names)
dimnames(mcov8) <- list(names, names)
mcov7; mcov8

names(vmean7) <- names   # OpenMx requires the means be named
names(vmean8) <- names

# Put data into lists - used in lavaan analysis
mcov <- list("Grade 7" = mcov7, "Grade 8" = mcov8)
vmean <- list(vmean7, vmean8)
n <- list(n7, n8)

# Get data into OpenMx format
data7 <- mxData(observed = mcov7, type = "cov", means = vmean7, numObs = n7)
data8 <- mxData(observed = mcov8, type = "cov", means = vmean8, numObs = n8)


### Method 1: Reference-Group Method
## Constrain latent variance to 1
## Constrain latent mean to 0
## These constraints apply to Grade 7 only.

## Collect the bits and pieces needed by OpenMx
# Loadings
loadings <- mxPath(from = "POS", to = names, arrows = 1,
   free = TRUE, values = 0.5,
   labels = c("lambda1", "lambda2", "lambda3"))

# Latent variances - Constrain Grade 7 variance to 1
varLatent7 <- mxPath(from = "POS", arrows = 2,
   free = FALSE, values = 1, labels = "phi7")

varLatent8 <- mxPath(from = "POS", arrows = 2,
   free = TRUE, values = 1, labels = "phi8")

# Latent means - Constrain Grade 7 mean to 0
means7 <- mxPath(from = "one", to = "POS", arrows = 1,
   free = FALSE, values = 0, labels = "kappa7")

means8 <- mxPath(from = "one", to = "POS", arrows = 1,
   free = TRUE, values = 1, labels = "kappa8") 

# Residual variances
varRes7 <- mxPath(from = names, arrows = 2,
   free = TRUE, values = 1,
   labels = c("theta71", "theta72", "theta73"))

varRes8 <- mxPath(from = names, arrows = 2,
   free = TRUE, values = 1, 
   labels = c("theta81", "theta82", "theta83"))

# Intercepts
intercepts <- mxPath(from = "one", to = names, arrows = 1,
   free = TRUE, values = 1,
   labels = c("tau1", "tau2", "tau3"))

## Setup models for each Grade
modGr7 <- mxModel("Grade7", type = "RAM",
   manifestVars = names, latentVars = "POS",
   data7, loadings, varLatent7, means7, varRes7, intercepts)

modGr8 <- mxModel("Grade8", type = "RAM",
   manifestVars = names, latentVars = "POS",
   data8, loadings, varLatent8, means8, varRes8, intercepts)

## Combine the two models using "mxFitFunctionMultigroup()"
fun <- mxFitFunctionMultigroup(c("Grade7.fitfunction", "Grade8.fitfunction"))
model1 <- mxModel("One Factor Two Group Model", modGr7, modGr8, fun)

## Run the model and get summary
fit1 <- mxRun(model1)
summary(fit1)

# Number of variables: 3
# Number of pieces of information in co/variance matrix: (3 X 4) / 2 = 6
# plus 3 means = 9 pieces of information for each group;
# that is, 18 for the model

# Number of parameters:
#   3 loadings (constrained to equality across groups)
#   1 latent mean
#   1 latent variance
#   6 residual variances (3 per group)
#   3 intercepts (constrained to equality across groups)
#   Total of 14 parameters

# Therefore, degrees of freedom = 4,
# and chi sq and fit indices are calculated.
# But make sure OpenMx estimates the reference models
# upon which to base calculations for chi sq and fit indices.

summary(fit1, refModels = mxRefModels(fit1, run = TRUE))
coef(fit1)

########################
## Check with lavaan
library(lavaan)

m1 <- "
  # Loadings
  POS =~ c(NA,NA)*c(lambda1, lambda1)*pos1 + c(lambda2, lambda2)*pos2 + c(lambda3, lambda3)*pos3

  # Latent variances - Constrain Grade 7 variance to 1
  POS ~~ c(1, NA)*c(phi7, phi8)*POS

  # Latent means - Constrain Grade 7 mean to 0
  POS ~ c(0, NA)*c(kappa7, kappa8)*1

  # Residual variances
  pos1 ~~ c(theta71, theta81)*pos1
  pos2 ~~ c(theta72, theta82)*pos2
  pos3 ~~ c(theta73, theta83)*pos3

  # Intercepts
  pos1 ~ c(tau1, tau1)*1
  pos2 ~ c(tau2, tau2)*1
  pos3 ~ c(tau3, tau3)*1
"

m1_fit <- sem(m1, sample.cov = mcov, sample.nobs = n, sample.mean = vmean)
summary(m1_fit)
OpenMx <- coef(fit1)
lavaan <- coef(m1_fit)

# Get coefs in same order
lavaan <- lavaan[match(names(OpenMx), names(lavaan))]
cbind(OpenMx, lavaan)
########################


### Method: Marker-Variable Method
## Constrain third loading in POS to 1
## Constrain third intercept to 0
## With strong measurement invariance,
## these constraints apply to both Grades.

## Collect the bits and pieces needed by OpenMx
# Loadings - Constrain 3rd loading to 1
loadings <- mxPath(from = "POS", to = names, arrows = 1,
   free = c(TRUE, TRUE, FALSE), values = c(0.5, 0.5, 1),
   labels = c("lambda1", "lambda2", "lambda3"))

# Latent variances
varLatent7 <- mxPath(from = "POS", arrows = 2,
   free = TRUE, values = 1, labels = "phi7")

varLatent8 <- mxPath(from = "POS", arrows = 2,
   free = TRUE, values = 1, labels = "phi8")

# Latent means
means7 <- mxPath(from = "one", to = "POS", arrows = 1,
   free = TRUE, values = 1, labels = "kappa7")

means8 <- mxPath(from = "one", to = "POS", arrows = 1,
   free = TRUE, values = 1, labels = "kappa8")

# Residual variances
varRes7 <- mxPath(from = names, arrows = 2,
   free = TRUE, values = 1,
   labels = c("theta71", "theta72", "theta73"))

varRes8 <- mxPath(from = names, arrows = 2,
   free = TRUE, values = 1,
   labels = c("theta81", "theta82", "theta83"))

# Intercepts - Constrain 3rd intercept to 0
intercepts <- mxPath(from = "one", to = names, arrows = 1,
   free = c(TRUE, TRUE, FALSE), values = c(1, 1, 0),
   labels = c("tau1", "tau2", "tau3"))

## Setup models for each Grade
modGr7 <- mxModel("Grade7", type = "RAM",
   manifestVars = names, latentVars = "POS",
   data7, loadings, varLatent7, means7, varRes7, intercepts)

modGr8 <- mxModel("Grade8", type = "RAM",
   manifestVars = names, latentVars = "POS",
   data8, loadings, varLatent8, means8, varRes8, intercepts)

## Combine the two models using "mxFitFunctionMultigroup()"
fun <- mxFitFunctionMultigroup(c("Grade7.fitfunction", "Grade8.fitfunction"))
model2 <- mxModel("One Factor Two Group Model", modGr7, modGr8, fun)

## Run the model and get summary
fit2 <- mxRun(model2)
summary(fit2, refModels = mxRefModels(fit2, run = TRUE))
coef(fit2)

########################
## Check with lavaan
library(lavaan)

m2 <- "
  # Loadings - Constrain 3rd loading to 1 in both Grades
  POS =~ c(NA,NA)*c(lambda1, lambda1)*pos1 + c(lambda2, lambda2)*pos2 + c(1,1)*c(lambda3, lambda3)*pos3

  # Latent variances
  POS ~~ c(phi7, phi8)*POS

  # Latent means
  POS ~ c(kappa7, kappa8)*1

  # Residual variances
  pos1 ~~ c(theta71, theta81)*pos1
  pos2 ~~ c(theta72, theta82)*pos2
  pos3 ~~ c(theta73, theta83)*pos3

  # Intercepts - Constrain 3rd intercept to 0 in both Grades
  pos1 ~ c(tau1, tau1)*1
  pos2 ~ c(tau2, tau2)*1
  pos3 ~ c(0,0)*c(tau3, tau3)*1
"

m2_fit <- sem(m2, sample.cov = mcov, sample.nobs = n, sample.mean = vmean)
summary(m2_fit)
OpenMx <- coef(fit2)
lavaan <- coef(m2_fit)

# Get coefs in same order
lavaan <- lavaan[match(names(OpenMx), names(lavaan))]
cbind(OpenMx, lavaan)
########################


### Method 3
## Constrain sum of loadings to equal number of loadings
## Constrain sum of intercepts to equal 0
## With strong measurement invariance,
## these constraints apply to both Grades.

## Collect the bits and pieces needed by OpenMx
# Loadings
loadings <- mxPath(from = "POS", to = names, arrows = 1,
   free = TRUE, values = 0.5,
   labels = c("lambda1", "lambda2", "lambda3"))

# Latent variances
varLatent7 <- mxPath(from = "POS", arrows = 2,
   free = TRUE, values = 1, labels = "phi7")

varLatent8 <- mxPath(from = "POS", arrows = 2,
   free = TRUE, values = 1, labels = "phi8")

# Latent means
means7 <- mxPath(from = "one", to = "POS", arrows = 1,
   free = TRUE, values = 0, labels = "kappa7")

means8 <- mxPath(from = "one", to = "POS", arrows = 1,
   free = TRUE, values = 1, labels = "kappa8")

# Residual variances
varRes7 <- mxPath(from = names, arrows = 2,
   free = TRUE, values = 1,
   labels = c("theta71", "theta72", "theta73"))

varRes8 <- mxPath(from = names, arrows = 2,
   free = TRUE, values = 1,
   labels = c("theta81", "theta82", "theta83"))

# Intercepts
intercepts <- mxPath(from = "one", to = names, arrows = 1,
   free = TRUE, values = 1,
   labels = c("tau1", "tau2", "tau3"))

## Setup models for each Grade
modGr7 <- mxModel("Grade7", type = "RAM",
   manifestVars = names, latentVars = "POS",
   data7, loadings, varLatent7, means7, varRes7, intercepts)

modGr8 <- mxModel("Grade8", type = "RAM",
   manifestVars = names, latentVars = "POS",
   data8, loadings, varLatent8, means8, varRes8, intercepts)

## Constraints
conLoad <- mxConstraint(lambda1 + lambda2 + lambda3 == 3)
conInter <- mxConstraint(tau1 + tau2 + tau3 == 0)

## Combine the two models using "mxFitFunctionMultigroup()"
fun <- mxFitFunctionMultigroup(c("Grade7.fitfunction", "Grade8.fitfunction"))
model3 <- mxModel("One Factor Two Group Model", modGr7, modGr8,
   conLoad, conInter, fun)

# Note: Constraints are added to final model, 
# not to each of the Grade 7 and Grade 8 models;
# otherwise, OpenMx will count them as 4 constraints, 
# accounting for 4 degrees of freedom, instead of 2.

## Run the model and get summary
fit3 <- mxRun(model3)
summary(fit3, refModels = mxRefModels(fit3, run = TRUE))

# Counting degrees of freedom.
# Number of variables: 3
# Number of pieces of information in the covariance matrix: 6
# plus 3 means = 9 per group;
# that is, 18 for the model

# Number of parameters:
# 3 loadings (constrained to equality across the groups)
# 3 intercepts (constrained to equality across the groups)
# 6 residual variances (3 per group)
# 2 latent means (1 per group)
# 2 latent variances (1 per group)
# Total of 16
# but take away 2 for the 2 constraints = 14,
# resulting in 18 - 14 = 4 degrees of freedom.
#
# If the constraints were added to each Grade's model,
# OpenMx would have counted them as 4 constraints (even thought they are identical),
# resulting in 6 degrees of freedom for the model.

########################
## Check with lavaan
library(lavaan)

m3 <- "
  # Loadings
  POS =~ c(NA,NA)*c(lambda1, lambda1)*pos1 + c(lambda2, lambda2)*pos2 + c(lambda3, lambda3)*pos3

  # Latent variances
  POS ~~ c(phi7, phi8)*POS

  # Latent means
  POS ~ c(kappa7, kappa8)*1

  # Residual variances
  pos1 ~~ c(theta71, theta81)*pos1
  pos2 ~~ c(theta72, theta82)*pos2
  pos3 ~~ c(theta73, theta83)*pos3

  # Intercepts
  pos1 ~ c(tau1, tau1)*1
  pos2 ~ c(tau2, tau2)*1
  pos3 ~ c(tau3, tau3)*1

  # Constraints
  lambda1 + lambda2 + lambda3 == 3
  tau1 + tau2 + tau3 == 0
"

m3_fit <- sem(m3, sample.cov = mcov, sample.nobs = n, sample.mean = vmean)
summary(m3_fit)
OpenMx <- coef(fit3)
lavaan <- coef(m3_fit)

# Get coefs in same order
lavaan <- lavaan[match(names(OpenMx), names(lavaan))]
cbind(OpenMx, lavaan)
########################