#
#
#
###############                 SIMULATING THE               ###############
############### Satorra-Bentler (chi-square) Difference Test ###############
###############           TWO GROUP Bi-Factor Model.         ###############
#
#
# The purpose of this simulation is to show how the Satorra-Bentler (SB) Difference Test performs 
# as the loadings of two groups diverge. All loadings (i.e. both groups' loadings) start at 0.50;
# meaning the NULL situation is present, the two groups do not differ. With each iteration 
# the loadings of Group 1 increase by 0.05 and the loadings of Group 2 decrease by 0.05. Therefore, 
# in the second iteration (or loop), the loadings for Group 1 are set to 0.55 and the loadings 
# for Group 2 are set to 0.45. In the tenth (final) iteration the loadings for Group 1 are set 
# at 0.95 and the loadings for Group 2 are set at 0.05; thus they are almost as different as 
# they can be. Keep in mind, this simulation does not return the results of each Confirmatory 
# Factor Analysis; only the results associated with the SB Difference Test. 

library(lavaan)

# The 'SB.Results' object (a list) will contain the results of the SB Difference Test for 
# each iteration [[1]] through [[10]].

SB.Results <- as.list(0)

# The following two objects (z1 & z2) are the initial values for the loadings; 0.05 will be added 
# to z1 and 0.05 will be subtracted from z2 prior to data generation (in each iteration).
 
z1 <- 0.45
z2 <- 0.55

# The simulation (or 'For Loop') begins below. You may want to read through it, prior to running 
# it, to get an idea of what is going on. Be sure to highlight/submit the entire function; it is 
# quite lengthy and contains: generation of each group's data, combining them into one data frame, 
# specifying the measurement model (CFA), fitting the NULL and ALTernative hypothesis models to 
# the specified measurement model, extracting pieces of output from each, then the Satorra-Bentler 
# (chi-square) Difference Test; with 3 Effect Size measures. 

# <<<<<<<<<< BEGIN FUNCTION >>>>>>>>>>

for (i in 1:10){

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

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

z1 <- z1 + .05

lambda.1.1<-z1
lambda.1.2<-z1
lambda.1.3<-z1
lambda.1.4<-z1

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)

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

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

z2 <- z2 - .05

lambda.2.1<-z2
lambda.2.2<-z2
lambda.2.3<-z2
lambda.2.4<-z2

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)

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

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

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'

### 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 in this model).

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

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

# 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).

SatorraBentlerDiffTest <- function(c0=NULL,c1=NULL,
                                   d0=NULL,d1=NULL,
                                   T0=NULL,T1=NULL,
                                tot.N=NULL,num.groups=NULL){
  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)
  T.RFI = ((T0 * c0 - T1 * c1)/cd.delta)/((T0 * c0)/cd.delta)
  cramers.v<-sqrt(T.statistic/(tot.N*(num.groups-1)))
  adj.cramers.v<-sqrt(T.statistic/(tot.N*(num.groups-1)*T.df))
  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 (in 
# this example, 2 groups), we can perform the Satorra-Bentler (chi-square) Difference Test
# and assign the results to 'SB.Results'. 

SB.Results[[i]] <- 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)
}

# <<<<<<<<<< END FUNCTION >>>>>>>>>>

# The 'SB.Results' object (below) contains 7 indices of the difference between the two groups at each 
# of the 10 iterations. The first number listed is the scaling factor, then the chi-square 
# for the SatorraBentler (SB) Difference Test, the degrees of freedom (df) for the SB Difference 
# Test, then the p-value for the SB Difference Test, then three Effect Size measures; Crammer's V, 
# then an adjusted Crammer's V (taking into account the df from the SB Difference Test), then 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. 

SB.Results

# Notice in the SB.Results, the SB difference chi-square displays large increases with each 
# iteration, growing from 2 digits, to 3 digits, to 4 digits, just in the first three iterations. 
# Also, the p-value becomes 0.0000 at the second iteration and the effect size measures (all 3) 
# reflect the growing differences between the two groups with each iteration. The p-value is 
# clearly benefiting from the large sample size (tot.N = 10000). 
# Below is a table showing the loadings for each group at each iteration in the simulation:
#
#                Group 1          Group 2          Difference
# Iteration 1 :   0.50             0.50             0.00
# Iteration 2 :   0.55             0.45             0.10
# Iteration 3 :   0.60             0.40             0.20
# Iteration 4 :   0.65             0.35             0.30
# Iteration 5 :   0.70             0.30             0.40
# Iteration 6 :   0.75             0.25             0.50
# Iteration 7 :   0.80             0.20             0.60
# Iteration 8 :   0.85             0.15             0.70
# Iteration 9 :   0.90             0.10             0.80
# Iteration 10:   0.95             0.05             0.90

# Showing the final loadings for Group 1 (z1) and Group 2 (z2).

z1
z2



# House cleaning:

ls()
rm(alt.sb,bifac.1,bifac.2,bifac.grps.df,i,m.model.3,m.modelfit.3.1,m.modelfit.3.2,null.sb,
   SatorraBentlerDiffTest,SB.Results,z1,z2)
ls()    

search()
detach("package:lavaan")
search()





# END; Last updated Feb. 25, 2011.