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