# R/Omega_Indices.R In BifactorIndicesCalculator: Bifactor Indices Calculator

#### Documented in cat_Omega_Hcat_Omega_SOmega_HOmega_S

#' Omega_S
#'
#' Computes an omega reliability estimate for all factors as described in Rodriguez, Reise, and
#' Haviland (2016).
#'
#' \code{Omega_S} is called by \code{\link{bifactorIndices}} and the various convenience functions
#' for exploratory models and/or Mplus output,
#' which are the only functions in this package intended for casual users.
#'
#' @param Lambda is a matrix of factor loadings
#' @param Theta is a vector of indicator error variances
#'
#' @return A \code{numeric}, the omega reliability estimate for all factors.
#'
#' @examples
#' Lambda <- matrix(c(.82, .10,   0,   0,
#'                    .77, .35,   0,   0,
#'                    .79, .32,   0,   0,
#'                    .66, .39,   0,   0,
#'                    .51,   0, .71,   0,
#'                    .56,   0, .43,   0,
#'                    .68,   0, .13,   0,
#'                    .60,   0, .50,   0,
#'                    .83,   0,   0, .47,
#'                    .60,   0,   0, .27,
#'                    .78,   0,   0, .28,
#'                    .55,   0,   0, .75),
#'                    ncol = 4, byrow = TRUE)
#' colnames(Lambda) <- c("General", "SF1", "SF2", "SF3")
#' Theta <- rep(1, nrow(Lambda)) - rowSums(Lambda^2)
#' Omega_S(Lambda, Theta)
#'
#' @references
#' Rodriguez, A., Reise, S. P., & Haviland, M. G. (2016). Evaluating bifactor models:
#' calculating and interpreting statistical indices. \emph{Psychological Methods, 21}(2),
#' 137 \doi{10.1037/met0000045}.
#'
#' @export
#'
#'

Omega_S <- function(Lambda, Theta) {
Omega_S_C <- function(Fac, Lambda, Theta) {
## Make a matrix of logical vectors for non-zero elements of Lambda.
inFactor <- Lambda[,Fac] != 0
## Compute the appropriate ratio of sums
sum(colSums(Lambda*inFactor)^2)/(sum(colSums(Lambda*inFactor)^2) + sum(Theta*inFactor))
}
if (is.null(Theta)) return(NULL)
omega_results <- sapply(1:ncol(Lambda), Omega_S_C, Lambda = Lambda, Theta = Theta)
names(omega_results) <- colnames(Lambda)
omega_results
}

#' cat_Omega_S
#'
#' Computes an omega reliability estimate for all factors as described in Rodriguez, Reise, and
#' Haviland (2016).
#'
#' \code{cat_Omega_S} is called by \code{\link{bifactorIndices}} and the various convenience functions
#' for exploratory models and/or Mplus output,
#' which are the only functions in this package intended for casual users.
#'
#' @param Lambda is a matrix of standardized factor loadings
#' @param Thresh is a list (indexed by items) of vectors of item thresholds (items must be
#' on a standardized metric).
#' @param Phi is the latent variable covariance matrix. Defaults to \code{NULL}, and
#' the identity matrix will be used. No other options are currently available.
#' @param Denom specifies how the variance of the total score will be computed. Defaults
#' to \code{NULL}, and the model implied total score variance will be used. No other options
#' are currently available.
#'
#' @return A \code{numeric}, the omega reliability estimate for all factors using the technique of
#' Green and Yang (2009).
#'
#' @examples
#'
#' Lambda <- matrix(c(.82, .10,   0,   0,
#'                    .77, .35,   0,   0,
#'                    .79, .32,   0,   0,
#'                    .66, .39,   0,   0,
#'                    .51,   0, .71,   0,
#'                    .56,   0, .43,   0,
#'                    .68,   0, .13,   0,
#'                    .60,   0, .50,   0,
#'                    .83,   0,   0, .47,
#'                    .60,   0,   0, .27,
#'                    .78,   0,   0, .28,
#'                    .55,   0,   0, .75),
#'                    ncol = 4, byrow = TRUE)
#' colnames(Lambda) <- c("General", "SF1", "SF2", "SF3")
#'
#' Thresh = list(c(-1, 0, 1),  c(-0.5, 0, 0.5),
#'               c(0, 1, 2),   c(0, 0.5, 1),
#'               c(-2, -1, 0), c(-1, -0.5, 0),
#'               c(-1, 0, 2),  c(-0.5, 0, 1),
#'               c(-2, 0, 1),  c(-1, 0, 0.5),
#'               c(-1, 0, 1),  c(-0.5, 0, 0.5))
#'
#' cat_Omega_S(Lambda, Thresh)
#'
#' @references
#' Rodriguez, A., Reise, S. P., & Haviland, M. G. (2016). Evaluating bifactor models:
#' calculating and interpreting statistical indices. \emph{Psychological Methods, 21}(2),
#' 137 \doi{10.1037/met0000045}.
#'
#' Green, S. B., & Yang, Y. (2009). Reliability of summed item scores using
#' structural equation modeling: An alternative to coefficient alpha.
#' \emph{Psychometrika, 74}(1), 155-167 \doi{10.1007/s11336-008-9099-3}.
#'
#' @export
#'
#'
#'

cat_Omega_S <- function(Lambda, Thresh, Phi = NULL, Denom = NULL) {

cat_Omega_S_C <- function(Fac, Lambda, Thresh, Phi = NULL, Denom = NULL) {
## Make a matrix of logical vectors for non-zero elements of Lambda.
inFactor <- Lambda[,Fac] != 0
## subset Lambda to only include appropriate items
## and make sure it is still a MATRIX!!
subLambda <- Lambda[inFactor,]
subLambda <- as.matrix(subLambda)
## create Phi matrix
if (is.null(Phi)) {Phi <- diag(ncol(subLambda))}

num_items <- nrow(subLambda)

## latent item covariances (correlations, because they're standardized!!)
lat_cov <- subLambda %*% Phi %*% t(subLambda)
if (is.null(Denom)) {
poly_cor <- lat_cov
diag(poly_cor) <- rep(1, num_items)
}

## Create covariance matrix of parallel items
par_cov_mat <- sapply(1:num_items, function (j) {
t_j <- Thresh[[j]]
sapply(1:num_items, function(jp) {
t_jp <- Thresh[[jp]]
## cov(x_j, x_jp)

## first, the left half of the expression in equation 19 in Green and Yang (2009)
left <- sum(sapply(1:length(t_j), function(c) {
sapply(1:length(t_jp), function(cp) {
r <- lat_cov[j, jp]
mnormt::pmnorm(c(t_j[c], t_jp[cp]), c(0, 0), matrix(c(1, r, r, 1), 2))
}) ## end cp
})) ## end c

## now the two expression in the right half of the expression in equation 19 in G&Y
right_j  <- sum(stats::pnorm(t_j))
right_jp <- sum(stats::pnorm(t_jp))

## put it together and what do you get? Bibbidi-Bobbidi-Boo
left - right_j * right_jp

}) ## end jp
}) ## end j

## create covariance matrix of items... copy/paste of par_cov_mat,
## but switch from lat_cov_mat to poly_cor
item_cov_mat <- sapply(1:num_items, function (j) {
t_j <- Thresh[[j]]
sapply(1:num_items, function(jp) {
t_jp <- Thresh[[jp]]
## cov(x_j, x_jp)

## first, the left half of the expression in equation 19 in Green and Yang (2009)
left <- sum(sapply(1:length(t_j), function(c) {
sapply(1:length(t_jp), function(cp) {
r <- poly_cor[j, jp]
mnormt::pmnorm(c(t_j[c], t_jp[cp]), c(0, 0), matrix(c(1, r, r, 1), 2))
}) ## end cp
})) ## end c

## now the two expression in the right half of the expression in equation 19 in G&Y
right_j  <- sum(stats::pnorm(t_j))
right_jp <- sum(stats::pnorm(t_jp))

## put it together and what do you get? Bibbidi-Bobbidi-Boo
left - right_j * right_jp

}) ## end jp
}) ## end j

numer <- sum(par_cov_mat)
denom <- sum(item_cov_mat)

numer/denom
}

omega_results <- sapply(1:ncol(Lambda), cat_Omega_S_C, Lambda = Lambda, Thresh = Thresh)
names(omega_results) <- colnames(Lambda)
omega_results
}

#' OmegaH
#'
#' Computes hierarchical omega reliability estimate for all factors as described in
#' Rodriguez, Reise, and Haviland (2016).
#'
#' \code{Omega_H} is called by \code{\link{bifactorIndices}} and the various convenience functions
#' for exploratory models and/or Mplus output,
#' which are the only functions in this package intended for casual users.
#'
#' @param Lambda is a matrix of factor loadings
#' @param Theta is a vector of indicator error variances
#'
#' @return A \code{numeric}, the omega reliability estimate for all factors.
#'
#' @examples
#' Lambda <- matrix(c(.82, .10,   0,   0,
#'                    .77, .35,   0,   0,
#'                    .79, .32,   0,   0,
#'                    .66, .39,   0,   0,
#'                    .51,   0, .71,   0,
#'                    .56,   0, .43,   0,
#'                    .68,   0, .13,   0,
#'                    .60,   0, .50,   0,
#'                    .83,   0,   0, .47,
#'                    .60,   0,   0, .27,
#'                    .78,   0,   0, .28,
#'                    .55,   0,   0, .75),
#'                    ncol = 4, byrow = TRUE)
#' colnames(Lambda) <- c("General", "SF1", "SF2", "SF3")
#' Theta <- rep(1, nrow(Lambda)) - rowSums(Lambda^2)
#' Omega_H(Lambda, Theta)
#'
#' @section References:
#' Rodriguez, A., Reise, S. P., & Haviland, M. G. (2016). Evaluating bifactor models:
#' Calculating and interpreting statistical indices. Psychological Methods, 21(2),
#' 137 \doi{10.1037/met0000045}.
#'
#' @export
#'
#'
#'

Omega_H <- function(Lambda, Theta) {
Omega_H_C <- function(Fac, Lambda, Theta) {
## Make a matrix of logical vectors for non-zero elements of Lambda.
inFactor <- Lambda[,Fac] != 0
## Compute the appropriate ratio of sums
sum(Lambda[,Fac])^2/(sum(colSums(Lambda*inFactor)^2) + sum(Theta*inFactor))
}
if (is.null(Theta)) return(NULL)
omega_results <- sapply(1:ncol(Lambda), Omega_H_C, Lambda = Lambda, Theta = Theta)
names(omega_results) <- colnames(Lambda)
omega_results
}

#' cat_Omega_H
#'
#' Computes hierarchical omega reliability estimate for all factors as described in Rodriguez, Reise, and
#' Haviland (2016).
#'
#' \code{cat_Omega_H} is called by \code{\link{bifactorIndices}} and the various convenience functions
#' for exploratory models and/or Mplus output,
#' which are the only functions in this package intended for casual users.
#'
#' @param Lambda is a matrix of standardized factor loadings
#' @param Thresh is a list (indexed by items) of vectors of item thresholds
#' @param Phi is the latent variable covariance matrix. Defaults to \code{NULL}, and
#' the identity matrix will be used. No other options are currently available.
#' @param Denom specifies how the variance of the total score will be computed. Defaults
#' to \code{NULL}, and the model implied total score variance will be used. No other options
#' are currently available.
#'
#' @return A \code{numeric}, the hierarchical omega reliability estimate for all factors using
#' the technique of Green and Yang (2009).
#'
#' @examples
#' Lambda <- matrix(c(.82, .10,   0,   0,
#'                    .77, .35,   0,   0,
#'                    .79, .32,   0,   0,
#'                    .66, .39,   0,   0,
#'                    .51,   0, .71,   0,
#'                    .56,   0, .43,   0,
#'                    .68,   0, .13,   0,
#'                    .60,   0, .50,   0,
#'                    .83,   0,   0, .47,
#'                    .60,   0,   0, .27,
#'                    .78,   0,   0, .28,
#'                    .55,   0,   0, .75),
#'                    ncol = 4, byrow = TRUE)
#' colnames(Lambda) <- c("General", "SF1", "SF2", "SF3")
#'
#' Thresh = list(c(-1, 0, 1),  c(-0.5, 0, 0.5),
#'               c(0, 1, 2),   c(0, 0.5, 1),
#'               c(-2, -1, 0), c(-1, -0.5, 0),
#'               c(-1, 0, 2),  c(-0.5, 0, 1),
#'               c(-2, 0, 1),  c(-1, 0, 0.5),
#'               c(-1, 0, 1),  c(-0.5, 0, 0.5))
#'
#' cat_Omega_H(Lambda, Thresh)
#'
#' @references
#' Rodriguez, A., Reise, S. P., & Haviland, M. G. (2016). Evaluating bifactor models:
#' calculating and interpreting statistical indices. \emph{Psychological Methods, 21}(2),
#' 137 \doi{10.1037/met0000045}.
#'
#' Green, S. B., & Yang, Y. (2009). Reliability of summed item scores using
#' structural equation modeling: An alternative to coefficient alpha.
#' \emph{Psychometrika, 74}(1), 155-167 \doi{10.1007/s11336-008-9099-3}.
#'
#' @export
#'
#'
#'

cat_Omega_H <- function(Lambda, Thresh, Phi = NULL, Denom = NULL) {
cat_Omega_H_C <- function(Fac, Lambda, Thresh, Phi = NULL, Denom = NULL) {
## Make a matrix of logical vectors for non-zero elements of Lambda.
inFactor <- Lambda[,Fac] != 0
## subset Lambda to only include appropriate items
subLambda <- Lambda[inFactor,]
subLambda <- as.matrix(subLambda)

## create Phi matrix
if (is.null(Phi)) {Phi <- diag(ncol(subLambda))}

num_items <- nrow(subLambda)

## latent item covariances (correlations, because they're standardized!!)
lat_cov <- subLambda %*% Phi %*% t(subLambda)
if (is.null(Denom)) {
poly_cor <- lat_cov
diag(poly_cor) <- rep(1, num_items)
}

# Now that we have the poly_cor matrix, we need to restrict to just the factor of interest
#for computing the numerator

subLambda <- subLambda[, Fac]
subLambda <- as.matrix(subLambda)

Phi <- Phi[Fac, Fac]
lat_cov <- subLambda %*% as.matrix(Phi) %*% t(subLambda)

## Create covariance matrix of parallel items
par_cov_mat <- sapply(1:num_items, function (j) {
t_j <- Thresh[[j]]
sapply(1:num_items, function(jp) {
t_jp <- Thresh[[jp]]
## cov(x_j, x_jp)

## first, the left half of the expression in equation 19 in Green and Yang (2009)
left <- sum(sapply(1:length(t_j), function(c) {
sapply(1:length(t_jp), function(cp) {
r <- lat_cov[j, jp]
mnormt::pmnorm(c(t_j[c], t_jp[cp]), c(0, 0), matrix(c(1, r, r, 1), 2))
}) ## end cp
})) ## end c

## now the two expression in the right half of the expression in equation 19 in G&Y
right_j  <- sum(stats::pnorm(t_j))
right_jp <- sum(stats::pnorm(t_jp))

## put it together and what do you get? Bibbidi-Bobbidi-Boo
left - right_j * right_jp

}) ## end jp
}) ## end j

## create covariance matrix of items... copy/paste of par_cov_mat,
## but switch from lat_cov_mat to poly_cor
item_cov_mat <- sapply(1:num_items, function (j) {
t_j <- Thresh[[j]]
sapply(1:num_items, function(jp) {
t_jp <- Thresh[[jp]]
## cov(x_j, x_jp)

## first, the left half of the expression in equation 19 in Green and Yang (2009)
left <- sum(sapply(1:length(t_j), function(c) {
sapply(1:length(t_jp), function(cp) {
r <- poly_cor[j, jp]
mnormt::pmnorm(c(t_j[c], t_jp[cp]), c(0, 0), matrix(c(1, r, r, 1), 2))
}) ## end cp
})) ## end c

## now the two expression in the right half of the expression in equation 19 in G&Y
right_j  <- sum(stats::pnorm(t_j))
right_jp <- sum(stats::pnorm(t_jp))

## put it together and what do you get? Bibbidi-Bobbidi-Boo
left - right_j * right_jp

}) ## end jp
}) ## end j

numer <- sum(par_cov_mat)
denom <- sum(item_cov_mat)

numer/denom
}

omega_results <- sapply(1:ncol(Lambda), cat_Omega_H_C, Lambda = Lambda, Thresh = Thresh)
names(omega_results) <- colnames(Lambda)
omega_results
}


## Try the BifactorIndicesCalculator package in your browser

Any scripts or data that you put into this service are public.

BifactorIndicesCalculator documentation built on May 13, 2021, 1:08 a.m.