Nothing
R/lav_test_yuan_bentler.R
In lavaan: Latent Variable Analysis
Defines functions
lav_test_yuan_bentler_mplus_trace
lav_test_yuan_bentler_trace
lav_test_yuan_bentler
# - 0.6-13: fix multiple-group UG^2 bug (reported by Gronneberg, Foldnes and
# Moss) when Satterthwaite = TRUE, ngroups > 1, and eq constraints.
#
# Note however that Satterthwaite = FALSE always (for now), so
# the fix has no (visible) effect
lav_test_yuan_bentler <- function(lavobject = NULL,
lavsamplestats = NULL,
lavmodel = NULL,
lavimplied = NULL,
lavh1 = NULL,
lavoptions = NULL,
lavdata = NULL,
TEST.unscaled = NULL,
E.inv = NULL,
B0.group = NULL,
test = "yuan.bentler",
mimic = "lavaan",
#method = "default",
ug2.old.approach = FALSE,
return.ugamma = FALSE) {
TEST <- list()
if(!is.null(lavobject)) {
lavsamplestats <- lavobject@SampleStats
lavmodel <- lavobject@Model
lavoptions <- lavobject@Options
lavpartable <- lavobject@ParTable
lavimplied <- lavobject@implied
lavh1 <- lavobject@h1
lavdata <- lavobject@Data
TEST$standard <- lavobject@test[[1]]
} else {
TEST$standard <- TEST.unscaled
}
# ug2.old.approach
if(missing(ug2.old.approach)) {
if(!is.null(lavoptions$ug2.old.approach)) {
ug2.old.approach <- lavoptions$ug2.old.approach
} else {
ug2.old.approach <- FALSE
}
}
# E.inv ok?
if( length(lavoptions$information) == 1L &&
length(lavoptions$h1.information) == 1L &&
length(lavoptions$observed.information) == 1L) {
E.inv.recompute <- FALSE
} else if( (lavoptions$information[1] == lavoptions$information[2]) &&
(lavoptions$h1.information[1] == lavoptions$h1.information[2]) &&
(lavoptions$information[2] == "expected" ||
lavoptions$observed.information[1] ==
lavoptions$observed.information[2]) ) {
E.inv.recompute <- FALSE
} else {
E.inv.recompute <- TRUE
# change information options
lavoptions$information[1] <- lavoptions$information[2]
lavoptions$h1.information[1] <- lavoptions$h1.information[2]
lavoptions$observed.information[1] <- lavoptions$observed.information[2]
}
if(!is.null(E.inv)) {
E.inv.recompute <- FALSE # user-provided
}
# check test
if(!all(test %in% c("yuan.bentler",
"yuan.bentler.mplus"))) {
warning("lavaan WARNING: test must be one of `yuan.bentler', or `yuan.bentler.mplus'; will use `yuan.bentler' only")
test <- "yuan.bentler"
}
# information
information <- lavoptions$information[1]
# ndat
ndat <- numeric(lavsamplestats@ngroups)
# do we have E.inv?
if(is.null(E.inv) || E.inv.recompute) {
E.inv <- try(lav_model_information(lavmodel = lavmodel,
lavsamplestats = lavsamplestats,
lavdata = lavdata,
lavimplied = lavimplied,
lavoptions = lavoptions,
extra = FALSE,
augmented = TRUE,
inverted = TRUE),
silent = TRUE)
if(inherits(E.inv, "try-error")) {
if(return.ugamma) {
warning("lavaan WARNING: could not invert information matrix needed for UGamma\n")
return(NULL)
} else {
TEST$standard$stat <- as.numeric(NA)
TEST$standard$stat.group <- rep(as.numeric(NA), lavdata@ngroups)
TEST$standard$pvalue <- as.numeric(NA)
TEST[[test[1]]] <- c(TEST$standard,
scaling.factor = as.numeric(NA),
shift.parameter = as.numeric(NA),
label = character(0))
warning("lavaan WARNING: could not invert information [matrix needed for robust test statistic\n")
TEST[[test[1]]]$test <- test[1] # to prevent lavTestLRT error when robust test is detected for some but not all models
return(TEST)
}
}
}
# catch df == 0
if(TEST$standard$df == 0L || TEST$standard$df < 0) {
TEST[[test[1]]] <- c(TEST$standard, scaling.factor = as.numeric(NA),
label = character(0))
TEST[[test[1]]]$test <- test[1] # to prevent lavTestLRT error when robust test is detected for some but not all models
return(TEST)
}
# mean and variance adjusted?
Satterthwaite <- FALSE # for now
#if(any(test %in% c("mean.var.adjusted", "scaled.shifted"))) {
# Satterthwaite <- TRUE
#}
# FIXME: should we not always use 'unstructured' here?
# if the model is, say, the independence model, the
# 'structured' information (A1) will be so far away from B1
# that we will end up with 'NA'
h1.options <- lavoptions
if(test == "yuan.bentler.mplus") {
# always 'unstructured' H1 information
h1.options$h1.information <- "unstructured"
}
# A1 is usually expected or observed
A1.group <- lav_model_h1_information(lavmodel = lavmodel,
lavsamplestats = lavsamplestats,
lavdata = lavdata,
lavimplied = lavimplied,
lavh1 = lavh1,
lavoptions = h1.options)
# B1 is always first.order
B1.group <- lav_model_h1_information_firstorder(lavmodel = lavmodel,
lavsamplestats = lavsamplestats,
lavdata = lavdata,
lavimplied = lavimplied,
lavh1 = lavh1,
lavoptions = h1.options)
if(test == "yuan.bentler.mplus") {
if(is.null(B0.group)) {
B0 <- lav_model_information_firstorder(lavmodel = lavmodel,
lavsamplestats = lavsamplestats,
lavdata = lavdata,
lavh1 = lavh1,
lavoptions = lavoptions,
extra = TRUE,
check.pd = FALSE,
augmented = FALSE,
inverted = FALSE)
B0.group <- attr(B0, "B0.group")
}
trace.UGamma <-
lav_test_yuan_bentler_mplus_trace(lavsamplestats = lavsamplestats,
A1.group = A1.group,
B1.group = B1.group,
B0.group = B0.group,
E.inv = E.inv,
meanstructure = lavmodel@meanstructure)
} else if(test == "yuan.bentler") {
# compute Delta
Delta <- computeDelta(lavmodel = lavmodel)
# compute Omega/Gamma
Omega <- lav_model_h1_omega(lavmodel = lavmodel,
lavsamplestats = lavsamplestats,
lavdata = lavdata,
lavimplied = lavimplied,
lavh1 = lavh1,
lavoptions = lavoptions)
# compute trace 'U %*% Gamma' (or 'U %*% Omega')
trace.UGamma <- lav_test_yuan_bentler_trace(
lavsamplestats = lavsamplestats,
meanstructure = lavmodel@meanstructure,
A1.group = A1.group,
B1.group = B1.group,
Delta = Delta,
Omega = Omega,
E.inv = E.inv,
ug2.old.approach = ug2.old.approach,
Satterthwaite = FALSE) # for now
}
# unscaled test
df <- TEST$standard$df
scaling.factor <- trace.UGamma / df
if(scaling.factor < 0) scaling.factor <- as.numeric(NA)
chisq.scaled <- TEST$standard$stat / scaling.factor
pvalue.scaled <- 1 - pchisq(chisq.scaled, df)
ndat <- sum(attr(trace.UGamma, "h1.ndat"))
npar <- lavmodel@nx.free
scaling.factor.h1 <- sum( attr(trace.UGamma, "h1") ) / ndat
scaling.factor.h0 <- sum( attr(trace.UGamma, "h0") ) / npar
trace.UGamma2 <- attr(trace.UGamma, "trace.UGamma2")
attributes(trace.UGamma) <- NULL
if("yuan.bentler" %in% test) {
TEST$yuan.bentler <-
list(test = test,
stat = chisq.scaled,
stat.group = (TEST$standard$stat.group /
scaling.factor),
df = df,
pvalue = pvalue.scaled,
scaling.factor = scaling.factor,
scaling.factor.h1 = scaling.factor.h1,
scaling.factor.h0 = scaling.factor.h0,
label = "Yuan-Bentler correction",
trace.UGamma = trace.UGamma,
trace.UGamma2 = trace.UGamma2)
} else if("yuan.bentler.mplus" %in% test) {
TEST$yuan.bentler.mplus <-
list(test = test,
stat = chisq.scaled,
stat.group = (TEST$standard$stat.group /
scaling.factor),
df = df,
pvalue = pvalue.scaled,
scaling.factor = scaling.factor,
scaling.factor.h1 = scaling.factor.h1,
scaling.factor.h0 = scaling.factor.h0,
label =
"Yuan-Bentler correction (Mplus variant)",
trace.UGamma = trace.UGamma,
trace.UGamma2 = as.numeric(NA))
}
TEST
}
lav_test_yuan_bentler_trace <- function(lavsamplestats = lavsamplestats,
meanstructure = TRUE,
A1.group = NULL,
B1.group = NULL,
Delta = NULL,
Omega = NULL,
E.inv = NULL,
ug2.old.approach = FALSE,
Satterthwaite = FALSE) {
# we always assume a meanstructure (nope, not any longer, since 0.6)
#meanstructure <- TRUE
ngroups <- lavsamplestats@ngroups
trace.h1 <- attr(Omega, "trace.h1")
h1.ndat <- attr(Omega, "h1.ndat")
if(ug2.old.approach || !Satterthwaite) {
trace.UGamma <- numeric( ngroups )
trace.UGamma2 <- numeric( ngroups )
trace.h0 <- numeric( ngroups )
for(g in 1:ngroups) {
fg <- lavsamplestats@nobs[[g]]/lavsamplestats@ntotal
A1 <- A1.group[[g]] * fg
B1 <- B1.group[[g]] * fg
DELTA <- Delta[[g]]
Gamma.g <- Omega[[g]] / fg
D.Einv.tD <- DELTA %*% tcrossprod(E.inv, DELTA)
#trace.h1[g] <- sum( B1 * t( A1.inv ) )
# fg cancels out: trace.h1[g] <- sum( fg*B1 * t( 1/fg*A1.inv ) )
trace.h0[g] <- sum( B1 * D.Einv.tD )
#trace.UGamma[g] <- trace.h1[g] - trace.h0[g]
U <- A1 - A1 %*% D.Einv.tD %*% A1
trace.UGamma[g] <- sum(U * Gamma.g)
if(Satterthwaite) {
UG <- U %*% Gamma.g
trace.UGamma2[g] <- sum(UG * t(UG))
}
} # g
trace.UGamma <- sum(trace.UGamma)
attr(trace.UGamma, "h1") <- trace.h1
attr(trace.UGamma, "h0") <- trace.h0
attr(trace.UGamma, "h1.ndat") <- h1.ndat
if(Satterthwaite) {
attr(trace.UGamma, "trace.UGamma2") <- sum(trace.UGamma2)
}
} else {
trace.UGamma <- trace.UGamma2 <- UG <- as.numeric(NA)
fg <- unlist(lavsamplestats@nobs)/lavsamplestats@ntotal
#if(Satterthwaite) {
A1.f <- A1.group
for(g in 1:ngroups) {
A1.f[[g]] <- A1.group[[g]] * fg[g]
}
A1.all <- lav_matrix_bdiag(A1.f)
B1.f <- B1.group
for(g in 1:ngroups) {
B1.f[[g]] <- B1.group[[g]] * fg[g]
}
B1.all <- lav_matrix_bdiag(B1.f)
Gamma.f <- Omega
for(g in 1:ngroups) {
Gamma.f[[g]] <- 1/fg[g] * Omega[[g]]
}
Gamma.all <- lav_matrix_bdiag(Gamma.f)
Delta.all <- do.call("rbind", Delta)
D.Einv.tD <- Delta.all %*% tcrossprod(E.inv, Delta.all)
trace.h0 <- sum( B1.all * D.Einv.tD )
U.all <- A1.all - A1.all %*% D.Einv.tD %*% A1.all
trace.UGamma <- sum(U.all * Gamma.all)
attr(trace.UGamma, "h1") <- sum(trace.h1)
attr(trace.UGamma, "h0") <- trace.h0
attr(trace.UGamma, "h1.ndat") <- sum(h1.ndat)
if(Satterthwaite) {
UG <- U.all %*% Gamma.all
trace.UGamma2 <- sum(UG * t(UG))
attr(trace.UGamma, "trace.UGamma2") <- trace.UGamma2
}
# } else {
}
trace.UGamma
}
lav_test_yuan_bentler_mplus_trace <- function(lavsamplestats=NULL,
A1.group = NULL,
B1.group = NULL,
B0.group=NULL,
E.inv=NULL,
meanstructure = TRUE) {
# typical for Mplus:
# - do NOT use the YB formula, but use an approximation
# relying on A0 ~= Delta' A1 Delta and the same for B0
#
# NOTE: if A0 is based on the hessian, then A0 only approximates
# Delta' A1 Delta
#
# - always use h1.information = "unstructured"!!!
ngroups <- lavsamplestats@ngroups
trace.UGamma <- numeric( lavsamplestats@ngroups )
trace.h1 <- numeric( lavsamplestats@ngroups )
trace.h0 <- numeric( lavsamplestats@ngroups )
h1.ndat <- numeric( lavsamplestats@ngroups )
for(g in 1:lavsamplestats@ngroups) {
# group weight
fg <- lavsamplestats@nobs[[g]]/lavsamplestats@ntotal
A1 <- A1.group[[g]]
B1 <- B1.group[[g]]
# mask independent 'fixed-x' variables
zero.idx <- which(diag(A1) == 0)
if(length(zero.idx) > 0L) {
A1.inv <- matrix(0, nrow(A1), ncol(A1))
a1 <- A1[-zero.idx, -zero.idx]
a1.inv <- solve(a1)
A1.inv[-zero.idx, -zero.idx] <- a1.inv
} else {
A1.inv <- solve(A1)
}
h1.ndat[g] <- ncol(A1) - length(zero.idx)
# if data is complete, why not just A1 %*% Gamma?
trace.h1[g] <- sum( B1 * t( A1.inv ) )
trace.h0[g] <- fg * sum( B0.group[[g]] * t(E.inv) )
trace.UGamma[g] <- (trace.h1[g] - trace.h0[g])
}
# we take the sum here
trace.UGamma <- sum(trace.UGamma)
attr(trace.UGamma, "h1") <- trace.h1
attr(trace.UGamma, "h0") <- trace.h0
attr(trace.UGamma, "h1.ndat") <- h1.ndat
trace.UGamma
}
Try the lavaan package in your browser
Any scripts or data that you put into this service are public.
lavaan documentation built on July 26, 2023, 5:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions lav_test_yuan_bentler_mplus_trace lav_test_yuan_bentler_trace lav_test_yuan_bentler
# - 0.6-13: fix multiple-group UG^2 bug (reported by Gronneberg, Foldnes and
# Moss) when Satterthwaite = TRUE, ngroups > 1, and eq constraints.
#
# Note however that Satterthwaite = FALSE always (for now), so
# the fix has no (visible) effect
lav_test_yuan_bentler <- function(lavobject = NULL,
lavsamplestats = NULL,
lavmodel = NULL,
lavimplied = NULL,
lavh1 = NULL,
lavoptions = NULL,
lavdata = NULL,
TEST.unscaled = NULL,
E.inv = NULL,
B0.group = NULL,
test = "yuan.bentler",
mimic = "lavaan",
#method = "default",
ug2.old.approach = FALSE,
return.ugamma = FALSE) {
TEST <- list()
if(!is.null(lavobject)) {
lavsamplestats <- lavobject@SampleStats
lavmodel <- lavobject@Model
lavoptions <- lavobject@Options
lavpartable <- lavobject@ParTable
lavimplied <- lavobject@implied
lavh1 <- lavobject@h1
lavdata <- lavobject@Data
TEST$standard <- lavobject@test[[1]]
} else {
TEST$standard <- TEST.unscaled
}
# ug2.old.approach
if(missing(ug2.old.approach)) {
if(!is.null(lavoptions$ug2.old.approach)) {
ug2.old.approach <- lavoptions$ug2.old.approach
} else {
ug2.old.approach <- FALSE
}
}
# E.inv ok?
if( length(lavoptions$information) == 1L &&
length(lavoptions$h1.information) == 1L &&
length(lavoptions$observed.information) == 1L) {
E.inv.recompute <- FALSE
} else if( (lavoptions$information[1] == lavoptions$information[2]) &&
(lavoptions$h1.information[1] == lavoptions$h1.information[2]) &&
(lavoptions$information[2] == "expected" ||
lavoptions$observed.information[1] ==
lavoptions$observed.information[2]) ) {
E.inv.recompute <- FALSE
} else {
E.inv.recompute <- TRUE
# change information options
lavoptions$information[1] <- lavoptions$information[2]
lavoptions$h1.information[1] <- lavoptions$h1.information[2]
lavoptions$observed.information[1] <- lavoptions$observed.information[2]
}
if(!is.null(E.inv)) {
E.inv.recompute <- FALSE # user-provided
}
# check test
if(!all(test %in% c("yuan.bentler",
"yuan.bentler.mplus"))) {
warning("lavaan WARNING: test must be one of `yuan.bentler', or `yuan.bentler.mplus'; will use `yuan.bentler' only")
test <- "yuan.bentler"
}
# information
information <- lavoptions$information[1]
# ndat
ndat <- numeric(lavsamplestats@ngroups)
# do we have E.inv?
if(is.null(E.inv) || E.inv.recompute) {
E.inv <- try(lav_model_information(lavmodel = lavmodel,
lavsamplestats = lavsamplestats,
lavdata = lavdata,
lavimplied = lavimplied,
lavoptions = lavoptions,
extra = FALSE,
augmented = TRUE,
inverted = TRUE),
silent = TRUE)
if(inherits(E.inv, "try-error")) {
if(return.ugamma) {
warning("lavaan WARNING: could not invert information matrix needed for UGamma\n")
return(NULL)
} else {
TEST$standard$stat <- as.numeric(NA)
TEST$standard$stat.group <- rep(as.numeric(NA), lavdata@ngroups)
TEST$standard$pvalue <- as.numeric(NA)
TEST[[test[1]]] <- c(TEST$standard,
scaling.factor = as.numeric(NA),
shift.parameter = as.numeric(NA),
label = character(0))
warning("lavaan WARNING: could not invert information [matrix needed for robust test statistic\n")
TEST[[test[1]]]$test <- test[1] # to prevent lavTestLRT error when robust test is detected for some but not all models
return(TEST)
}
}
}
# catch df == 0
if(TEST$standard$df == 0L || TEST$standard$df < 0) {
TEST[[test[1]]] <- c(TEST$standard, scaling.factor = as.numeric(NA),
label = character(0))
TEST[[test[1]]]$test <- test[1] # to prevent lavTestLRT error when robust test is detected for some but not all models
return(TEST)
}
# mean and variance adjusted?
Satterthwaite <- FALSE # for now
#if(any(test %in% c("mean.var.adjusted", "scaled.shifted"))) {
# Satterthwaite <- TRUE
#}
# FIXME: should we not always use 'unstructured' here?
# if the model is, say, the independence model, the
# 'structured' information (A1) will be so far away from B1
# that we will end up with 'NA'
h1.options <- lavoptions
if(test == "yuan.bentler.mplus") {
# always 'unstructured' H1 information
h1.options$h1.information <- "unstructured"
}
# A1 is usually expected or observed
A1.group <- lav_model_h1_information(lavmodel = lavmodel,
lavsamplestats = lavsamplestats,
lavdata = lavdata,
lavimplied = lavimplied,
lavh1 = lavh1,
lavoptions = h1.options)
# B1 is always first.order
B1.group <- lav_model_h1_information_firstorder(lavmodel = lavmodel,
lavsamplestats = lavsamplestats,
lavdata = lavdata,
lavimplied = lavimplied,
lavh1 = lavh1,
lavoptions = h1.options)
if(test == "yuan.bentler.mplus") {
if(is.null(B0.group)) {
B0 <- lav_model_information_firstorder(lavmodel = lavmodel,
lavsamplestats = lavsamplestats,
lavdata = lavdata,
lavh1 = lavh1,
lavoptions = lavoptions,
extra = TRUE,
check.pd = FALSE,
augmented = FALSE,
inverted = FALSE)
B0.group <- attr(B0, "B0.group")
}
trace.UGamma <-
lav_test_yuan_bentler_mplus_trace(lavsamplestats = lavsamplestats,
A1.group = A1.group,
B1.group = B1.group,
B0.group = B0.group,
E.inv = E.inv,
meanstructure = lavmodel@meanstructure)
} else if(test == "yuan.bentler") {
# compute Delta
Delta <- computeDelta(lavmodel = lavmodel)
# compute Omega/Gamma
Omega <- lav_model_h1_omega(lavmodel = lavmodel,
lavsamplestats = lavsamplestats,
lavdata = lavdata,
lavimplied = lavimplied,
lavh1 = lavh1,
lavoptions = lavoptions)
# compute trace 'U %*% Gamma' (or 'U %*% Omega')
trace.UGamma <- lav_test_yuan_bentler_trace(
lavsamplestats = lavsamplestats,
meanstructure = lavmodel@meanstructure,
A1.group = A1.group,
B1.group = B1.group,
Delta = Delta,
Omega = Omega,
E.inv = E.inv,
ug2.old.approach = ug2.old.approach,
Satterthwaite = FALSE) # for now
}
# unscaled test
df <- TEST$standard$df
scaling.factor <- trace.UGamma / df
if(scaling.factor < 0) scaling.factor <- as.numeric(NA)
chisq.scaled <- TEST$standard$stat / scaling.factor
pvalue.scaled <- 1 - pchisq(chisq.scaled, df)
ndat <- sum(attr(trace.UGamma, "h1.ndat"))
npar <- lavmodel@nx.free
scaling.factor.h1 <- sum( attr(trace.UGamma, "h1") ) / ndat
scaling.factor.h0 <- sum( attr(trace.UGamma, "h0") ) / npar
trace.UGamma2 <- attr(trace.UGamma, "trace.UGamma2")
attributes(trace.UGamma) <- NULL
if("yuan.bentler" %in% test) {
TEST$yuan.bentler <-
list(test = test,
stat = chisq.scaled,
stat.group = (TEST$standard$stat.group /
scaling.factor),
df = df,
pvalue = pvalue.scaled,
scaling.factor = scaling.factor,
scaling.factor.h1 = scaling.factor.h1,
scaling.factor.h0 = scaling.factor.h0,
label = "Yuan-Bentler correction",
trace.UGamma = trace.UGamma,
trace.UGamma2 = trace.UGamma2)
} else if("yuan.bentler.mplus" %in% test) {
TEST$yuan.bentler.mplus <-
list(test = test,
stat = chisq.scaled,
stat.group = (TEST$standard$stat.group /
scaling.factor),
df = df,
pvalue = pvalue.scaled,
scaling.factor = scaling.factor,
scaling.factor.h1 = scaling.factor.h1,
scaling.factor.h0 = scaling.factor.h0,
label =
"Yuan-Bentler correction (Mplus variant)",
trace.UGamma = trace.UGamma,
trace.UGamma2 = as.numeric(NA))
}
TEST
}
lav_test_yuan_bentler_trace <- function(lavsamplestats = lavsamplestats,
meanstructure = TRUE,
A1.group = NULL,
B1.group = NULL,
Delta = NULL,
Omega = NULL,
E.inv = NULL,
ug2.old.approach = FALSE,
Satterthwaite = FALSE) {
# we always assume a meanstructure (nope, not any longer, since 0.6)
#meanstructure <- TRUE
ngroups <- lavsamplestats@ngroups
trace.h1 <- attr(Omega, "trace.h1")
h1.ndat <- attr(Omega, "h1.ndat")
if(ug2.old.approach || !Satterthwaite) {
trace.UGamma <- numeric( ngroups )
trace.UGamma2 <- numeric( ngroups )
trace.h0 <- numeric( ngroups )
for(g in 1:ngroups) {
fg <- lavsamplestats@nobs[[g]]/lavsamplestats@ntotal
A1 <- A1.group[[g]] * fg
B1 <- B1.group[[g]] * fg
DELTA <- Delta[[g]]
Gamma.g <- Omega[[g]] / fg
D.Einv.tD <- DELTA %*% tcrossprod(E.inv, DELTA)
#trace.h1[g] <- sum( B1 * t( A1.inv ) )
# fg cancels out: trace.h1[g] <- sum( fg*B1 * t( 1/fg*A1.inv ) )
trace.h0[g] <- sum( B1 * D.Einv.tD )
#trace.UGamma[g] <- trace.h1[g] - trace.h0[g]
U <- A1 - A1 %*% D.Einv.tD %*% A1
trace.UGamma[g] <- sum(U * Gamma.g)
if(Satterthwaite) {
UG <- U %*% Gamma.g
trace.UGamma2[g] <- sum(UG * t(UG))
}
} # g
trace.UGamma <- sum(trace.UGamma)
attr(trace.UGamma, "h1") <- trace.h1
attr(trace.UGamma, "h0") <- trace.h0
attr(trace.UGamma, "h1.ndat") <- h1.ndat
if(Satterthwaite) {
attr(trace.UGamma, "trace.UGamma2") <- sum(trace.UGamma2)
}
} else {
trace.UGamma <- trace.UGamma2 <- UG <- as.numeric(NA)
fg <- unlist(lavsamplestats@nobs)/lavsamplestats@ntotal
#if(Satterthwaite) {
A1.f <- A1.group
for(g in 1:ngroups) {
A1.f[[g]] <- A1.group[[g]] * fg[g]
}
A1.all <- lav_matrix_bdiag(A1.f)
B1.f <- B1.group
for(g in 1:ngroups) {
B1.f[[g]] <- B1.group[[g]] * fg[g]
}
B1.all <- lav_matrix_bdiag(B1.f)
Gamma.f <- Omega
for(g in 1:ngroups) {
Gamma.f[[g]] <- 1/fg[g] * Omega[[g]]
}
Gamma.all <- lav_matrix_bdiag(Gamma.f)
Delta.all <- do.call("rbind", Delta)
D.Einv.tD <- Delta.all %*% tcrossprod(E.inv, Delta.all)
trace.h0 <- sum( B1.all * D.Einv.tD )
U.all <- A1.all - A1.all %*% D.Einv.tD %*% A1.all
trace.UGamma <- sum(U.all * Gamma.all)
attr(trace.UGamma, "h1") <- sum(trace.h1)
attr(trace.UGamma, "h0") <- trace.h0
attr(trace.UGamma, "h1.ndat") <- sum(h1.ndat)
if(Satterthwaite) {
UG <- U.all %*% Gamma.all
trace.UGamma2 <- sum(UG * t(UG))
attr(trace.UGamma, "trace.UGamma2") <- trace.UGamma2
}
# } else {
}
trace.UGamma
}
lav_test_yuan_bentler_mplus_trace <- function(lavsamplestats=NULL,
A1.group = NULL,
B1.group = NULL,
B0.group=NULL,
E.inv=NULL,
meanstructure = TRUE) {
# typical for Mplus:
# - do NOT use the YB formula, but use an approximation
# relying on A0 ~= Delta' A1 Delta and the same for B0
#
# NOTE: if A0 is based on the hessian, then A0 only approximates
# Delta' A1 Delta
#
# - always use h1.information = "unstructured"!!!
ngroups <- lavsamplestats@ngroups
trace.UGamma <- numeric( lavsamplestats@ngroups )
trace.h1 <- numeric( lavsamplestats@ngroups )
trace.h0 <- numeric( lavsamplestats@ngroups )
h1.ndat <- numeric( lavsamplestats@ngroups )
for(g in 1:lavsamplestats@ngroups) {
# group weight
fg <- lavsamplestats@nobs[[g]]/lavsamplestats@ntotal
A1 <- A1.group[[g]]
B1 <- B1.group[[g]]
# mask independent 'fixed-x' variables
zero.idx <- which(diag(A1) == 0)
if(length(zero.idx) > 0L) {
A1.inv <- matrix(0, nrow(A1), ncol(A1))
a1 <- A1[-zero.idx, -zero.idx]
a1.inv <- solve(a1)
A1.inv[-zero.idx, -zero.idx] <- a1.inv
} else {
A1.inv <- solve(A1)
}
h1.ndat[g] <- ncol(A1) - length(zero.idx)
# if data is complete, why not just A1 %*% Gamma?
trace.h1[g] <- sum( B1 * t( A1.inv ) )
trace.h0[g] <- fg * sum( B0.group[[g]] * t(E.inv) )
trace.UGamma[g] <- (trace.h1[g] - trace.h0[g])
}
# we take the sum here
trace.UGamma <- sum(trace.UGamma)
attr(trace.UGamma, "h1") <- trace.h1
attr(trace.UGamma, "h0") <- trace.h0
attr(trace.UGamma, "h1.ndat") <- h1.ndat
trace.UGamma
}
Try the lavaan package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.