#
#
#
############### Testing the TWO GROUP Bi-Factor Model with the ###############
###############  Satorra-Bentler (chi-square) Difference Test  ###############
#
#
# This script simulates and demonstrates grouped data (e.g. males & females) in 
# the context of a Bi-Factor Model which has a General factor (g) and 
# three Sub-factors (f1, f2, & f3). All four latent factors are orthogonal. 
# The "Bi-Factor" nomenclature refers to the fact that each manifest variable 
# is influenced by two orthogonal latent factors. 
# To see a 'general' graphical representation of the Bi-Factor model (without 
# groups), paste the following URL into your browser and hit the enter key:
# http://www.unt.edu/rss/class/Jon/R_SC/Module8/BiFactor_Figure.pdf

# Generate the data for GROUP 1 (all identified as "something.1"; such as Factor 2 
# in Group 1 = f2.1).

n.1  <- 5000
f1.1 <- rnorm(n.1)
f2.1 <- rnorm(n.1)
f3.1 <- rnorm(n.1)
g.1 <- rnorm(n.1)

lambda.1.1<-.9
lambda.1.2<-.8
lambda.1.3<-.7
lambda.1.4<-.6

x1  <- lambda.1.1*f1.1 + lambda.1.3*g.1 + rnorm(n.1,0,sqrt(1 - (lambda.1.1*lambda.1.3)^2))
x2  <- lambda.1.2*f1.1 + lambda.1.3*g.1 + rnorm(n.1,0,sqrt(1 - (lambda.1.2*lambda.1.3)^2))
x3  <- lambda.1.3*f1.1 + lambda.1.3*g.1 + rnorm(n.1,0,sqrt(1 - (lambda.1.3*lambda.1.3)^2))
x4  <- lambda.1.4*f1.1 + lambda.1.3*g.1 + rnorm(n.1,0,sqrt(1 - (lambda.1.4*lambda.1.3)^2))

x5  <- lambda.1.1*f2.1 + lambda.1.3*g.1 + rnorm(n.1,0,sqrt(1 - (lambda.1.1*lambda.1.3)^2))
x6  <- lambda.1.2*f2.1 + lambda.1.3*g.1 + rnorm(n.1,0,sqrt(1 - (lambda.1.2*lambda.1.3)^2))
x7  <- lambda.1.3*f2.1 + lambda.1.3*g.1 + rnorm(n.1,0,sqrt(1 - (lambda.1.3*lambda.1.3)^2))
x8  <- lambda.1.4*f2.1 + lambda.1.3*g.1 + rnorm(n.1,0,sqrt(1 - (lambda.1.4*lambda.1.3)^2))

x9  <- lambda.1.1*f3.1 + lambda.1.3*g.1 + rnorm(n.1,0,sqrt(1 - (lambda.1.1*lambda.1.3)^2))
x10 <- lambda.1.2*f3.1 + lambda.1.3*g.1 + rnorm(n.1,0,sqrt(1 - (lambda.1.2*lambda.1.3)^2))
x11 <- lambda.1.3*f3.1 + lambda.1.3*g.1 + rnorm(n.1,0,sqrt(1 - (lambda.1.3*lambda.1.3)^2))
x12 <- lambda.1.4*f3.1 + lambda.1.3*g.1 + rnorm(n.1,0,sqrt(1 - (lambda.1.4*lambda.1.3)^2))

bifac.1 <- data.frame(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12)
rm(f1.1,f2.1,f3.1,g.1,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,
   lambda.1.1,lambda.1.2,lambda.1.3,lambda.1.4)

head(bifac.1)
summary(bifac.1)
nrow(bifac.1)

####### Using the 'lavaan' package to estimate the bi-factor model for GROUP 1. 

library(lavaan)

### Specify the Confirmatory Factor Analysis (CFA); which is to say, 
### Specify the Measurement Model.

m.model.1 <- '
  general.factor  =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 + x12
  factor1         =~ x1 + x2 + x3 + x4
  factor2         =~ x5 + x6 + x7 + x8
  factor3         =~ x9 + x10 + x11 + x12
  general.factor  ~~ 0*factor1
  general.factor  ~~ 0*factor2
  general.factor  ~~ 0*factor3
  factor1         ~~ 0*factor2
  factor1         ~~ 0*factor3
  factor2         ~~ 0*factor3'

### Fit the measurement model and get the summary. 

m.modelfit.1 <- cfa(m.model.1, data = bifac.1, std.lv = TRUE, information="observed")
summary(m.modelfit.1, fit.measures = TRUE, standardize = FALSE)

###########################################################################

# Generate the data for GROUP 2 (all identified as "something.2"; such as Factor 1 in Group 2 = f1.2).

n.2  <- 5000
f1.2 <- rnorm(n.2)
f2.2 <- rnorm(n.2)
f3.2 <- rnorm(n.2)
g.2 <- rnorm(n.2)

lambda.2.1<-.6 # In GROUP 1 this lambda was 0.9
lambda.2.2<-.8
lambda.2.3<-.7
lambda.2.4<-.9 # In GROUP 1 this lambda was 0.6

x1  <- lambda.2.1*f1.2 + lambda.2.3*g.2 + rnorm(n.2,0,sqrt(1 - (lambda.2.1*lambda.2.3)^2))
x2  <- lambda.2.2*f1.2 + lambda.2.3*g.2 + rnorm(n.2,0,sqrt(1 - (lambda.2.2*lambda.2.3)^2))
x3  <- lambda.2.3*f1.2 + lambda.2.3*g.2 + rnorm(n.2,0,sqrt(1 - (lambda.2.3*lambda.2.3)^2))
x4  <- lambda.2.4*f1.2 + lambda.2.3*g.2 + rnorm(n.2,0,sqrt(1 - (lambda.2.4*lambda.2.3)^2))

x5  <- lambda.2.1*f2.2 + lambda.2.3*g.2 + rnorm(n.2,0,sqrt(1 - (lambda.2.1*lambda.2.3)^2))
x6  <- lambda.2.2*f2.2 + lambda.2.3*g.2 + rnorm(n.2,0,sqrt(1 - (lambda.2.2*lambda.2.3)^2))
x7  <- lambda.2.3*f2.2 + lambda.2.3*g.2 + rnorm(n.2,0,sqrt(1 - (lambda.2.3*lambda.2.3)^2))
x8  <- lambda.2.4*f2.2 + lambda.2.3*g.2 + rnorm(n.2,0,sqrt(1 - (lambda.2.4*lambda.2.3)^2))

x9  <- lambda.2.1*f3.2 + lambda.2.3*g.2 + rnorm(n.2,0,sqrt(1 - (lambda.2.1*lambda.2.3)^2))
x10 <- lambda.2.2*f3.2 + lambda.2.3*g.2 + rnorm(n.2,0,sqrt(1 - (lambda.2.2*lambda.2.3)^2))
x11 <- lambda.2.3*f3.2 + lambda.2.3*g.2 + rnorm(n.2,0,sqrt(1 - (lambda.2.3*lambda.2.3)^2))
x12 <- lambda.2.4*f3.2 + lambda.2.3*g.2 + rnorm(n.2,0,sqrt(1 - (lambda.2.4*lambda.2.3)^2))

bifac.2 <- data.frame(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12)
rm(f1.2,f2.2,f3.2,g.2,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,
   lambda.2.1,lambda.2.2,lambda.2.3,lambda.2.4)

head(bifac.2)
summary(bifac.2)
nrow(bifac.2)

####### Using the 'lavaan' package to estimate the bi-factor model for GROUP 2. 

### Specify the Confirmatory Factor Analysis (CFA); which is to say, 
### Specify the Measurement Model.

m.model.2 <- '
  general.factor  =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 + x12
  factor1         =~ x1 + x2 + x3 + x4
  factor2         =~ x5 + x6 + x7 + x8
  factor3         =~ x9 + x10 + x11 + x12
  general.factor  ~~ 0*factor1
  general.factor  ~~ 0*factor2
  general.factor  ~~ 0*factor3
  factor1         ~~ 0*factor2
  factor1         ~~ 0*factor3
  factor2         ~~ 0*factor3'

### Fit the measurement model. 

m.modelfit.2 <- cfa(m.model.2, data = bifac.2, std.lv = TRUE, information="observed")
summary(m.modelfit.2, fit.measures = TRUE, standardize = FALSE)

###############################################################################################

####### Combining the two group's data into one 'data.frame' (df).

group <- as.factor(c(rep(1,n.1),rep(2,n.2)))
bifac.grps <- rbind(bifac.1,bifac.2)
bifac.grps <- cbind(group, bifac.grps)
bifac.grps.df<-data.frame(bifac.grps)
head(bifac.grps.df)
nrow(bifac.grps.df)

rm(group,bifac.grps,n.1,n.2)

####### Using 'lavaan' package to estimate the MULTI-GROUP bi-factor model.

### SPECIFY the measurement model (CFA).

m.model.3 <- '
general.factor  =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 + x12
factor1         =~ x1 + x2 + x3 + x4
factor2         =~ x5 + x6 + x7 + x8
factor3         =~ x9 + x10 + x11 + x12
general.factor  ~~ 0*factor1
general.factor  ~~ 0*factor2
general.factor  ~~ 0*factor3
factor1         ~~ 0*factor2
factor1         ~~ 0*factor3
factor2         ~~ 0*factor3
factor1         ~~ 1*factor1
factor2         ~~ 1*factor2
factor3         ~~ 1*factor3
general.factor  ~~ 1*general.factor'

### Testing measurement invariance across the groups.

#   Note: the group="group" command requests a multiple group analysis
#   where "group" is simply the name of the grouping variable which 
#   identifies groups '1' and '2'.
#   This command can be used in 'cfa' and 'sem' functions of 'lavaan'.

m.model.3.invariance.test <- measurementInvariance(m.model.3, data=bifac.grps.df, group="group")
m.model.3.invariance.test

### FIT the null and alternative models to the measurement model specified above.

# The NULL model (i.e. hypothesize the groups are the same; group "loadings" are constrained 
# to be equal).

m.modelfit.3.1 <- cfa(m.model.3, data = bifac.grps.df, std.lv = FALSE,
                      meanstructure=TRUE, std.ov=FALSE, se="robust",
                      orthogonal=TRUE,
                      estimator="MLR",
                      group="group", group.equal=c("loadings"))

#  Note: The 'standardize = TRUE' argument (below) returns standardized parameter
#  estimates. Also note: the loadings (in the summary output below) for each group
#  are the same. 

summary(m.modelfit.3.1, fit.measures = TRUE, standardize = TRUE)

fitMeasures(m.modelfit.3.1)

# Extract the model's robust chi-square value (4), degrees of freedom (5), & scaling 
# factor (7) from the 'fitMeasures' and assign those 3 parameters to a vector for 
# later use with the Satorra-Bentler (chi-square) Difference Test. 

null.sb <- as.vector(fitMeasures(m.modelfit.3.1)[c(4,5,7)])
null.sb

# The Alternative model (i.e. hypothesize the groups are different; groups uncontrained, 
# not equal or allowed to vary [group.equal=NULL]).

m.modelfit.3.2 <- cfa(m.model.3, data = bifac.grps.df, std.lv = FALSE,
                      meanstructure=TRUE, std.ov=FALSE, se="robust",
                      orthogonal=TRUE,
                      estimator="MLR",
                      group="group", group.equal=NULL)

#  Note: in this summary output (below), the groups have different loadings. 

summary(m.modelfit.3.2, fit.measures = TRUE, standardize = TRUE)

# Extract the model's robust chi-square value, degrees of freedom (df), & scaling 
# factor from the 'fitMeasures' and assign those 3 parameters to a vector for 
# later use with the Satorra-Bentler (chi-square) Difference Test. 

alt.sb <- as.vector(fitMeasures(m.modelfit.3.2)[c(4,5,7)])
alt.sb

##############################################################################

# Create the Satorra-Bentler (chi-square) Difference Test function to determine 
# if there is a difference in model fit between the NULL model (m.model.fit.3.1) 
# and the Alternative model (m.model.fit.3.2).
# Make sure to highlight & submit the entire function from { to } which 
# includes "Notes" (comments).

SatorraBentlerDiffTest <- function(c0=NULL,c1=NULL,
                                   d0=NULL,d1=NULL,
                                   T0=NULL,T1=NULL,
                                tot.N=NULL,num.groups=NULL){
# Notes:
#   Calculation of the correct chi-squared difference test statistic for
#   models estimated with the MLM or MLR estimator
#
#          c0 - scaling correction factor for the null model (m.modelfit.3.1)
#          c1 - scaling correction factor for the alternative model. (m.modelfit.3.2)
#          d0 - degrees of freedom for the null model
#          d1 - degrees of freedom for the alternative model.
#          T0 - Satorra-Bentler scaled chi-squared value for the null model.
#          T1 - Satorra-Bentler scaled chi-squared value for the alternative model.
#       tot.N - Total sample size
#   cramers.v - Effects size based on chi-square change statistic (T.statistic)
#               and total sample size (N)
#
  cd.delta <- (d0 * c0 - d1 * c1)/(d0 - d1)
  T.statistic <- (T0 * c0 - T1 * c1)/cd.delta
  T.df <- d0 - d1
  T.statistic.df.prob <- 1-pchisq(T.statistic, df=T.df)
  T.statistic.df.prob.results<-c(cd.delta,T.statistic,T.df,T.statistic.df.prob)
  cramers.v<-sqrt(T.statistic/(tot.N*(num.groups-1)))
  adj.cramers.v<-sqrt(T.statistic/(tot.N*(num.groups-1)*T.df))
  T.RFI = ((T0 * c0 - T1 * c1)/cd.delta)/((T0 * c0)/cd.delta)
  return(c(T.statistic.df.prob.results,cramers.v,adj.cramers.v,T.RFI))
}

# Using the extracted chi-square, df, and scaling correction from eariler, as well 
# as the total sample size (nrow of the combined data) and the number of groups, 
# we can perform the Satorra-Bentler (chi-square) Difference Test. 

SatorraBentlerDiffTest(T0 = null.sb[1], T1 = alt.sb[1], d0 = null.sb[2], d1 = alt.sb[2], 
                       c0 = null.sb[3], c1 = alt.sb[3], tot.N = nrow(bifac.grps.df), num.groups = 2)

# The function above returns: scaling factor, chi-square change or difference, degrees of freedom, 
# p-value, Cramer's V effect size, an adjusted Cramer's V effect size (adjusted for the df of the 
# Satorra-Bentler Difference Test), and a Relative Fit Index (T.RFI) for the SB Difference 
# Test. The T.RFI is a Relative Fit Index, common examples of RFI are the TLI (Tucker Lewis 
# Index) and the CFI (Comparative Fit Index). Here the T.RFI can be thought of as a Proportional 
# Model Fit improvement metric (larger is better).
# If the p-value is less than 0.001 then reject Null model and conclude that the groups are 
# significantly different. The T.RFI ranges from (approximately) zero to 1.00. It can occasionally 
# be negative when there is no difference between groups (i.e. random variation around zero with very 
# small values; e.g. +/-0.002). 








# END; Last updated Feb. 24, 2011.