# # # ############### 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. # This function has been moved to another package: "semTools" # 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'. library(semTools) 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 (6), degrees of freedom (7), & scaling # factor (9) 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(6,7,9)]) 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.partial=c("loadings")) # 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(6,7,9)]) 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.