Package 'semTools'

Implement Saturated Correlates with FIML

Description

Automatically add auxiliary variables to a lavaan model when using full information maximum likelihood (FIML) to handle missing data

Usage

auxiliary(model, data, aux, fun, ..., envir = getNamespace("lavaan"),
  return.syntax = FALSE)

lavaan.auxiliary(model, data, aux, ..., envir = getNamespace("lavaan"))

cfa.auxiliary(model, data, aux, ..., envir = getNamespace("lavaan"))

sem.auxiliary(model, data, aux, ..., envir = getNamespace("lavaan"))

growth.auxiliary(model, data, aux, ..., envir = getNamespace("lavaan"))

Arguments

model	The analysis model can be specified with 1 of 2 objects: lavaan lavaan::model.syntax() specifying a hypothesized model without mention of auxiliary variables in aux a parameter table, as returned by lavaan::parTable(), specifying the target model without auxiliary variables. This option requires these columns (and silently ignores all others): c("lhs","op","rhs","user","group","free","label","plabel","start")
data	data.frame that includes auxiliary variables as well as any observed variables in the model
aux	character. Names of auxiliary variables to add to model
fun	character. Name of a specific lavaan function used to fit model to data (i.e., "lavaan", "cfa", "sem", or "growth"). Only required for auxiliary.
...	Additional arguments to pass to ⁠fun=⁠.
envir	Passed to do.call().
return.syntax	logical indicating whether to return a character string of lavaan::model.syntax() that can be added to a target ⁠model=⁠ that is also a character string. This can be advantageous, for example, to use add saturated correlates to a blavaan model.

Details

These functions are wrappers around the corresponding lavaan functions. You can use them the same way you use lavaan::lavaan(), but you must pass your full data.frame to the data argument. Because the saturated-correlates approaches (Enders, 2008) treats exogenous variables as random, fixed.x must be set to FALSE. Because FIML requires continuous data (although nonnormality corrections can still be requested), no variables in the model nor auxiliary variables specified in aux can be declared as ordered.

Value

a fitted lavaan::lavaan object. Additional information is stored as a list in the ⁠@external⁠ slot:

• baseline.model. a fitted lavaan::lavaan object. Results of fitting an appropriate independence model for the calculation of incremental fit indices (e.g., CFI, TLI) in which the auxiliary variables remain saturated, so only the target variables are constrained to be orthogonal. See Examples for how to send this baseline model to lavaan::fitMeasures().

• aux. The character vector of auxiliary variable names.

• baseline.syntax. A character vector generated within the auxiliary function, specifying the baseline.model syntax.

Author(s)

Terrence D. Jorgensen (University of Amsterdam; [email protected])

References

Enders, C. K. (2008). A note on the use of missing auxiliary variables in full information maximum likelihood-based structural equation models. Structural Equation Modeling, 15(3), 434–448. doi:10.1080/10705510802154307

Examples

dat1 <- lavaan::HolzingerSwineford1939
set.seed(12345)
dat1$z <- rnorm(nrow(dat1))
dat1$x5 <- ifelse(dat1$z < quantile(dat1$z, .3), NA, dat1$x5)
dat1$x9 <- ifelse(dat1$z > quantile(dat1$z, .8), NA, dat1$x9)

targetModel <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"

## works just like cfa(), but with an extra "aux" argument
fitaux1 <- cfa.auxiliary(targetModel, data = dat1, aux = "z",
                         missing = "fiml", estimator = "mlr")

## with multiple auxiliary variables and multiple groups
fitaux2 <- cfa.auxiliary(targetModel, data = dat1, aux = c("z","ageyr","grade"),
                         group = "school", group.equal = "loadings")

## calculate correct incremental fit indices (e.g., CFI, TLI)
fitMeasures(fitaux2, fit.measures = c("cfi","tli"))
## NOTE: lavaan will use the internally stored baseline model, which
##       is the independence model plus saturated auxiliary parameters
lavInspect(fitaux2@external$baseline.model, "free")

---

Calculate average variance extracted

Description

Calculate average variance extracted (AVE) per factor from lavaan object

Usage

AVE(object, obs.var = TRUE, omit.imps = c("no.conv", "no.se"),
  omit.factors = character(0), dropSingle = TRUE, return.df = TRUE)

Arguments

object	A lavaan::lavaan or lavaan.mi::lavaan.mi object, expected to contain only exogenous common factors (i.e., a CFA model). Cross-loadings are not allowed and will result in NA for any factor with indicator(s) that cross-load.
obs.var	logical indicating whether to compute AVE using observed variances in the denominator. Setting FALSE triggers using model-implied variances in the denominator.
omit.imps	character vector specifying criteria for omitting imputations from pooled results. Can include any of c("no.conv", "no.se", "no.npd"), the first 2 of which are the default setting, which excludes any imputations that did not converge or for which standard errors could not be computed. The last option ("no.npd") would exclude any imputations which yielded a nonpositive definite covariance matrix for observed or latent variables, which would include any "improper solutions" such as Heywood cases. NPD solutions are not excluded by default because they are likely to occur due to sampling error, especially in small samples. However, gross model misspecification could also cause NPD solutions, users can compare pooled results with and without this setting as a sensitivity analysis to see whether some imputations warrant further investigation.
omit.factors	character vector naming any common factors modeled in object whose indicators' AVE is not of interest.
dropSingle	logical indicating whether to exclude factors defined by a single indicator from the returned results. If TRUE (default), single indicators will still be included in the total column when return.total = TRUE.
return.df	logical indicating whether to return reliability coefficients in a data.frame (one row per group/level), which is possible when every model block includes the same factors (after excluding those in omit.factors and applying dropSingle).

Details

The average variance extracted (AVE) can be calculated by

$AVE = \frac{\bold{1}^\prime \textrm{diag}\left(\Lambda\Psi\Lambda^\prime\right)\bold{1}}{\bold{1}^\prime \textrm{diag}\left(\hat{\Sigma}\right) \bold{1}},$

Note that this formula is modified from Fornell & Larcker (1981) in the case that factor variances are not 1. The proposed formula from Fornell & Larcker (1981) assumes that the factor variances are 1. Note that AVE will not be provided for factors consisting of items with dual loadings. AVE is the property of items but not the property of factors. AVE is calculated with polychoric correlations when ordinal indicators are used.

Value

numeric vector of average variance extracted from indicators per factor. For models with multiple "blocks" (any combination of groups and levels), vectors may be returned as columns in a data.frame with additional columns indicating the group/level (see ⁠return.df=⁠ argument description for caveat).

Author(s)

Terrence D. Jorgensen (University of Amsterdam; [email protected])

References

Fornell, C., & Larcker, D. F. (1981). Evaluating structural equation models with unobservable variables and measurement errors. Journal of Marketing Research, 18(1), 39–50. doi:10.2307/3151312

See Also

compRelSEM() for composite reliability estimates

Examples

data(HolzingerSwineford1939)
HS9 <- HolzingerSwineford1939[ , c("x7","x8","x9")]
HSbinary <- as.data.frame( lapply(HS9, cut, 2, labels=FALSE) )
names(HSbinary) <- c("y7","y8","y9")
HS <- cbind(HolzingerSwineford1939, HSbinary)

HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ y7 + y8 + y9 '

fit <- cfa(HS.model, data = HS, ordered = c("y7","y8","y9"), std.lv = TRUE)

## works for factors with exclusively continuous OR categorical indicators
AVE(fit) # uses observed (or unconstrained polychoric/polyserial) by default
AVE(fit, obs.var = FALSE)


## works for multigroup models and for multilevel models (and both)
data(Demo.twolevel)
## assign clusters to arbitrary groups
Demo.twolevel$g <- ifelse(Demo.twolevel$cluster %% 2L, "type1", "type2")
model2 <- ' group: type1
  level: within
    fac =~ y1 + L2*y2 + L3*y3
  level: between
    fac =~ y1 + L2*y2 + L3*y3

group: type2
  level: within
    fac =~ y1 + L2*y2 + L3*y3
  level: between
    fac =~ y1 + L2*y2 + L3*y3
'
fit2 <- sem(model2, data = Demo.twolevel, cluster = "cluster", group = "g")
AVE(fit2)

---

Class For the Results of Bollen-Stine Bootstrap with Incomplete Data

Description

This class contains the results of Bollen-Stine bootstrap with missing data.

Usage

## S4 method for signature 'BootMiss'
show(object)

## S4 method for signature 'BootMiss'
summary(object)

## S4 method for signature 'BootMiss'
hist(x, ..., alpha = 0.05, nd = 2,
  printLegend = TRUE, legendArgs = list(x = "topleft"))

Arguments

object, x	object of class BootMiss
...	Additional arguments to pass to graphics::hist()
alpha	alpha level used to draw confidence limits
nd	number of digits to display
printLegend	logical. If TRUE (default), a legend will be printed with the histogram
legendArgs	list of arguments passed to the graphics::legend() function. The default argument is a list placing the legend at the top-left of the figure.

Value

The hist method returns a list of length == 2, containing the arguments for the call to hist and the arguments to the call for legend, respectively.

Slots

time

A list containing 2 difftime objects (transform and fit), indicating the time elapsed for data transformation and for fitting the model to bootstrap data sets, respectively.

transData

Transformed data

bootDist

The vector of $chi^2$ values from bootstrap data sets fitted by the target model

origChi

The $chi^2$ value from the original data set

df

The degree of freedom of the model

bootP

The p value comparing the original $chi^2$ with the bootstrap distribution

Objects from the Class

Objects can be created via the bsBootMiss() function.

Author(s)

Terrence D. Jorgensen (University of Amsterdam; [email protected])

See Also

bsBootMiss()

Examples

# See the example from the bsBootMiss function

---

Bollen-Stine Bootstrap with the Existence of Missing Data

Description

Implement the Bollen and Stine's (1992) Bootstrap when missing observations exist. The implemented method is proposed by Savalei and Yuan (2009). This can be used in two ways. The first and easiest option is to fit the model to incomplete data in lavaan using the FIML estimator, then pass that lavaan object to bsBootMiss.

Usage

bsBootMiss(x, transformation = 2, nBoot = 500, model, rawData, Sigma, Mu,
  group, ChiSquared, EMcov, writeTransData = FALSE, transDataOnly = FALSE,
  writeBootData = FALSE, bootSamplesOnly = FALSE, writeArgs, seed = NULL,
  suppressWarn = TRUE, showProgress = TRUE, ...)

Arguments

x	A target lavaan object used in the Bollen-Stine bootstrap
transformation	The transformation methods in Savalei and Yuan (2009). There are three methods in the article, but only the first two are currently implemented here. Use transformation = 1 when there are few missing data patterns, each of which has a large size, such as in a planned-missing-data design. Use transformation = 2 when there are more missing data patterns. The currently unavailable transformation = 3 would be used when several missing data patterns have n = 1.
nBoot	The number of bootstrap samples.
model	Optional. The target model if x is not provided.
rawData	Optional. The target raw data set if x is not provided.
Sigma	Optional. The model-implied covariance matrix if x is not provided.
Mu	Optional. The model-implied mean vector if x is not provided.
group	Optional character string specifying the name of the grouping variable in rawData if x is not provided.
ChiSquared	Optional. The model's $\chi^2$ test statistic if x is not provided.
EMcov	Optional, if x is not provided. The EM (or Two-Stage ML) estimated covariance matrix used to speed up Transformation 2 algorithm.
writeTransData	Logical. If TRUE, the transformed data set is written to a text file, transDataOnly is set to TRUE, and the transformed data is returned invisibly.
transDataOnly	Logical. If TRUE, the result will provide the transformed data only.
writeBootData	Logical. If TRUE, the stacked bootstrap data sets are written to a text file, bootSamplesOnly is set to TRUE, and the list of bootstrap data sets are returned invisibly.
bootSamplesOnly	Logical. If TRUE, the result will provide bootstrap data sets only.
writeArgs	Optional list. If writeBootData = TRUE or writeBootData = TRUE, user can pass arguments to the utils::write.table() function as a list. Some default values are provided: file = "bootstrappedSamples.dat", row.names = FALSE, and na = "-999", but the user can override all of these by providing other values for those arguments in the writeArgs list.
seed	The seed number used in randomly drawing bootstrap samples.
suppressWarn	Logical. If TRUE, warnings from lavaan function will be suppressed when fitting the model to each bootstrap sample.
showProgress	Logical. Indicating whether to display a progress bar while fitting models to bootstrap samples.
...	The additional arguments in the lavaan::lavaan() function. See also lavaan::lavOptions()

Details

The second is designed for users of other software packages (e.g., LISREL, EQS, Amos, or Mplus). Users can import their data, $\chi^2$ value, and model-implied moments from another package, and they have the option of saving (or writing to a file) either the transformed data or bootstrapped samples of that data, which can be analyzed in other programs. In order to analyze the bootstrapped samples and return a p value, users of other programs must still specify their model using lavaan syntax.

Value

As a default, this function returns a BootMiss object containing the results of the bootstrap samples. Use show, summary, or hist to examine the results. Optionally, the transformed data set is returned if transDataOnly = TRUE. Optionally, the bootstrap data sets are returned if bootSamplesOnly = TRUE.

Author(s)

Terrence D. Jorgensen (University of Amsterdam; [email protected])

Syntax for transformations borrowed from http://www2.psych.ubc.ca/~vsavalei/

References

Bollen, K. A., & Stine, R. A. (1992). Bootstrapping goodness-of-fit measures in structural equation models. Sociological Methods & Research, 21(2), 205–229. doi:10.1177/0049124192021002004

Savalei, V., & Yuan, K.-H. (2009). On the model-based bootstrap with missing data: Obtaining a p-value for a test of exact fit. Multivariate Behavioral Research, 44(6), 741–763. doi:10.1080/00273170903333590

See Also

BootMiss

Examples

dat1 <- HolzingerSwineford1939
dat1$x5 <- ifelse(dat1$x1 <= quantile(dat1$x1, .3), NA, dat1$x5)
dat1$x9 <- ifelse(is.na(dat1$x5), NA, dat1$x9)

targetModel <- "
visual  =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed   =~ x7 + x8 + x9
"
targetFit <- sem(targetModel, dat1, meanstructure = TRUE, std.lv = TRUE,
                 missing = "fiml", group = "school")
summary(targetFit, fit = TRUE, standardized = TRUE)


## The number of bootstrap samples should be much higher than this example
temp <- bsBootMiss(targetFit, transformation = 1, nBoot = 10, seed = 31415)

temp
summary(temp)
hist(temp)
hist(temp, printLegend = FALSE) # suppress the legend
## user can specify alpha level (default: alpha = 0.05), and the number of
## digits to display (default: nd = 2).  Pass other arguments to hist(...),
## or a list of arguments to legend() via "legendArgs"
hist(temp, alpha = .01, nd = 3, xlab = "something else", breaks = 25,
     legendArgs = list("bottomleft", box.lty = 2))

---

Small-N correction for $chi^2$ test statistic

Description

Calculate small-N corrections for $chi^2$ model-fit test statistic to adjust for small sample size (relative to model size).

Usage

chisqSmallN(fit0, fit1 = NULL, smallN.method = if (is.null(fit1))
  c("swain", "yuan.2015") else "yuan.2005", ..., omit.imps = c("no.conv",
  "no.se"))

Arguments

fit0, fit1	lavaan::lavaan or lavaan.mi::lavaan.mi object(s)
smallN.method	character indicating the small-N correction method to use. Multiple may be chosen (all of which assume normality), as described in Shi et al. (2018): c("swain","yuan.2015","yuan.2005","bartlett"). Users may also simply select "all".
...	Additional arguments to the lavaan::lavTestLRT() or lavaan.mi::lavTestLRT.mi() functions. Ignored when is.null(fit1).
omit.imps	character vector specifying criteria for omitting imputations from pooled results. Ignored unless fit0 (and optionally fit1) is a lavaan.mi::lavaan.mi object. See lavaan.mi::lavTestLRT.mi() for a description of options and defaults.

Details

Four finite-sample adjustments to the chi-squared statistic are currently available, all of which are described in Shi et al. (2018). These all assume normally distributed data, and may not work well with severely nonnormal data. Deng et al. (2018, section 4) review proposed small-N adjustments that do not assume normality, which rarely show promise, so they are not implemented here. This function currently will apply small-N adjustments to scaled test statistics with a warning that they do not perform well (Deng et al., 2018).

Value

A list of numeric vectors: one for the originally requested statistic(s), along with one per requested smallN.method. All include the the (un)adjusted test statistic, its df, and the p value for the test under the null hypothesis that the model fits perfectly (or that the 2 models have equivalent fit). The adjusted chi-squared statistic(s) also include(s) the scaling factor for the small-N adjustment.

Author(s)

Terrence D. Jorgensen (University of Amsterdam; [email protected])

References

Deng, L., Yang, M., & Marcoulides, K. M. (2018). Structural equation modeling with many variables: A systematic review of issues and developments. Frontiers in Psychology, 9, 580. doi:10.3389/fpsyg.2018.00580

Shi, D., Lee, T., & Terry, R. A. (2018). Revisiting the model size effect in structural equation modeling. Structural Equation Modeling, 25(1), 21–40. doi:10.1080/10705511.2017.1369088

Examples

HS.model <- '
    visual  =~ x1 + b1*x2 + x3
    textual =~ x4 + b2*x5 + x6
    speed   =~ x7 + b3*x8 + x9
'
fit1 <- cfa(HS.model, data = HolzingerSwineford1939[1:50,])
## test a single model (implicitly compared to a saturated model)
chisqSmallN(fit1)

## fit a more constrained model
fit0 <- cfa(HS.model, data = HolzingerSwineford1939[1:50,],
            orthogonal = TRUE)
## compare 2 models
chisqSmallN(fit1, fit0)

---

Copy or save the result of lavaan or FitDiff objects into a clipboard or a file

Description

Copy or save the result of lavaan or FitDiff object into a clipboard or a file. From the clipboard, users may paste the result into the Microsoft Excel or spreadsheet application to create a table of the output.

Usage

clipboard(object, what = "summary", ...)

saveFile(object, file, what = "summary", tableFormat = FALSE,
  fit.measures = "default", writeArgs = list(), ...)

Arguments

object	An object of class lavaan::lavaan or FitDiff.
what	The attributes of the lavaan object to be copied in the clipboard. "summary" is to copy the screen provided from the summary function. "mifit" is to copy the result from the miPowerFit() function. Other attributes listed in the inspect method in the lavaan::lavaan could also be used, such as "coef", "se", "fit", "samp", and so on. Ignored for FitDiff-class objects.
...	Additional arguments when passing a lavaan object to the summary or miPowerFit() function.
file	A file name used for saving the result.
tableFormat	If TRUE, save the result in the table format using tabs for separation. Otherwise, save the result as the output screen printed in the R console.
fit.measures	character vector specifying names of fit measures returned by lavaan::fitMeasures() to be copied/saved. Only relevant if object is class FitDiff.
writeArgs	list of additional arguments to be passed to utils::write.table()

Value

The resulting output will be saved into a clipboard or a file. If using the clipboard function, users may paste it in the other applications.

Author(s)

Sunthud Pornprasertmanit ([email protected])

Terrence D. Jorgensen (University of Amsterdam; [email protected])

Examples

library(lavaan)
HW.model <- ' visual  =~ x1 + c1*x2 + x3
              textual =~ x4 + c1*x5 + x6
              speed   =~ x7 +    x8 + x9 '

fit <- cfa(HW.model, data = HolzingerSwineford1939, group = "school")

if(interactive()){
# Copy the summary of the lavaan object
clipboard(fit)

# pass additional arguments to summary() method for class?lavaan
clipboard(fit, rsquare = TRUE, standardized = TRUE, fit.measures = TRUE)

# Copy modification indices and fit stats from the miPowerFit() function
clipboard(fit, "mifit")

# Copy the parameter estimates
clipboard(fit, "coef")

# Copy the standard errors
clipboard(fit, "se")

# Copy the sample statistics
clipboard(fit, "samp")

# Copy the fit measures
clipboard(fit, "fit")

# Save the summary of the lavaan object
saveFile(fit, "out.txt")

# Save modification indices and fit stats from the miPowerFit() function
saveFile(fit, "out.txt", "mifit")

# Save the parameter estimates
saveFile(fit, "out.txt", "coef")

# Save the standard errors
saveFile(fit, "out.txt", "se")

# Save the sample statistics
saveFile(fit, "out.txt", "samp")

# Save the fit measures
saveFile(fit, "out.txt", "fit")
}

---

Combine the results from the quark function

Description

This function builds upon the quark() function to provide a final dataset comprised of the original dataset provided to quark() and enough principal components to be able to account for a certain level of variance in the data.

Usage

combinequark(quark, percent)

Arguments

quark	Provide the quark() object that was returned. It should be a list of objects. Make sure to include it in its entirety.
percent	Provide a percentage of variance that you would like to have explained. That many components (columns) will be extracted and kept with the output dataset. Enter this variable as a number WITHOUT a percentage sign.

Value

The output of this function is the original dataset used in quark combined with enough principal component scores to be able to account for the amount of variance that was requested.

Author(s)

Steven R. Chesnut (University of Southern Mississippi [email protected])

See Also

quark()

Examples

set.seed(123321)
dat <- HolzingerSwineford1939[,7:15]
misspat <- matrix(runif(nrow(dat) * 9) < 0.3, nrow(dat))
dat[misspat] <- NA
dat <- cbind(HolzingerSwineford1939[,1:3], dat)

quark.list <- quark(data = dat, id = c(1, 2))

final.data <- combinequark(quark = quark.list, percent = 80)

---

Build an object summarizing fit indices across multiple models

Description

This function will create the template to compare fit indices across multiple fitted lavaan objects. The results can be exported to a clipboard or a file later.

Usage

compareFit(..., nested = TRUE, argsLRT = list(), indices = TRUE,
  moreIndices = FALSE, baseline.model = NULL, nPrior = 1)

Arguments

...	fitted lavaan models or list(s) of lavaan objects. lavaan.mi::lavaan.mi objects are also accepted, but all models must belong to the same class.
nested	logical indicating whether the models in ... are nested. See net() for an empirical test of nesting.
argsLRT	list of arguments to pass to lavaan::lavTestLRT(), as well as to lavaan.mi::lavTestLRT.mi() and lavaan::fitMeasures() when comparing lavaan.mi::lavaan.mi models.
indices	logical indicating whether to return fit indices from the lavaan::fitMeasures() function. Selecting particular indices is controlled in the summary method; see FitDiff.
moreIndices	logical indicating whether to return fit indices from the moreFitIndices() function. Selecting particular indices is controlled in the summary method; see FitDiff.
baseline.model	optional fitted lavaan::lavaan model passed to lavaan::fitMeasures() to calculate incremental fit indices.
nPrior	passed to moreFitIndices(), if relevant

Value

A FitDiff object that saves model fit comparisons across multiple models. If the models are not nested, only fit indices for each model are returned. If the models are nested, the differences in fit indices are additionally returned, as well as test statistics comparing each sequential pair of models (ordered by their degrees of freedom).

Author(s)

Terrence D. Jorgensen (University of Amsterdam; [email protected])

Sunthud Pornprasertmanit ([email protected])

See Also

FitDiff, clipboard()

Examples

HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '

## non-nested models
fit1 <- cfa(HS.model, data = HolzingerSwineford1939)

m2 <- ' f1 =~ x1 + x2 + x3 + x4
        f2 =~ x5 + x6 + x7 + x8 + x9 '
fit2 <- cfa(m2, data = HolzingerSwineford1939)

(out1 <- compareFit(fit1, fit2, nested = FALSE))
summary(out1)


## nested model comparisons: measurement equivalence/invariance
fit.config <- cfa(HS.model, data = HolzingerSwineford1939, group = "school")
fit.metric <- cfa(HS.model, data = HolzingerSwineford1939, group = "school",
                  group.equal = "loadings")
fit.scalar <- cfa(HS.model, data = HolzingerSwineford1939, group = "school",
                  group.equal = c("loadings","intercepts"))
fit.strict <- cfa(HS.model, data = HolzingerSwineford1939, group = "school",
                  group.equal = c("loadings","intercepts","residuals"))

measEqOut <- compareFit(fit.config, fit.metric, fit.scalar, fit.strict,
                        moreIndices = TRUE) # include moreFitIndices()
summary(measEqOut)
summary(measEqOut, fit.measures = "all")
summary(measEqOut, fit.measures = c("aic", "bic", "sic", "ibic"))

---

Composite Reliability using SEM

Description

Calculate composite reliability from estimated factor-model parameters

Usage

compRelSEM(object, obs.var = TRUE, tau.eq = FALSE, ord.scale = TRUE,
  config = character(0), shared = character(0), higher = character(0),
  return.total = FALSE, dropSingle = TRUE, omit.factors = character(0),
  omit.indicators = character(0), omit.imps = c("no.conv", "no.se"),
  return.df = TRUE)

Arguments

object	A lavaan::lavaan or lavaan.mi::lavaan.mi object, expected to contain only exogenous common factors (i.e., a CFA model).
obs.var	logical indicating whether to compute reliability using observed variances in the denominator. Setting FALSE triggers using model-implied variances in the denominator.
tau.eq	logical indicating whether to assume (essential) tau-equivalence, yielding a coefficient analogous to $\alpha$. Setting FALSE yields an $\omega$-type coefficient.
ord.scale	logical indicating whether to apply Green and Yang's (2009, formula 21) correction, so that reliability is calculated for the actual ordinal response scale (ignored for factors with continuous indicators). Setting FALSE yields coefficients that are only applicable to the continuous latent-response scale.
config	character vector naming any configural constructs in a multilevel CFA. For these constructs (and optional total composite), Lai's (2021) coefficients $\omega^\textrm{W}$ and $\omega^\textrm{2L}$ are returned (or corresponding $\alpha$ coefficients when tau.eq=TRUE), rather than Geldhof et al.'s (2014) coefficients for hypothetical composites of latent components (although the same formula is used for $\omega^\textrm{W}$ in either case). Note that the same name must be used for the factor component represented at each level of the model.
shared	character vector naming any shared constructs in a multilevel CFA. For these constructs (and optional total composite), Lai's (2021) coefficient $\omega^\textrm{B}$ or $\alpha^\textrm{B}$ is returned, rather than Geldhof et al.'s (2014) between-level coefficient for hypothetical composites of latent cluster means. Lai's (2021) coefficient quantifies reliability relative to error associated with both indicators (measurement error) and subjects (sampling error), like a generalizability coefficient. Given that subjects can be considered as raters of their cluster's shared construct, an interrater reliability (IRR) coefficient is also returned, quantifying reliability relative to rater/sampling error alone. To quantify reliability relative to indicator/measurement error alone (i.e., $\omega^\textrm{2L}$), the ⁠shared=⁠ construct name(s) can additionally be included in ⁠config=⁠ argument.
higher	character vector naming any higher-order constructs in object for which composite reliability should be calculated. Ignored when tau.eq=TRUE because alpha is not based on a CFA model; instead, users must fit a CFA with tau-equivalence constraints. To obtain Lai's (2021) multilevel composite-reliability indices for a higher-order factor, do not use this argument; instead, specify the higher-order factor(s) using the ⁠shared=⁠ or ⁠config=⁠ argument (compRelSEM will automatically check whether it includes latent indicators and apply the appropriate formula).
return.total	logical indicating whether to return a final column containing the reliability of a composite of all indicators (not listed in omit.indicators) of first-order factors not listed in omit.factors. Ignored in 1-factor models, and should only be set TRUE if all factors represent scale dimensions that could be meaningfully collapsed to a single composite (scale sum or scale mean). Setting a negative value (e.g., -1 returns only the total-composite reliability (excluding coefficients per factor), except when requesting Lai's (2021) coefficients for multilevel configural or ⁠shared=⁠ constructs.
dropSingle	logical indicating whether to exclude factors defined by a single indicator from the returned results. If TRUE (default), single indicators will still be included in the total column when return.total = TRUE.
omit.factors	character vector naming any common factors modeled in object whose composite reliability is not of interest. For example, higher-order or method factors. Note that reliabilityL2() should be used to calculate composite reliability of a higher-order factor.
omit.indicators	character vector naming any observed variables that should be omitted from the composite whose reliability is calculated.
omit.imps	character vector specifying criteria for omitting imputations from pooled results. Can include any of c("no.conv", "no.se", "no.npd"), the first 2 of which are the default setting, which excludes any imputations that did not converge or for which standard errors could not be computed. The last option ("no.npd") would exclude any imputations which yielded a nonpositive definite covariance matrix for observed or latent variables, which would include any "improper solutions" such as Heywood cases. NPD solutions are not excluded by default because they are likely to occur due to sampling error, especially in small samples. However, gross model misspecification could also cause NPD solutions, users can compare pooled results with and without this setting as a sensitivity analysis to see whether some imputations warrant further investigation.
return.df	logical indicating whether to return reliability coefficients in a data.frame (one row per group/level), which is possible when every model block includes the same factors (after excluding those in omit.factors and applying dropSingle).

Details

Several coefficients for factor-analysis reliability have been termed "omega", which Cho (2021) argues is a misleading misnomer and argues for using $\rho$ to represent them all, differentiated by descriptive subscripts. In our package, we strive to provide unlabeled coefficients, leaving it to the user to decide on a label in their report. But we do use the symbols $\alpha$ and $\omega$ in the formulas below in order to distinguish coefficients that do (not) assume essential tau-equivalence. For higher-order constructs with latent indicators, only $\omega$ is available. Lai's (2021) multilevel coefficients are labeled in accordance with the symbols used in that article (more details below).

Bentler (1968) first introduced factor-analysis reliability for a unidimensional factor model with congeneric indicators, labeling the coeficients $\alpha$. McDonald (1999) later referred to this and other reliability coefficients, first as $\theta$ (in 1970), then as $\omega$, which is a source of confusion when reporting coefficients (Cho, 2021). Coefficients based on factor models were later generalized to account for multidimenisionality (possibly with cross-loadings) and correlated errors. The general $\omega$ formula implemented in this function is:

$\omega = \frac{\left( \sum^{k}_{i = 1} \lambda_i \right)^{2} Var\left( \psi \right)}{\bold{1}^\prime \hat{\Sigma} \bold{1}},$

where $\hat{\Sigma}$ can be the model-implied covariance matrix from either the saturated model (i.e., the "observed" covariance matrix, used by default) or from the hypothesized CFA model, controlled by the obs.var argument. A $k$-dimensional vector $\bold{1}$ is used to sum elements in the matrix. Note that if the model includes any directed effects (latent regression slopes), all coefficients are calculated from total factor variances: lavInspect(object, "cov.lv").

Assuming (essential) tau-equivalence (tau.eq=TRUE) makes $\omega$ equivalent to coefficient $\alpha$ from classical test theory (Cronbach, 1951):

$\alpha = \frac{k}{k - 1}\left[ 1 - \frac{\sum^{k}_{i = 1} \sigma_{ii}}{\sum^{k}_{i = 1} \sigma_{ii} + 2\sum_{i < j} \sigma_{ij}} \right],$

where $k$ is the number of items in a factor's composite, $\sigma_{ii}$ signifies item i's variance, and $\sigma_{ij}$ signifies the covariance between items i and j. Again, the obs.var argument controls whether $\alpha$ is calculated using the observed or model-implied covariance matrix.

By setting return.total=TRUE, one can estimate reliability for a single composite calculated using all indicators in a multidimensional CFA (Bentler, 1972, 2009). Setting return.total = -1 will return only the total-composite reliability (not per factor).

Higher-Order Factors: The reliability of a composite that represents a higher-order construct requires partitioning the model-implied factor covariance matrix $\Phi$ in order to isolate the common-factor variance associated only with the higher-order factor. Using a second-order factor model, the model-implied covariance matrix of observed indicators $\hat{\Sigma}$ can be partitioned into 3 sources:

1. the second-order common-factor (co)variance: $\Lambda \bold{B} \Phi_2 \bold{B}^{\prime} \Lambda^{\prime}$

2. the residual variance of the first-order common factors (i.e., not accounted for by the second-order factor): $\Lambda \Psi_{u} \Lambda^{\prime}$

3. the measurement error of observed indicators: $\Theta$

where $\Lambda$ contains first-order factor loadings, $\bold{B}$ contains second-order factor loadings, $\Phi_2$ is the model-implied covariance matrix of the second-order factor(s), and $\Psi_{u}$ is the covariance matrix of first-order factor disturbances. In practice, we can use the full $\bold{B}$ matrix and full model-implied $\Phi$ matrix (i.e., including all latent factors) because the zeros in $\bold{B}$ will cancel out unwanted components of $\Phi$. Thus, we can calculate the proportion of variance of a composite score calculated from the observed indicators (e.g., a total score or scale mean) that is attributable to the second-order factor (i.e., coefficient $\omega$):

$\omega=\frac{\bold{1}^{\prime} \Lambda \bold{B} \Phi \bold{B}^{\prime} \Lambda^{\prime} \bold{1} }{ \bold{1}^{\prime} \hat{\Sigma} \bold{1}},$

where $\bold{1}$ is the k-dimensional vector of 1s and k is the number of observed indicators in the composite. Note that if a higher-order factor also has observed indicators, it is necessary to model the observed indicators as single-indicator constructs, so that all of the higher-order factor indicators are latent (with loadings in the Beta matrix, not Lambda).

Categorical Indicators: When all indicators (per composite) are ordinal, the ord.scale argument controls whether the coefficient is calculated on the latent-response scale (FALSE) or on the observed ordinal scale (TRUE, the default). For $\omega$-type coefficients (tau.eq=FALSE), Green and Yang's (2009, formula 21) approach is used to transform factor-model results back to the ordinal response scale. When ord.scale=TRUE and tau.eq=TRUE, coefficient $\alpha$ is calculated using the covariance matrix calculated from the integer-valued numeric weights for ordinal categories, consistent with its definition (Chalmers, 2018) and the alpha function in the psych package; this implies obs.var=TRUE, so obs.var=FALSE will be ignored When ord.scale=FALSE, the standard $\alpha$ formula is applied to the polychoric correlation matrix ("ordinal $\alpha$"; Zumbo et al., 2007), estimated from the saturated or hypothesized model (see obs.var), and $\omega$ is calculated from CFA results without applying Green and Yang's (2009) correction (see Zumbo & Kroc, 2019, for a rationalization). No method analogous to Green and Yang (2009) has been proposed for calculating reliability with a mixture of categorical and continuous indicators, so an error is returned if object includes factors with a mixture of indicator types (unless omitted using omit.factors). If categorical indicators load on a different factor(s) than continuous indicators, then reliability will still be calculated separately for those factors, but return.total must be FALSE (unless omit.factors is used to isolate factors with indicators of the same type).

Multilevel Measurement Models: Under the default settings, compRelSEM() will apply the same formula in each "block" (group and/or level of analysis). In the case of multilevel (ML-)SEMs, this yields "reliability" for latent within- and between-level components, as proposed by Geldhof et al. (2014). Although this works fine to calculate reliability per group, this is not recommended for ML-SEMs because the coefficients do not correspond to actual composites that would be calculated from the observed data. Lai (2021) proposed coefficients for reliability of actual composites, depending on the type of construct, which requires specifying the names of constructs for which reliability is desired (or multiple constructs whose indicators would compose a multidimensional composite). Configural (⁠config=⁠) and/or ⁠shared=⁠ constructs can be specified; the same construct can be specified in both arguments, so that overall scale-reliability can be estimated for a shared construct by including it in config. Instead of organizing the output by block (the default), specifying ⁠config=⁠ and/or ⁠shared=⁠ will prompt organizing the list of output by ⁠$config⁠ and/or ⁠$shared⁠.

• The overall (⁠_2L⁠) scale reliability for configural constructs is returned, along with the reliability of a purely individual-level composite (⁠_W⁠, calculated by cluster-mean centering).

• The reliability for a shared construct quantifies generalizability across both indicators and raters (i.e., subjects rating their cluster's construct). Lüdtke et al. (2011) refer to these as measurement error and sampling error, respectively. An interrater reliability (IRR) coefficient is also returned, quantifying generalizability across rater/sampling-error only. To obtain a scale-reliability coefficient (quantifying a shared construct's generalizability across indicator/measurement-error only), include the same factor name in ⁠config=⁠. Jak et al. (2021) recommended modeling components of the same construct at both levels, but users may also saturate the within-level model (Lai, 2021).

Be careful about including Level-2 variables in the model, especially whether it makes sense to include them in a total composite for a Level-2 construct. dropSingle=TRUE only prevents estimating reliability for a single-indicator construct, not from including such an indicator in a total composite. It is permissible for ⁠shared=⁠ constructs to have additional indicators at Level-2 only. If it is necessary to model other Level-2 variables (e.g., to justify the missing-at-random assumption when using ⁠missing="FIML" estimation⁠), they should be placed in the ⁠omit.indicators=⁠ argument to exclude them from total composites.

Value

A numeric vector of composite reliability coefficients per factor, or a list of vectors per "block" (group and/or level of analysis), optionally returned as a data.frame when possible (see ⁠return.df=⁠ argument description for caveat). If there are multiple factors, whose multidimensional indicators combine into a single composite, users can request return.total=TRUE to add a column including a reliability coefficient for the total composite, or return.total = -1 to return only the total-composite reliability (ignored when ⁠config=⁠ or ⁠shared=⁠ is specified because each factor's specification must be checked across levels).

Author(s)

Terrence D. Jorgensen (University of Amsterdam; [email protected])

Uses hidden functions written by Sunthud Pornprasertmanit ([email protected]) for the old reliability() function.

References

Bentler, P. M. (1968). Alpha-maximized factor analysis (alphamax): Its relation to alpha and canonical factor analysis. Psychometrika, 33(3), 335–345. doi:10.1007/BF02289328

Bentler, P. M. (1972). A lower-bound method for the dimension-free measurement of internal consistency. Social Science Research, 1(4), 343–357. doi:10.1016/0049-089X(72)90082-8

Bentler, P. M. (2009). Alpha, dimension-free, and model-based internal consistency reliability. Psychometrika, 74(1), 137–143. doi:10.1007/s11336-008-9100-1

Chalmers, R. P. (2018). On misconceptions and the limited usefulness of ordinal alpha. Educational and Psychological Measurement, 78(6), 1056–1071. doi:10.1177/0013164417727036

Cho, E. (2021) Neither Cronbach’s alpha nor McDonald’s omega: A commentary on Sijtsma and Pfadt. Psychometrika, 86(4), 877–886. doi:10.1007/s11336-021-09801-1

Cronbach, L. J. (1951). Coefficient alpha and the internal structure of tests. Psychometrika, 16(3), 297–334. doi:10.1007/BF02310555

Geldhof, G. J., Preacher, K. J., & Zyphur, M. J. (2014). Reliability estimation in a multilevel confirmatory factor analysis framework. Psychological Methods, 19(1), 72–91. doi:10.1037/a0032138

Green, S. B., & Yang, Y. (2009). Reliability of summed item scores using structural equation modeling: An alternative to coefficient alpha. Psychometrika, 74(1), 155–167. doi:10.1007/s11336-008-9099-3

Jak, S., Jorgensen, T. D., & Rosseel, Y. (2021). Evaluating cluster-level factor models with lavaan and Mplus. Psych, 3(2), 134–152. doi:10.3390/psych3020012

Lai, M. H. C. (2021). Composite reliability of multilevel data: It’s about observed scores and construct meanings. Psychological Methods, 26(1), 90–102. doi:10.1037/met0000287

Lüdtke, O., Marsh, H. W., Robitzsch, A., & Trautwein, U. (2011). A 2 $\times$ 2 taxonomy of multilevel latent contextual models: Accuracy–bias trade-offs in full and partial error correction models. Psychological Methods, 16(4), 444–467. doi:10.1037/a0024376

McDonald, R. P. (1999). Test theory: A unified treatment. Mahwah, NJ: Erlbaum.

Zumbo, B. D., Gadermann, A. M., & Zeisser, C. (2007). Ordinal versions of coefficients alpha and theta for Likert rating scales. Journal of Modern Applied Statistical Methods, 6(1), 21–29. doi:10.22237/jmasm/1177992180

Zumbo, B. D., & Kroc, E. (2019). A measurement is a choice and Stevens’ scales of measurement do not help make it: A response to Chalmers. Educational and Psychological Measurement, 79(6), 1184–1197. doi:10.1177/0013164419844305

See Also

maximalRelia() for the maximal reliability of weighted composite

Examples

data(HolzingerSwineford1939)
HS9 <- HolzingerSwineford1939[ , c("x7","x8","x9")]
HSbinary <- as.data.frame( lapply(HS9, cut, 2, labels=FALSE) )
names(HSbinary) <- c("y7","y8","y9")
HS <- cbind(HolzingerSwineford1939, HSbinary)

HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ y7 + y8 + y9 '

fit <- cfa(HS.model, data = HS, ordered = c("y7","y8","y9"), std.lv = TRUE)

## works for factors with exclusively continuous OR categorical indicators
compRelSEM(fit)

## reliability for ALL indicators only available when they are
## all continuous or all categorical
compRelSEM(fit, omit.factors = "speed", return.total = TRUE)


## loop over visual indicators to calculate alpha if one indicator is removed
for (i in paste0("x", 1:3)) {
  cat("Drop ", i, ":\n", sep = "")
  print(compRelSEM(fit, omit.factors = c("textual","speed"),
                   omit.indicators = i, tau.eq = TRUE))
}
## item-total correlations obtainable by adding a composite to the data
HS$Visual <- HS$x1 + HS$x2 + HS$x3
cor(HS$Visual, y = HS[paste0("x", 1:3)])
## comparable to psych::alpha(HS[paste0("x", 1:3)])

## Reliability of a composite that represents a higher-order factor
mod.hi <- ' visual  =~ x1 + x2 + x3
            textual =~ x4 + x5 + x6
            speed   =~ x7 + x8 + x9
            general =~ visual + textual + speed '

fit.hi <- cfa(mod.hi, data = HolzingerSwineford1939)
compRelSEM(fit.hi, higher = "general")
## reliabilities for lower-order composites also returned


## works for multigroup models and for multilevel models (and both)
data(Demo.twolevel)
## assign clusters to arbitrary groups
Demo.twolevel$g <- ifelse(Demo.twolevel$cluster %% 2L, "type1", "type2")
model2 <- ' group: type1
  level: 1
    f1 =~ y1 + L2*y2 + L3*y3
    f2 =~ y4 + L5*y5 + L6*y6
  level: 2
    f1 =~ y1 + L2*y2 + L3*y3
    f2 =~ y4 + L5*y5 + L6*y6

group: type2
  level: 1
    f1 =~ y1 + L2*y2 + L3*y3
    f2 =~ y4 + L5*y5 + L6*y6
  level: 2
    f1 =~ y1 + L2*y2 + L3*y3
    f2 =~ y4 + L5*y5 + L6*y6
'
fit2 <- sem(model2, data = Demo.twolevel, cluster = "cluster", group = "g")
compRelSEM(fit2) # Geldhof's indices (hypothetical, for latent components)

## Lai's (2021) indices for Level-1 and configural constructs
compRelSEM(fit2, config = c("f1","f2"))
## Lai's (2021) indices for shared (Level-2) constructs
## (also an interrater reliability coefficient)
compRelSEM(fit2, shared = c("f1","f2"))


## Shared construct using saturated within-level model
mod.sat1 <- ' level: 1
  y1 ~~ y1 + y2 + y3 + y4 + y5 + y6
  y2 ~~ y2 + y3 + y4 + y5 + y6
  y3 ~~ y3 + y4 + y5 + y6
  y4 ~~ y4 + y5 + y6
  y5 ~~ y5 + y6
  y6 ~~ y6

  level: 2
  f1 =~ y1 + L2*y2 + L3*y3
  f2 =~ y4 + L5*y5 + L6*y6
'
fit.sat1 <- sem(mod.sat1, data = Demo.twolevel, cluster = "cluster")
compRelSEM(fit.sat1, shared = c("f1","f2"))


## Simultaneous shared-and-configural model (Stapleton et al, 2016, 2019),
## not recommended, but possible by omitting shared or configural factor.
mod.both <- ' level: 1
    fc =~ y1 + L2*y2 + L3*y3 + L4*y4 + L5*y5 + L6*y6
  level: 2
  ## configural construct
    fc =~ y1 + L2*y2 + L3*y3 + L4*y4 + L5*y5 + L6*y6
  ## orthogonal shared construct
    fs =~ NA*y1 + y2 + y3 + y4 + y5 + y6
    fs ~~ 1*fs + 0*fc
'
fit.both <- sem(mod.both, data = Demo.twolevel, cluster = "cluster")
compRelSEM(fit.both, shared = "fs", config = "fc")

---

Simulated Dataset to Demonstrate Two-way Latent Interaction

Description

A simulated data set with 2 independent factors and 1 dependent factor where each factor has three indicators

Usage

dat2way

Format

A data.frame with 500 observations of 9 variables.

x1

The first indicator of the first independent factor

x2

The second indicator of the first independent factor

x3

The third indicator of the first independent factor

x4

The first indicator of the second independent factor

x5

The second indicator of the second independent factor

x6

The third indicator of the second independent factor

x7

The first indicator of the dependent factor

x8

The second indicator of the dependent factor

x9

The third indicator of the dependent factor

Source

Data were generated by the MASS::mvrnorm() function in the MASS package.

Examples

head(dat2way)

---

Simulated Dataset to Demonstrate Three-way Latent Interaction

Description

A simulated data set with 3 independent factors and 1 dependent factor where each factor has three indicators

Usage

dat3way

Format

A data.frame with 500 observations of 12 variables.

x1

The first indicator of the first independent factor

x2

The second indicator of the first independent factor

x3

The third indicator of the first independent factor

x4

The first indicator of the second independent factor

x5

The second indicator of the second independent factor

x6

The third indicator of the second independent factor

x7

The first indicator of the third independent factor

x8

The second indicator of the third independent factor

x9

The third indicator of the third independent factor

x10

The first indicator of the dependent factor

x11

The second indicator of the dependent factor

x12

The third indicator of the dependent factor

Source

Data were generated by the MASS::mvrnorm() function in the MASS package.

Examples

head(dat3way)

---

Simulated Data set to Demonstrate Categorical Measurement Invariance

Description

A simulated data set with 2 factors with 4 indicators each separated into two groups

Usage

datCat

Format

A data.frame with 200 observations of 9 variables.

g

Sex of respondents

u1

Indicator 1

u2

Indicator 2

u3

Indicator 3

u4

Indicator 4

u5

Indicator 5

u6

Indicator 6

u7

Indicator 7

u8

Indicator 8

Source

Data were generated using the lavaan package.

Examples

head(datCat)

---

Calculate discriminant validity statistics

Description

Calculate discriminant validity statistics based on a fitted lavaan object

Usage

discriminantValidity(object, cutoff = 0.9, merge = FALSE, level = 0.95,
  boot.ci.type = "perc")

Arguments

object	The lavaan::lavaan model object returned by the lavaan::cfa() function.
cutoff	A cutoff to be used in the constrained models in likelihood ratio tests.
merge	Whether the constrained models should be constructed by merging two factors as one. Implies cutoff = 1.
level	The confidence level required.
boot.ci.type	If bootstrapping was used, the type of interval required. The value should be one of "norm", "basic", "perc", or "bca.simple". For the first three options, see the help page of the boot.ci function in the boot package. The "bca.simple" option produces intervals using the adjusted bootstrap percentile (BCa) method, but with no correction for acceleration (only for bias). Note that the p-value is still computed assuming that the z-statistic follows a standard normal distribution.

Details

Evaluated on the measurement scale level, discriminant validity is commonly evaluated by checking if each pair of latent correlations is sufficiently below one (in absolute value) that the latent variables can be thought of representing two distinct constructs.

discriminantValidity function calculates two sets of statistics that are commonly used in discriminant validity evaluation. The first set are factor correlation estimates and their confidence intervals. The second set is a series of nested model tests, where the baseline model is compared against a set of constrained models that are constructed by constraining each factor correlation to the specified cutoff one at a time.

The function assume that the object is set of confirmatory factor analysis results where the latent variables are scaled by fixing their variances to 1s. If the model is not a CFA model, the function will calculate the statistics for the correlations among exogenous latent variables, but for the residual variances with endogenous variables. If the latent variables are scaled in some other way (e.g. fixing the first loadings), the function issues a warning and re-estimates the model by fixing latent variances to 1 (and estimating all loadings) so that factor covariances are already estimated as correlations.

The likelihood ratio tests are done by comparing the original baseline model against more constrained alternatives. By default, these alternatives are constructed by fixing each correlation at a time to a cutoff value. The typical purpose of this test is to demonstrate that the estimated factor correlation is well below the cutoff and a significant $chi^2$ statistic thus indicates support for discriminant validity. In some cases, the original correlation estimate may already be greater than the cutoff, making it redundant to fit a "restricted" model. When this happens, the likelihood ratio test will be replaced by comparing the baseline model against itself. For correlations that are estimated to be negative, a negation of the cutoff is used in the constrained model.

Another alternative is to do a nested model comparison against a model where two factors are merged as one by setting the merge argument to TRUE. In this comparison, the constrained model is constructed by removing one of the correlated factors from the model and assigning its indicators to the factor that remains in the model.

Value

A data.frame of latent variable correlation estimates, their confidence intervals, and a likelihood ratio tests against constrained models. with the following attributes:

baseline

The baseline model after possible rescaling.

constrained

A list of the fitted constrained models used in the likelihood ratio test.

Author(s)

Mikko Rönkkö (University of Jyväskylä; [email protected]):

References

Rönkkö, M., & Cho, E. (2022). An updated guideline for assessing discriminant validity. Organizational Research Methods, 25(1), 6–14. doi:10.1177/1094428120968614

Examples

library(lavaan)

HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '

fit <- cfa(HS.model, data = HolzingerSwineford1939)
discriminantValidity(fit)
discriminantValidity(fit, merge = TRUE)

---

Class For Rotated Results from EFA

Description

This class contains the results of rotated exploratory factor analysis

Usage

## S4 method for signature 'EFA'
show(object)

## S4 method for signature 'EFA'
summary(object, suppress = 0.1, sort = TRUE)

Arguments

object	object of class EFA
suppress	any standardized loadings less than the specified value will not be printed to the screen
sort	logical. If TRUE (default), factor loadings will be sorted by their size in the console output

Slots

loading

Rotated standardized factor loading matrix

rotate

Rotation matrix

gradRotate

gradient of the objective function at the rotated loadings

convergence

Convergence status

phi:

Factor correlation matrix. Will be an identity matrix if orthogonal rotation is used.

se

Standard errors of the rotated standardized factor loading matrix

method

Method of rotation

call

The command used to generate this object

Objects from the Class

Objects can be created via the orthRotate or oblqRotate function.

Author(s)

Sunthud Pornprasertmanit ([email protected])

See Also

efaUnrotate; orthRotate; oblqRotate

Examples

unrotated <- efaUnrotate(HolzingerSwineford1939, nf = 3,
                         varList = paste0("x", 1:9), estimator = "mlr")
summary(unrotated, std = TRUE)
lavInspect(unrotated, "std")

# Rotated by Quartimin
rotated <- oblqRotate(unrotated, method = "quartimin")
summary(rotated)

---

Empirical Kaiser criterion

Description

Identify the number of factors to extract based on the Empirical Kaiser Criterion (EKC). The analysis can be run on a data.frame or data matrix (data), or on a correlation or covariance matrix (sample.cov) and the sample size (sample.nobs). A data.frame is returned with two columns: the eigenvalues from your data or covariance matrix and the reference eigenvalues. The number of factors suggested by the Empirical Kaiser Criterion (i.e. the sample eigenvalues greater than the reference eigenvalues), and the number of factors suggested by the original Kaiser Criterion (i.e. sample eigenvalues > 1) is printed above the output.

Usage

efa.ekc(data = NULL, sample.cov = NULL, sample.nobs = NULL,
  missing = "default", ordered = NULL, plot = TRUE)

Arguments

data	A data.frame or data matrix containing columns of variables to be factor-analyzed.
sample.cov	A covariance or correlation matrix can be used, instead of data, to estimate the eigenvalues.
sample.nobs	Number of observations (i.e. sample size) if is.null(data) and ⁠sample.cov=⁠ is used.
missing	If "listwise", incomplete cases are removed listwise from the data.frame. If "direct" or "ml" or "fiml" and the ⁠estimator=⁠ is maximum likelihood, an EM algorithm is used to estimate an unrestricted covariance matrix (and mean vector). If "pairwise", pairwise deletion is used. If '"default"“, the value is set depending on the estimator and the mimic option (see lavaan::lavCor() for details).
ordered	character vector. Only used if object is a data.frame. Treat these variables as ⁠ordered=⁠ (ordinal) variables. Importantly, all other variables will be treated as numeric (unless is.ordered == TRUE in data). (see also lavCor)
plot	logical. Whether to print a scree plot comparing the sample eigenvalues with the reference eigenvalues.

Value

A data.frame showing the sample and reference eigenvalues.

The number of factors suggested by the Empirical Kaiser Criterion (i.e. the sample eigenvalues greater than the reference eigenvalues) is returned as an attribute (see Examples).

The number of factors suggested by the original Kaiser Criterion (i.e. sample eigenvalues > 1) is also printed as a header to the data.frame

Author(s)

Ylenio Longo (University of Nottingham; [email protected])

Terrence D. Jorgensen (University of Amsterdam; [email protected])

References

Braeken, J., & van Assen, M. A. L. M. (2017). An empirical Kaiser criterion. Psychological Methods, 22(3), 450–466. doi:10.1037/met0000074

Examples

## Simulate data with 3 factors
model <- '
  f1 =~ .3*x1 + .5*x2 + .4*x3
  f2 =~ .3*x4 + .5*x5 + .4*x6
  f3 =~ .3*x7 + .5*x8 + .4*x9
'
dat <- simulateData(model, seed = 123)
## save summary statistics
myCovMat <- cov(dat)
myCorMat <- cor(dat)
N <- nrow(dat)

## Run the EKC function
(out <- efa.ekc(dat))

## To extract the recommended number of factors using the EKC:
attr(out, "nfactors")

## If you do not have raw data, you can use summary statistics
(x1 <- efa.ekc(sample.cov = myCovMat, sample.nobs = N, plot = FALSE))
(x2 <- efa.ekc(sample.cov = myCorMat, sample.nobs = N, plot = FALSE))

---

Simulated Data set to Demonstrate Longitudinal Measurement Invariance

Description

A simulated data set with 1 factors with 3 indicators in three timepoints

Usage

exLong

Format

A data.frame with 200 observations of 10 variables.

sex

Sex of respondents

y1t1

Indicator 1 in Time 1

y2t1

Indicator 2 in Time 1

y3t1

Indicator 3 in Time 1

y1t2

Indicator 1 in Time 2

y2t2

Indicator 2 in Time 2

y3t2

Indicator 3 in Time 2

y1t3

Indicator 1 in Time 3

y2t3

Indicator 2 in Time 3

y3t3

Indicator 3 in Time 3

Source

Data were generated using the simsem package.

Examples

head(exLong)

---

Find the statistical power based on population RMSEA

Description

Find the proportion of the samples from the sampling distribution of RMSEA in the alternative hypothesis rejected by the cutoff dervied from the sampling distribution of RMSEA in the null hypothesis. This function can be applied for both test of close fit and test of not-close fit (MacCallum, Browne, & Suguwara, 1996)

Usage

findRMSEApower(rmsea0, rmseaA, df, n, alpha = 0.05, group = 1)

Arguments

rmsea0	Null RMSEA
rmseaA	Alternative RMSEA
df	Model degrees of freedom
n	Sample size of a dataset
alpha	Alpha level used in power calculations
group	The number of group that is used to calculate RMSEA.

Details

This function find the proportion of sampling distribution derived from the alternative RMSEA that is in the critical region derived from the sampling distribution of the null RMSEA. If rmseaA is greater than rmsea0, the test of close fit is used and the critical region is in the right hand side of the null sampling distribution. On the other hand, if rmseaA is less than rmsea0, the test of not-close fit is used and the critical region is in the left hand side of the null sampling distribution (MacCallum, Browne, & Suguwara, 1996).

There is also a Shiny app called "power4SEM" that provides a graphical user interface for this functionality (Jak et al., in press). It can be accessed at https://sjak.shinyapps.io/power4SEM/.

Author(s)

Sunthud Pornprasertmanit ([email protected])

References

MacCallum, R. C., Browne, M. W., & Sugawara, H. M. (1996). Power analysis and determination of sample size for covariance structure modeling. Psychological Methods, 1(2), 130–149. doi:10.1037/1082-989X.1.2.130

Jak, S., Jorgensen, T. D., Verdam, M. G., Oort, F. J., & Elffers, L. (2021). Analytical power calculations for structural equation modeling: A tutorial and Shiny app. Behavior Research Methods, 53, 1385–1406. doi:10.3758/s13428-020-01479-0

See Also

• plotRMSEApower() to plot the statistical power based on population RMSEA given the sample size

• plotRMSEAdist() to visualize the RMSEA distributions

• findRMSEAsamplesize() to find the minium sample size for a given statistical power based on population RMSEA

Examples

findRMSEApower(rmsea0 = .05, rmseaA = .08, df = 20, n = 200)

---

Find power given a sample size in nested model comparison

Description

Find the sample size that the power in rejection the samples from the alternative pair of RMSEA is just over the specified power.

Usage

findRMSEApowernested(rmsea0A = NULL, rmsea0B = NULL, rmsea1A,
  rmsea1B = NULL, dfA, dfB, n, alpha = 0.05, group = 1)

Arguments

rmsea0A	The $H_0$ baseline RMSEA
rmsea0B	The $H_0$ alternative RMSEA (trivial misfit)
rmsea1A	The $H_1$ baseline RMSEA
rmsea1B	The $H_1$ alternative RMSEA (target misfit to be rejected)
dfA	degree of freedom of the more-restricted model
dfB	degree of freedom of the less-restricted model
n	Sample size
alpha	The alpha level
group	The number of group in calculating RMSEA

Author(s)

Bell Clinton

Pavel Panko (Texas Tech University; [email protected])

Sunthud Pornprasertmanit ([email protected])

References

MacCallum, R. C., Browne, M. W., & Cai, L. (2006). Testing differences between nested covariance structure models: Power analysis and null hypotheses. Psychological Methods, 11(1), 19–35. doi:10.1037/1082-989X.11.1.19

See Also

• plotRMSEApowernested() to plot the statistical power for nested model comparison based on population RMSEA given the sample size

• findRMSEAsamplesizenested() to find the minium sample size for a given statistical power in nested model comparison based on population RMSEA

Examples

findRMSEApowernested(rmsea0A = 0.06, rmsea0B = 0.05, rmsea1A = 0.08,
                     rmsea1B = 0.05, dfA = 22, dfB = 20, n = 200,
                     alpha = 0.05, group = 1)

---

Find the minimum sample size for a given statistical power based on population RMSEA

Description

Find the minimum sample size for a specified statistical power based on population RMSEA. This function can be applied for both test of close fit and test of not-close fit (MacCallum, Browne, & Suguwara, 1996)

Usage

findRMSEAsamplesize(rmsea0, rmseaA, df, power = 0.8, alpha = 0.05,
  group = 1)

Arguments

rmsea0	Null RMSEA
rmseaA	Alternative RMSEA
df	Model degrees of freedom
power	Desired statistical power to reject misspecified model (test of close fit) or retain good model (test of not-close fit)
alpha	Alpha level used in power calculations
group	The number of group that is used to calculate RMSEA.

Details

This function find the minimum sample size for a specified power based on an iterative routine. The sample size keep increasing until the calculated power from findRMSEApower() function is just over the specified power. If group is greater than 1, the resulting sample size is the sample size per group.

Author(s)

Sunthud Pornprasertmanit ([email protected])

References

MacCallum, R. C., Browne, M. W., & Sugawara, H. M. (1996). Power analysis and determination of sample size for covariance structure modeling. Psychological Methods, 1(2), 130–149. doi:10.1037/1082-989X.1.2.130

Jak, S., Jorgensen, T. D., Verdam, M. G., Oort, F. J., & Elffers, L. (2021). Analytical power calculations for structural equation modeling: A tutorial and Shiny app. Behavior Research Methods, 53, 1385–1406. doi:10.3758/s13428-020-01479-0

See Also

• plotRMSEApower() to plot the statistical power based on population RMSEA given the sample size

• plotRMSEAdist() to visualize the RMSEA distributions

• findRMSEApower() to find the statistical power based on population RMSEA given a sample size

Examples

findRMSEAsamplesize(rmsea0 = .05, rmseaA = .08, df = 20, power = 0.80)

---

Find sample size given a power in nested model comparison

Description

Find the sample size that the power in rejection the samples from the alternative pair of RMSEA is just over the specified power.

Usage

findRMSEAsamplesizenested(rmsea0A = NULL, rmsea0B = NULL, rmsea1A,
  rmsea1B = NULL, dfA, dfB, power = 0.8, alpha = 0.05, group = 1)

Arguments

rmsea0A	The $H_0$ baseline RMSEA
rmsea0B	The $H_0$ alternative RMSEA (trivial misfit)
rmsea1A	The $H_1$ baseline RMSEA
rmsea1B	The $H_1$ alternative RMSEA (target misfit to be rejected)
dfA	degree of freedom of the more-restricted model.
dfB	degree of freedom of the less-restricted model.
power	The desired statistical power.
alpha	The alpha level.
group	The number of group in calculating RMSEA.

Author(s)

Bell Clinton

Pavel Panko (Texas Tech University; [email protected])

Sunthud Pornprasertmanit ([email protected])

References

MacCallum, R. C., Browne, M. W., & Cai, L. (2006). Testing differences between nested covariance structure models: Power analysis and null hypotheses. Psychological Methods, 11(1), 19–35. doi:10.1037/1082-989X.11.1.19

See Also

• plotRMSEApowernested() to plot the statistical power for nested model comparison based on population RMSEA given the sample size

• findRMSEApowernested() to find the power for a given sample size in nested model comparison based on population RMSEA

Examples

findRMSEAsamplesizenested(rmsea0A = 0, rmsea0B = 0, rmsea1A = 0.06,
                          rmsea1B = 0.05, dfA = 22, dfB = 20, power = 0.80,
                          alpha = .05, group = 1)

---

Class For Representing A Template of Model Fit Comparisons

Description

This class contains model fit measures and model fit comparisons among multiple models

Usage

## S4 method for signature 'FitDiff'
show(object)

## S4 method for signature 'FitDiff'
summary(object, fit.measures = "default", nd = 3,
  tag = "†")

Arguments

object	object of class FitDiff
fit.measures	character vector naming fit indices the user can request from lavaan::fitMeasures(). If "default", the fit measures will be c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr", "aic", "bic"). If "all", all available fit measures will be returned.
nd	number of digits printed
tag	single character used to flag the model preferred by each fit index. To omit tags, set to NULL or NA.

Slots

name

character. The name of each model

model.class

character. One class to which each model belongs

nested

data.frame. Model fit comparisons between adjacently nested models that are ordered by their degrees of freedom (df)

fit

data.frame. Fit measures of all models specified in the name slot, ordered by their df

fit.diff

data.frame. Sequential differences in fit measures in the fit slot

Objects from the Class

Objects can be created via the compareFit() function.

Author(s)

Terrence D. Jorgensen (University of Amsterdam; [email protected])

Sunthud Pornprasertmanit ([email protected])

See Also

compareFit(); clipboard()

Examples

HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '
fit.config <- cfa(HS.model, data = HolzingerSwineford1939, group = "school")
## invariance constraints
fit.metric <- cfa(HS.model, data = HolzingerSwineford1939, group = "school",
                  group.equal = "loadings")
fit.scalar <- cfa(HS.model, data = HolzingerSwineford1939, group = "school",
                  group.equal = c("loadings","intercepts"))
fit.strict <- cfa(HS.model, data = HolzingerSwineford1939, group = "school",
                  group.equal = c("loadings","intercepts","residuals"))
measEqOut <- compareFit(fit.config, fit.metric, fit.scalar, fit.strict)
summary(measEqOut)
summary(measEqOut, fit.measures = "all")
summary(measEqOut, fit.measures = c("aic", "bic"))

if(interactive()){
## Save results to a file
saveFile(measEqOut, file = "measEq.txt")

## Copy to a clipboard
clipboard(measEqOut)
}

---

Fraction of Missing Information.

Description

This function estimates the Fraction of Missing Information (FMI) for summary statistics of each variable, using either an incomplete data set or a list of imputed data sets.

Usage

fmi(data, method = "saturated", group = NULL, ords = NULL,
  varnames = NULL, exclude = NULL, return.fit = FALSE)

Arguments

data	Either a single data.frame with incomplete observations, or a list of imputed data sets.
method	character. If "saturated" or "sat" (default), the model used to estimate FMI is a freely estimated covariance matrix and mean vector for numeric variables, and/or polychoric correlations and thresholds for ordered categorical variables, for each group (if applicable). If "null", only means and variances are estimated for numeric variables, and/or thresholds for ordered categorical variables (i.e., covariances and/or polychoric/polyserial correlations are constrained to zero). See Details for more information.
group	character. The optional name of a grouping variable, to request FMI in each group.
ords	Optional character vector naming ordered-categorical variables, if they are not already stored as class ordered in data.
varnames	Optional character vector of variable names, to calculate FMI for a subset of variables in data. By default, all numeric and ⁠ordered=⁠ variables will be included, unless ⁠data=⁠ is a single incomplete data.frame, in which case only numeric variables can be used with FIML estimation. Other variable types will be removed.
exclude	Optional character vector naming variables to exclude from the analysis.
return.fit	logical. If TRUE, the fitted lavaan::lavaan or lavaan.mi::lavaan.mi model is returned, so FMI can be found from summary(..., fmi=TRUE).

Details

The function estimates a saturated model with lavaan::lavaan() for a single incomplete data set using FIML, or with lavaan.mi::lavaan.mi() for a list of imputed data sets. If method = "saturated", FMI will be estiamted for all summary statistics, which could take a lot of time with big data sets. If method = "null", FMI will only be estimated for univariate statistics (e.g., means, variances, thresholds). The saturated model gives more reliable estimates, so it could also help to request a subset of variables from a large data set.

Value

fmi() returns a list with at least 2 of the following:

Covariances	A list of symmetric matrices: (1) the estimated/pooled covariance matrix, or a list of group-specific matrices (if applicable) and (2) a matrix of FMI, or a list of group-specific matrices (if applicable). Only available if method = "saturated". When method="cor", this element is replaced by Correlations.
Variances	The estimated/pooled variance for each numeric variable. Only available if method = "null" (otherwise, it is on the diagonal of Covariances).
Means	The estimated/pooled mean for each numeric variable.
Thresholds	The estimated/pooled threshold(s) for each ordered-categorical variable.

Author(s)

Mauricio Garnier Villarreal (Vrije Universiteit Amsterdam; [email protected])

Terrence Jorgensen (University of Amsterdam; [email protected])

References

Rubin, D. B. (1987). Multiple imputation for nonresponse in surveys. New York, NY: Wiley.

Savalei, V. & Rhemtulla, M. (2012). On obtaining estimates of the fraction of missing information from full information maximum likelihood. Structural Equation Modeling, 19(3), 477–494. doi:10.1080/10705511.2012.687669

Wagner, J. (2010). The fraction of missing information as a tool for monitoring the quality of survey data. Public Opinion Quarterly, 74(2), 223–243. doi:10.1093/poq/nfq007

Examples

HSMiss <- HolzingerSwineford1939[ , c(paste("x", 1:9, sep = ""),
                                      "ageyr","agemo","school")]
set.seed(12345)
HSMiss$x5 <- ifelse(HSMiss$x5 <= quantile(HSMiss$x5, .3), NA, HSMiss$x5)
age <- HSMiss$ageyr + HSMiss$agemo/12
HSMiss$x9 <- ifelse(age <= quantile(age, .3), NA, HSMiss$x9)

## calculate FMI (using FIML, provide partially observed data set)
(out1 <- fmi(HSMiss, exclude = "school"))
(out2 <- fmi(HSMiss, exclude = "school", method = "null"))
(out3 <- fmi(HSMiss, varnames = c("x5","x6","x7","x8","x9")))
(out4 <- fmi(HSMiss, method = "cor", group = "school")) # correlations by group

## significance tests in lavaan(.mi) object
out5 <- fmi(HSMiss, method = "cor", return.fit = TRUE)
summary(out5) # factor loading == SD, covariance = correlation

if(requireNamespace("lavaan.mi")){
  ## ordered-categorical data
  data(binHS5imps, package = "lavaan.mi")

  ## calculate FMI, using list of imputed data sets
  fmi(binHS5imps, group = "school")
}

---

Wrapper for goric.lavaan() from the restriktor package

Description

The goricaSEM() function is an interface to restriktor::goric.lavaan(), allowing users to perform generalized order-restricted information criterion approximation (GORICA) analysis specifically for structural equation models fitted using the lavaan package.

Usage

goricaSEM(object, ..., hypotheses = NULL, comparison = NULL,
  type = "gorica", standardized = FALSE, debug = FALSE)

Arguments

object	A lavaan::lavaan object.
...	Additional arguments passed to restriktor::goric.lavaan().
hypotheses	A named list of hypotheses to test. See Details for information on how to specify hypotheses.
comparison	A character string specifying the type of comparison. Options are "unconstrained", "complement", or "none". Default behavior depends on the number of hypotheses.
type	A character string indicating the type of analysis, either "gorica" (default) or "goricac".
standardized	logical indicating whether standardized estimates are used in the analysis. Defaults to FALSE.
debug	logical indicating whether to print debugging information. Defaults to FALSE.

Details

This function is designed as a wrapper for the restriktor::goric.lavaan() function. It calculates GORICA values and weights, which can be used to compare models or hypotheses under inequality constraints.

The ⁠hypotheses=⁠ argument allows users to specify constraints in text-based syntax or matrix notation. For text-based syntax, constraints are specified as a string (e.g., "a1 > a2"). For matrix notation, a named list with ⁠$constraints⁠, ⁠$rhs⁠, and ⁠$neq⁠ elements can be provided.

The ⁠comparison=⁠ argument determines whether the specified hypothesis is compared against its "complement", the "unconstrained" model, or neither ("none").

Value

A list containing the results of the goric.lavaan function, including:

• The log-likelihood.

• Penalty term.

• GORIC(A) values and weights.

• Relative GORIC(A) weights.

Author(s)

Leonard Vanbrabant and Rebecca Kuiper

References

Kuiper, R. M., Hoijtink, H., & Silvapulle, M. J. (2011). An Akaike-type information criterion for model selection under inequality constraints. Biometrika, 98(2), 495–501. doi:10.1093/biomet/asr002

Vanbrabant, L., Van Loey, N., & Kuiper, R. M. (2020). Evaluating a theory-based hypothesis against its complement using an AIC-type information criterion with an application to facial burn injury. Psychological Methods, 25(2), 129–142. doi:10.1037/met0000238

See Also

restriktor::goric.lavaan()

Examples

## Example: Perform GORICA analysis on a lavaan model
library(lavaan)
library(restriktor)

## Define the SEM model
model <- '
  ind60 =~ x1 + x2 + x3
  dem60 =~ y1 + a1*y2 + b1*y3 + c1*y4
  dem65 =~ y5 + a2*y6 + b2*y7 + c2*y8
  dem60 ~ ind60
  dem65 ~ ind60 + dem60
  y1 ~~ y5
  y2 ~~ y4 + y6
  y3 ~~ y7
  y4 ~~ y8
  y6 ~~ y8
'

## Fit the model
data(PoliticalDemocracy)
fit <- sem(model, data = PoliticalDemocracy)

## Define hypotheses
myHypothesis <- 'a1 > a2, b1 > b2, c1 > c2'

## Perform GORICA analysis
result <- goricaSEM(fit, hypotheses = list(H1 = myHypothesis),
                    standardized = FALSE, comparison = "complement",
                    type = "gorica")

## Print result
print(result)

---

Assessing Discriminant Validity using Heterotrait–Monotrait Ratio

Description

This function assesses discriminant validity through the heterotrait-monotrait ratio (HTMT) of the correlations (Henseler, Ringlet & Sarstedt, 2015). Specifically, it assesses the arithmetic (Henseler et al., ) or geometric (Roemer et al., 2021) mean correlation among indicators across constructs (i.e. heterotrait–heteromethod correlations) relative to the geometric-mean correlation among indicators within the same construct (i.e. monotrait–heteromethod correlations). The resulting HTMT(2) values are interpreted as estimates of inter-construct correlations. Absolute values of the correlations are recommended to calculate the HTMT matrix, and are required to calculate HTMT2. Correlations are estimated using the lavaan::lavCor() function.

Usage

htmt(model, data = NULL, sample.cov = NULL, missing = "listwise",
  ordered = NULL, absolute = TRUE, htmt2 = TRUE)

Arguments

model	lavaan lavaan::model.syntax() of a confirmatory factor analysis model where at least two factors are required for indicators measuring the same construct.
data	A data.frame or data matrix
sample.cov	A covariance or correlation matrix can be used, instead of ⁠data=⁠, to estimate the HTMT values.
missing	If "listwise", cases with missing values are removed listwise from the data frame. If "direct" or "ml" or "fiml" and the estimator is maximum likelihood, an EM algorithm is used to estimate the unrestricted covariance matrix (and mean vector). If "pairwise", pairwise deletion is used. If "default", the value is set depending on the estimator and the mimic option (see details in lavaan::lavCor()).
ordered	Character vector. Only used if object is a data.frame. Treat these variables as ordered (ordinal) variables. Importantly, all other variables will be treated as numeric (unless is.ordered in ⁠data=⁠). See also lavaan::lavCor().
absolute	logical indicating whether HTMT values should be estimated based on absolute correlations (default is TRUE). This is recommended for HTMT but required for HTMT2 (so silently ignored).
htmt2	logical indicating whether to use the geometric mean (default, appropriate for congeneric indicators) or arithmetic mean (which assumes tau-equivalence).

Value

A matrix showing HTMT(2) values (i.e., discriminant validity) between each pair of factors.

Author(s)

Ylenio Longo (University of Nottingham; [email protected])

Terrence D. Jorgensen (University of Amsterdam; [email protected])

References

Henseler, J., Ringle, C. M., & Sarstedt, M. (2015). A new criterion for assessing discriminant validity in variance-based structural equation modeling. Journal of the Academy of Marketing Science, 43(1), 115–135. doi:10.1007/s11747-014-0403-8

Roemer, E., Schuberth, F., & Henseler, J. (2021). HTMT2—An improved criterion for assessing discriminant validity in structural equation modeling. Industrial Management & Data Systems, 121(21), 2637–2650. doi:10.1108/IMDS-02-2021-0082

Voorhees, C. M., Brady, M. K., Calantone, R., & Ramirez, E. (2016). Discriminant validity testing in marketing: An analysis, causes for concern, and proposed remedies. Journal of the Academy of Marketing Science, 44(1), 119–134. doi:10.1007/s11747-015-0455-4

Examples

HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '

dat <- HolzingerSwineford1939[, paste0("x", 1:9)]
htmt(HS.model, dat)

## save covariance matrix
HS.cov <- cov(HolzingerSwineford1939[, paste0("x", 1:9)])
## HTMT using arithmetic mean
htmt(HS.model, sample.cov = HS.cov, htmt2 = FALSE)

---

Specify starting values from a lavaan output

Description

This function will save the parameter estimates of a lavaan output and impose those parameter estimates as starting values for another analysis model. The free parameters with the same names or the same labels across two models will be imposed the new starting values. This function may help to increase the chance of convergence in a complex model (e.g., multitrait-multimethod model or complex longitudinal invariance model).

Usage

imposeStart(out, expr, silent = TRUE)

Arguments

out	The lavaan output that users wish to use the parameter estimates as staring values for an analysis model
expr	The original code that users use to run a lavaan model
silent	Logical to print the parameter table with new starting values

Value

A fitted lavaan model

Author(s)

Sunthud Pornprasertmanit ([email protected])

Examples

## The following example show that the longitudinal weak invariance model
## using effect coding was not convergent with three time points but convergent
## with two time points. Thus, the parameter estimates from the model with
## two time points are used as starting values of the three time points.
## The model with new starting values is convergent properly.

weak2time <- '
	# Loadings
	f1t1 =~ LOAD1*y1t1 + LOAD2*y2t1 + LOAD3*y3t1
    f1t2 =~ LOAD1*y1t2 + LOAD2*y2t2 + LOAD3*y3t2

	# Factor Variances
	f1t1 ~~ f1t1
	f1t2 ~~ f1t2

	# Factor Covariances
	f1t1 ~~ f1t2

	# Error Variances
	y1t1 ~~ y1t1
	y2t1 ~~ y2t1
	y3t1 ~~ y3t1
	y1t2 ~~ y1t2
	y2t2 ~~ y2t2
	y3t2 ~~ y3t2

	# Error Covariances
	y1t1 ~~ y1t2
	y2t1 ~~ y2t2
	y3t1 ~~ y3t2

	# Factor Means
	f1t1 ~ NA*1
	f1t2 ~ NA*1

	# Measurement Intercepts
	y1t1 ~ INT1*1
	y2t1 ~ INT2*1
	y3t1 ~ INT3*1
	y1t2 ~ INT4*1
	y2t2 ~ INT5*1
	y3t2 ~ INT6*1

	# Constraints for Effect-coding Identification
	LOAD1 == 3 - LOAD2 - LOAD3
	INT1 == 0 - INT2 - INT3
	INT4 == 0 - INT5 - INT6
'
model2time <- lavaan(weak2time, data = exLong)

weak3time <- '
	# Loadings
	f1t1 =~ LOAD1*y1t1 + LOAD2*y2t1 + LOAD3*y3t1
    f1t2 =~ LOAD1*y1t2 + LOAD2*y2t2 + LOAD3*y3t2
    f1t3 =~ LOAD1*y1t3 + LOAD2*y2t3 + LOAD3*y3t3

	# Factor Variances
	f1t1 ~~ f1t1
	f1t2 ~~ f1t2
	f1t3 ~~ f1t3

	# Factor Covariances
	f1t1 ~~ f1t2 + f1t3
	f1t2 ~~ f1t3

	# Error Variances
	y1t1 ~~ y1t1
	y2t1 ~~ y2t1
	y3t1 ~~ y3t1
	y1t2 ~~ y1t2
	y2t2 ~~ y2t2
	y3t2 ~~ y3t2
	y1t3 ~~ y1t3
	y2t3 ~~ y2t3
	y3t3 ~~ y3t3

	# Error Covariances
	y1t1 ~~ y1t2
	y2t1 ~~ y2t2
	y3t1 ~~ y3t2
	y1t1 ~~ y1t3
	y2t1 ~~ y2t3
	y3t1 ~~ y3t3
	y1t2 ~~ y1t3
	y2t2 ~~ y2t3
	y3t2 ~~ y3t3

	# Factor Means
	f1t1 ~ NA*1
	f1t2 ~ NA*1
	f1t3 ~ NA*1

	# Measurement Intercepts
	y1t1 ~ INT1*1
	y2t1 ~ INT2*1
	y3t1 ~ INT3*1
	y1t2 ~ INT4*1
	y2t2 ~ INT5*1
	y3t2 ~ INT6*1
	y1t3 ~ INT7*1
	y2t3 ~ INT8*1
	y3t3 ~ INT9*1

	# Constraints for Effect-coding Identification
	LOAD1 == 3 - LOAD2 - LOAD3
	INT1 == 0 - INT2 - INT3
	INT4 == 0 - INT5 - INT6
	INT7 == 0 - INT8 - INT9
'
### The following command does not provide convergent result
# model3time <- lavaan(weak3time, data = exLong)

### Use starting values from the model with two time points
model3time <- imposeStart(model2time, lavaan(weak3time, data = exLong))
summary(model3time)

---

Make products of indicators using no centering, mean centering, double-mean centering, or residual centering

Description

The indProd function will make products of indicators using no centering, mean centering, double-mean centering, or residual centering. The orthogonalize function is the shortcut of the indProd function to make the residual-centered indicators products.

Usage

indProd(data, var1, var2, var3 = NULL, match = TRUE, meanC = TRUE,
  residualC = FALSE, doubleMC = TRUE, namesProd = NULL)

orthogonalize(data, var1, var2, var3 = NULL, match = TRUE,
  namesProd = NULL)

Arguments

data	The desired data to be transformed.
var1	Names or indices of the variables loaded on the first factor
var2	Names or indices of the variables loaded on the second factor
var3	Names or indices of the variables loaded on the third factor (for three-way interaction)
match	Specify TRUE to use match-paired approach (Marsh, Wen, & Hau, 2004). If FALSE, the resulting products are all possible products.
meanC	Specify TRUE for mean centering the main effect indicator before making the products
residualC	Specify TRUE for residual centering the products by the main effect indicators (Little, Bovaird, & Widaman, 2006).
doubleMC	Specify TRUE for centering the resulting products (Lin et. al., 2010)
namesProd	The names of resulting products

Value

The original data attached with the products.

Author(s)

Sunthud Pornprasertmanit ([email protected]) Alexander Schoemann (East Carolina University; [email protected])

References

Marsh, H. W., Wen, Z. & Hau, K. T. (2004). Structural equation models of latent interactions: Evaluation of alternative estimation strategies and indicator construction. Psychological Methods, 9(3), 275–300. doi:10.1037/1082-989X.9.3.275

Lin, G. C., Wen, Z., Marsh, H. W., & Lin, H. S. (2010). Structural equation models of latent interactions: Clarification of orthogonalizing and double-mean-centering strategies. Structural Equation Modeling, 17(3), 374–391. doi:10.1080/10705511.2010.488999

Little, T. D., Bovaird, J. A., & Widaman, K. F. (2006). On the merits of orthogonalizing powered and product terms: Implications for modeling interactions among latent variables. Structural Equation Modeling, 13(4), 497–519. doi:10.1207/s15328007sem1304_1

See Also

• probe2WayMC() For probing the two-way latent interaction when the results are obtained from mean-centering, or double-mean centering.

• probe3WayMC() For probing the three-way latent interaction when the results are obtained from mean-centering, or double-mean centering.

• probe2WayRC() For probing the two-way latent interaction when the results are obtained from residual-centering approach.

• probe3WayRC() For probing the two-way latent interaction when the results are obtained from residual-centering approach.

• plotProbe() Plot the simple intercepts and slopes of the latent interaction.

Examples

## Mean centering / two-way interaction / match-paired
dat <- indProd(attitude[ , -1], var1 = 1:3, var2 = 4:6)

## Residual centering / two-way interaction / match-paired
dat2 <- indProd(attitude[ , -1], var1 = 1:3, var2 = 4:6, match = FALSE,
                meanC = FALSE, residualC = TRUE, doubleMC = FALSE)

## Double-mean centering / two-way interaction / match-paired
dat3 <- indProd(attitude[ , -1], var1 = 1:3, var2 = 4:6, match = FALSE,
                meanC = TRUE, residualC = FALSE, doubleMC = TRUE)

## Mean centering / three-way interaction / match-paired
dat4 <- indProd(attitude[ , -1], var1 = 1:2, var2 = 3:4, var3 = 5:6)

## Residual centering / three-way interaction / match-paired
dat5 <- orthogonalize(attitude[ , -1], var1 = 1:2, var2 = 3:4, var3 = 5:6,
                      match = FALSE)

## Double-mean centering / three-way interaction / match-paired
dat6 <- indProd(attitude[ , -1], var1 = 1:2, var2 = 3:4, var3 = 5:6,
                match = FALSE, meanC = TRUE, residualC = TRUE,
                doubleMC = TRUE)


## To add product-indicators to multiple-imputed data sets

HSMiss <- HolzingerSwineford1939[ , c(paste0("x", 1:9), "ageyr","agemo")]
set.seed(12345)
HSMiss$x5 <- ifelse(HSMiss$x5 <= quantile(HSMiss$x5, .3), NA, HSMiss$x5)
age <- HSMiss$ageyr + HSMiss$agemo/12
HSMiss$x9 <- ifelse(age <= quantile(age, .3), NA, HSMiss$x9)
library(Amelia)
set.seed(12345)
HS.amelia <- amelia(HSMiss, m = 3, p2s = FALSE)
imps <- HS.amelia$imputations # extract a list of imputations
## apply indProd() to the list of data.frames
imps2 <- lapply(imps, indProd,
                var1 = c("x1","x2","x3"), var2 = c("x4","x5","x6"))
## verify:
lapply(imps2, head)

---

Generate data via the Kaiser-Dickman (1962) algorithm.

Description

Given a covariance matrix and sample size, generate raw data that correspond to the covariance matrix. Data can be generated to match the covariance matrix exactly, or to be a sample from the population covariance matrix.

Usage

kd(covmat, n, type = c("exact", "sample"))

Arguments

covmat	a symmetric, positive definite covariance matrix
n	the sample size for the data that will be generated
type	type of data generation. exact generates data that exactly correspond to covmat. sample treats covmat as a poulation covariance matrix, generating a sample of size n.

Details

By default, R's cov() function divides by n-1. The data generated by this algorithm result in a covariance matrix that matches covmat, but you must divide by n instead of n-1.

Value

kd returns a data matrix of dimension n by nrow(covmat).

Author(s)

Ed Merkle (University of Missouri; [email protected])

References

Kaiser, H. F. and Dickman, K. (1962). Sample and population score matrices and sample correlation matrices from an arbitrary population correlation matrix. Psychometrika, 27(2), 179–182. doi:10.1007/BF02289635

Examples

#### First Example

## Get data
dat <- HolzingerSwineford1939[ , 7:15]
hs.n <- nrow(dat)

## Covariance matrix divided by n
hscov <- ((hs.n-1)/hs.n) * cov(dat)

## Generate new, raw data corresponding to hscov
newdat <- kd(hscov, hs.n)

## Difference between new covariance matrix and hscov is minimal
newcov <- (hs.n-1)/hs.n * cov(newdat)
summary(as.numeric(hscov - newcov))

## Generate sample data, treating hscov as population matrix
newdat2 <- kd(hscov, hs.n, type = "sample")

#### Another example

## Define a covariance matrix
covmat <- matrix(0, 3, 3)
diag(covmat) <- 1.5
covmat[2:3,1] <- c(1.3, 1.7)
covmat[3,2] <- 2.1
covmat <- covmat + t(covmat)

## Generate data of size 300 that have this covariance matrix
rawdat <- kd(covmat, 300)

## Covariances are exact if we compute sample covariance matrix by
## dividing by n (vs by n - 1)
summary(as.numeric((299/300)*cov(rawdat) - covmat))

## Generate data of size 300 where covmat is the population covariance matrix
rawdat2 <- kd(covmat, 300)

---

Finding excessive kurtosis

Description

Finding excessive kurtosis ($g_{2}$) of an object

Usage

kurtosis(object, population = FALSE)

Arguments

object	A vector used to find a excessive kurtosis
population	TRUE to compute the parameter formula. FALSE to compute the sample statistic formula.

Details

The excessive kurtosis computed by default is $g_{2}$, the fourth standardized moment of the empirical distribution of object. The population parameter excessive kurtosis $\gamma_{2}$ formula is

$\gamma_{2} = \frac{\mu_{4}}{\mu^{2}_{2}} - 3,$

where $\mu_{i}$ denotes the $i$ order central moment.

The excessive kurtosis formula for sample statistic $g_{2}$ is

$g_{2} = \frac{k_{4}}{k^{2}_{2}} - 3,$

where $k_{i}$ are the $i$ order k-statistic.

The standard error of the excessive kurtosis is

$Var(\hat{g}_{2}) = \frac{24}{N}$

where $N$ is the sample size.

Value

A value of an excessive kurtosis with a test statistic if the population is specified as FALSE

Author(s)

Sunthud Pornprasertmanit ([email protected])

References

Weisstein, Eric W. (n.d.). Kurtosis. Retrieved from MathWorld–A Wolfram Web Resource: http://mathworld.wolfram.com/Kurtosis.html

See Also

• skew() Find the univariate skewness of a variable

• mardiaSkew() Find the Mardia's multivariate skewness of a set of variables

• mardiaKurtosis() Find the Mardia's multivariate kurtosis of a set of variables

Examples

kurtosis(1:5)

---

emmeans Support Functions for lavaan Models

Description

Provide emmeans support for lavaan objects

Usage

recover_data.lavaan(object, lavaan.DV, data = NULL, ...)

emm_basis.lavaan(object, trms, xlev, grid, lavaan.DV, ...)

Arguments

object	An object of class lavaan::lavaan(). See Details.
lavaan.DV	character string naming the variable(s) for which expected marginal means / trends should be produced. A vector of names indicates a multivariate outcome, treated by default as repeated measures.
data	An optional data.frame without missing values, to be passed when missing="FIML" estimation was useed, thus avoiding a reference-grid with missing values.
...	Further arguments passed to emmeans::recover_data.lm or emmeans::emm_basis.lm
trms, xlev, grid	See emmeans::emm_basis

Details

Supported DVs

lavaan.DV must be an endogenous variable, by appearing on the left-hand side of either a regression operator ("~") or an intercept operator ("~1"), or both.

lavaan.DV can also be a vector of endogenous variable, in which case they will be treated by emmeans as a multivariate outcome (often, this indicates repeated measures) represented by an additional factor named rep.meas by default. The ⁠mult.name=⁠ argument can be used to overwrite this default name.

Unsupported Models

This functionality does not support the following models:

• Multi-level models are not supported.

• Models not fit to a data.frame (i.e., models fit to a covariance matrix).

Dealing with Fixed Parameters

Fixed parameters (set with lavaan's modifiers) are treated as-is: their values are set by the users, and they have a SE of 0 (as such, they do not co-vary with any other parameter).

Dealing with Multigroup Models

If a multigroup model is supplied, a factor is added to the reference grid, the name matching the group argument supplied when fitting the model. Note that you must set nesting = NULL.

Dealing with Missing Data

Limited testing suggests that these functions do work when the model was fit to incomplete data.

Dealing with Factors

By default emmeans recognizes binary variables (0,1) as a "factor" with two levels (and not a continuous variable). With some clever contrast defenitions it should be possible to get the desired emmeans / contasts. See example below.

Author(s)

Mattan S. Ben-Shachar (Ben-Gurion University of the Negev; [email protected])

Examples

## Not run: 

  library(lavaan)
  library(emmeans)

  #### Moderation Analysis ####

  mean_sd <- function(x) mean(x) + c(-sd(x), 0, sd(x))

  model <- '
  # regressions
  Sepal.Length ~ b1 * Sepal.Width + b2 * Petal.Length + b3 * Sepal.Width:Petal.Length


  # define mean parameter label for centered math for use in simple slopes
  Sepal.Width ~ Sepal.Width.mean * 1

  # define variance parameter label for centered math for use in simple slopes
  Sepal.Width ~~ Sepal.Width.var * Sepal.Width

  # simple slopes for condition effect
  SD.below := b2 + b3 * (Sepal.Width.mean - sqrt(Sepal.Width.var))
  mean     := b2 + b3 * (Sepal.Width.mean)
  SD.above := b2 + b3 * (Sepal.Width.mean + sqrt(Sepal.Width.var))
  '

  semFit <- sem(model = model,
                data = iris)

  ## Compare simple slopes
  # From `emtrends`
  test(
    emtrends(semFit, ~ Sepal.Width, "Petal.Length",
             lavaan.DV = "Sepal.Length",
             cov.red = mean_sd)
  )

  # From lavaan
  parameterEstimates(semFit, output = "pretty")[13:15, ]
  # Identical slopes.
  # SEs differ due to lavaan estimating uncertainty of the mean / SD
  # of Sepal.Width, whereas emmeans uses the mean+-SD as is (fixed).


  #### Latent DV ####

  model <- '
  LAT1 =~ Sepal.Length + Sepal.Width

  LAT1 ~ b1 * Petal.Width + 1 * Petal.Length

  Petal.Length ~ Petal.Length.mean * 1

  V1 := 1 * Petal.Length.mean + 1 * b1
  V2 := 1 * Petal.Length.mean + 2 * b1
  '

  semFit <- sem(model = model,
                data = iris, std.lv = TRUE)

  ## Compare emmeans
  # From emmeans
  test(
    emmeans(semFit, ~ Petal.Width,
            lavaan.DV = "LAT1",
            at = list(Petal.Width = 1:2))
  )

  # From lavaan
  parameterEstimates(semFit, output = "pretty")[15:16, ]
  # Identical means.
  # SEs differ due to lavaan estimating uncertainty of the mean
  # of Petal.Length, whereas emmeans uses the mean as is.

  #### Multi-Variate DV ####

  model <- '
  ind60 =~ x1 + x2 + x3

  # metric invariance
  dem60 =~ y1 + a*y2 + b*y3 + c*y4
  dem65 =~ y5 + a*y6 + b*y7 + c*y8

  # scalar invariance
  y1 + y5 ~ d*1
  y2 + y6 ~ e*1
  y3 + y7 ~ f*1
  y4 + y8 ~ g*1

  # regressions (slopes differ: interaction with time)
  dem60 ~ b1*ind60
  dem65 ~ b2*ind60 + NA*1 + Mean.Diff*1

  # residual correlations
  y1 ~~ y5
  y2 ~~ y4 + y6
  y3 ~~ y7
  y4 ~~ y8
  y6 ~~ y8

  # conditional mean differences (besides mean(ind60) == 0)
   low := (-1*b2 + Mean.Diff) - (-1*b1) # 1 SD below M
  high := (b2 + Mean.Diff) - b1         # 1 SD above M
'

  semFit <- sem(model, data = PoliticalDemocracy)


  ## Compare contrasts
  # From emmeans
  emmeans(semFit, pairwise ~ rep.meas|ind60,
          lavaan.DV = c("dem60","dem65"),
          at = list(ind60 = c(-1,1)))[[2]]

  # From lavaan
  parameterEstimates(semFit, output = "pretty")[49:50, ]


  #### Multi Group ####

  model <- 'x1 ~ c(int1, int2)*1 + c(b1, b2)*ageyr
  diff_11 := (int2 + b2*11) - (int1 + b1*11)
  diff_13 := (int2 + b2*13) - (int1 + b1*13)
  diff_15 := (int2 + b2*15) - (int1 + b1*15)
'
  semFit <- sem(model, group = "school", data = HolzingerSwineford1939)


  ## Compare contrasts
  # From emmeans (note `nesting = NULL`)
  emmeans(semFit, pairwise ~ school | ageyr, lavaan.DV = "x1",
          at = list(ageyr = c(11, 13, 15)), nesting = NULL)[[2]]

  # From lavaan
  parameterEstimates(semFit, output = "pretty")

  #### Dealing with factors ####

  warpbreaks <- cbind(warpbreaks,
                      model.matrix(~ wool + tension, data = warpbreaks))

  model <- "
  # Split for convenience
  breaks ~ 1
  breaks ~ woolB
  breaks ~ tensionM + tensionH
  breaks ~ woolB:tensionM + woolB:tensionH
  "

  semFit <- sem(model, warpbreaks)

  ## Compare contrasts
  # From lm -> emmeans
  lmFit <- lm(breaks ~ wool * tension, data = warpbreaks)
  lmEM <- emmeans(lmFit, ~ tension + wool)
  contrast(lmEM, method = data.frame(L_all = c(-1, .05, 0.5),
                                     M_H   = c(0, 1, -1)), by = "wool")

  # From lavaan -> emmeans
  lavEM <- emmeans(semFit, ~ tensionM + tensionH + woolB,
                   lavaan.DV = "breaks")
  contrast(lavEM,
           method = list(
             "L_all|A" = c(c(-1, .05, 0.5, 0), rep(0, 4)),
             "M_H  |A" = c(c(0, 1, -1, 0),     rep(0, 4)),
             "L_all|A" = c(rep(0, 4),          c(-1, .05, 0.5, 0)),
             "M_H  |A" = c(rep(0, 4),          c(0, 1, -1, 0))
           ))

## End(Not run)

---

Find standardized factor loading from coefficient alpha

Description

Find standardized factor loading from coefficient alpha assuming that all items have equal loadings.

Usage

loadingFromAlpha(alpha, ni)

Arguments

alpha	A desired coefficient alpha value.
ni	A desired number of items.

Value

result

The standardized factor loadings that make desired coefficient alpha with specified number of items.

Author(s)

Sunthud Pornprasertmanit ([email protected])

Examples

loadingFromAlpha(0.8, 4)

---

Calculate Population Moments for Ordinal Data Treated as Numeric

Description

This function calculates ordinal-scale moments implied by LRV-scale moments

Usage

lrv2ord(Sigma, Mu, thresholds, cWts)

Arguments

Sigma	Population covariance matrix(), with variable names saved in the dimnames() attribute.
Mu	Optional numeric vector of population means. If missing, all means will be set to zero.
thresholds	Either a single numeric vector of population thresholds used to discretize each normally distributed variable, or a named list of each discretized variable's vector of thresholds. The discretized variables may be a subset of all variables in Sigma if the remaining variables are intended to be observed rather than latent normally distributed variables.
cWts	Optional (default when missing is to use 0 for the lowest category, followed by successive integers for each higher category). Either a single numeric vector of category weights (if they are identical across all variables) or a named list of each discretized variable's vector of category weights.

Details

Binary and ordinal data are frequently accommodated in SEM by incorporating a threshold model that links each observed categorical response variable to a corresponding latent response variable that is typically assumed to be normally distributed (Kamata & Bauer, 2008; Wirth & Edwards, 2007). This function can be useful for real-data analysis or for designing Monte Carlo simulations, as described by Jorgensen and Johnson (2022).

Value

A list including the LRV-scale population moments (means, covariance matrix, correlation matrix, and thresholds), the category weights, a data.frame of implied univariate moments (means, SDs, skewness, and excess kurtosis (i.e., in excess of 3, which is the kurtosis of the normal distribution) for discretized data treated as numeric, and the implied covariance and correlation matrix of discretized data treated as numeric.

Author(s)

Terrence D. Jorgensen (University of Amsterdam; [email protected])

Andrew R. Johnson (Curtin University; [email protected])

References

Jorgensen, T. D., & Johnson, A. R. (2022). How to derive expected values of structural equation model parameters when treating discrete data as continuous. Structural Equation Modeling, 29(4), 639–650. doi:10.1080/10705511.2021.1988609

Kamata, A., & Bauer, D. J. (2008). A note on the relation between factor analytic and item response theory models. Structural Equation Modeling, 15(1), 136–153. doi:10.1080/10705510701758406

Wirth, R. J., & Edwards, M. C. (2007). Item factor analysis: Current approaches and future directions. Psychological Methods, 12(1), 58–79. doi:10.1037/1082-989X.12.1.58

Examples

## SCENARIO 1: DIRECTLY SPECIFY POPULATION PARAMETERS

## specify population model in LISREL matrices
Nu <- rep(0, 4)
Alpha <- c(1, -0.5)
Lambda <- matrix(c(1, 1, 0, 0, 0, 0, 1, 1), nrow = 4, ncol = 2,
                 dimnames = list(paste0("y", 1:4), paste0("eta", 1:2)))
Psi <- diag(c(1, .75))
Theta <- diag(4)
Beta <- matrix(c(0, .5, 0, 0), nrow = 2, ncol = 2)

## calculate model-implied population means and covariance matrix
## of latent response variables (LRVs)
IB <- solve(diag(2) - Beta) # to save time and space
Mu_LRV <- Nu + Lambda %*% IB %*% Alpha
Sigma_LRV <- Lambda %*% IB %*% Psi %*% t(IB) %*% t(Lambda) + Theta

## Specify (unstandardized) thresholds to discretize normally distributed data
## generated from Mu_LRV and Sigma_LRV, based on marginal probabilities
PiList <- list(y1 = c(.25, .5, .25),
               y2 = c(.17, .33, .33, .17),
               y3 = c(.1, .2, .4, .2, .1),
               ## make final variable highly asymmetric
               y4 = c(.33, .25, .17, .12, .08, .05))
sapply(PiList, sum) # all sum to 100%
CumProbs <- sapply(PiList, cumsum)
## unstandardized thresholds
TauList <- mapply(qnorm, p = lapply(CumProbs, function(x) x[-length(x)]),
                  m = Mu_LRV, sd = sqrt(diag(Sigma_LRV)))
for (i in 1:4) names(TauList[[i]]) <- paste0(names(TauList)[i], "|t",
                                             1:length(TauList[[i]]))

## assign numeric weights to each category (optional, see default)
NumCodes <- list(y1 = c(-0.5, 0, 0.5), y2 = 0:3, y3 = 1:5, y4 = 1:6)


## Calculate Population Moments for Numerically Coded Ordinal Variables
lrv2ord(Sigma = Sigma_LRV, Mu = Mu_LRV, thresholds = TauList, cWts = NumCodes)


## SCENARIO 2: USE ESTIMATED PARAMETERS AS POPULATION

data(datCat) # already stored as c("ordered","factor")
fit <- cfa(' f =~ 1*u1 + 1*u2 + 1*u3 + 1*u4 ', data = datCat)
lrv2ord(Sigma = fit, thresholds = fit) # use same fit for both
## or use estimated thresholds with specified parameters, but note that
## lrv2ord() will only extract standardized thresholds
dimnames(Sigma_LRV) <- list(paste0("u", 1:4), paste0("u", 1:4))
lrv2ord(Sigma = cov2cor(Sigma_LRV), thresholds = fit)

---

Finding Mardia's multivariate kurtosis

Description

Finding Mardia's multivariate kurtosis of multiple variables

Usage

mardiaKurtosis(dat, use = "everything")

Arguments

dat	The target matrix or data frame with multiple variables
use	Missing data handling method from the stats::cov() function.

Details

The Mardia's multivariate kurtosis formula (Mardia, 1970) is

$b_{2, d} = \frac{1}{n}\sum^n_{i=1}\left[ \left(\bold{X}_i - \bold{\bar{X}} \right)^{'} \bold{S}^{-1} \left(\bold{X}_i - \bold{\bar{X}} \right) \right]^2,$

where $d$ is the number of variables, $X$ is the target dataset with multiple variables, $n$ is the sample size, $\bold{S}$ is the sample covariance matrix of the target dataset, and $\bold{\bar{X}}$ is the mean vectors of the target dataset binded in $n$ rows. When the population multivariate kurtosis is normal, the $b_{2,d}$ is asymptotically distributed as normal distribution with the mean of $d(d + 2)$ and variance of $8d(d + 2)/n$.

Value

A value of a Mardia's multivariate kurtosis with a test statistic

Author(s)

Sunthud Pornprasertmanit ([email protected])

References

Mardia, K. V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika, 57(3), 519–530. doi:10.2307/2334770

See Also

• skew() Find the univariate skewness of a variable

• kurtosis() Find the univariate excessive kurtosis of a variable

• mardiaSkew() Find the Mardia's multivariate skewness of a set of variables

Examples

library(lavaan)
mardiaKurtosis(HolzingerSwineford1939[ , paste0("x", 1:9)])

---

Finding Mardia's multivariate skewness

Description

Finding Mardia's multivariate skewness of multiple variables

Usage

mardiaSkew(dat, use = "everything")

Arguments

dat	The target matrix or data frame with multiple variables
use	Missing data handling method from the stats::cov() function.

Details

The Mardia's multivariate skewness formula (Mardia, 1970) is

$b_{1, d} = \frac{1}{n^2}\sum^n_{i=1}\sum^n_{j=1}\left[ \left(\bold{X}_i - \bold{\bar{X}} \right)^{'} \bold{S}^{-1} \left(\bold{X}_j - \bold{\bar{X}} \right) \right]^3,$

where $d$ is the number of variables, $X$ is the target dataset with multiple variables, $n$ is the sample size, $\bold{S}$ is the sample covariance matrix of the target dataset, and $\bold{\bar{X}}$ is the mean vectors of the target dataset binded in $n$ rows. When the population multivariate skewness is normal, the $\frac{n}{6}b_{1,d}$ is asymptotically distributed as $\chi^2$ distribution with $d(d + 1)(d + 2)/6$ degrees of freedom.

Value

A value of a Mardia's multivariate skewness with a test statistic

Author(s)

Sunthud Pornprasertmanit ([email protected])

References

Mardia, K. V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika, 57(3), 519–530. doi:10.2307/2334770

See Also

• skew() Find the univariate skewness of a variable

• kurtosis() Find the univariate excessive kurtosis of a variable

• mardiaKurtosis() Find the Mardia's multivariate kurtosis of a set of variables

Examples

library(lavaan)
mardiaSkew(HolzingerSwineford1939[ , paste0("x", 1:9)])

---