# R/find.R In simsem: SIMulated Structural Equation Modeling

#### Documented in findFactorInterceptfindFactorMeanfindFactorResidualVarfindFactorTotalCovfindFactorTotalVarfindIndInterceptfindIndMeanfindIndResidualVarfindIndTotalVarfindPossibleFactorCorfindRecursiveSet

# findFactorIntercept: Find the factor intercept if regression coefficients and
# factor means are specified

findFactorIntercept <- function(beta, factorMean = NULL, gamma = NULL, covmean = NULL) {
if(!is.null(gamma)) {
beta <- parseGammaToBeta(beta, gamma)
factorMean <- c(covmean, factorMean)
}
ni <- nrow(beta)
set <- findRecursiveSet(beta)
intercept <- rep(0, ni)
if (is.null(factorMean))
factorMean <- rep(0, ni)
intercept[set[[1]]] <- factorMean[set[[1]]]
iv <- NULL
iv.mean <- factorMean[set[[1]]]
for (i in seq_len(length(set) - 1)) {
iv <- c(iv, set[[i]])
dv <- set[[i + 1]]
temp.path <- matrix(beta[dv, iv], nrow = length(dv), ncol = length(iv))
mean.reg <- (temp.path %*% iv.mean)
dv.mean <- factorMean[dv]
intercept[dv] <- dv.mean - mean.reg
if (i < (length(set) - 1)) {
agg <- c(iv, dv)
iv.mean <- factorMean[agg]
}
}
if(!is.null(gamma)) {
intercept <- intercept[(length(covmean) + 1):length(intercept)]
}
return(as.vector(intercept))
}

# findFactorResidualVar: Find the factor residual variance if total variances,
# correlation, and regression coefficients are specified.

findFactorResidualVar <- function(beta, corPsi, totalVarPsi = NULL, gamma = NULL, covcov = NULL) {
if(!is.null(gamma)) {
beta <- parseGammaToBeta(beta, gamma)
corPsi <- parseCovCovToPsi(corPsi, cov2cor(covcov))
totalVarPsi <- c(diag(covcov), totalVarPsi)
}
if (sum(diag(corPsi)) == 0)
diag(corPsi) <- 1
ni <- nrow(beta)
set <- findRecursiveSet(beta)
errorVar <- rep(1, ni)
if (is.null(totalVarPsi))
totalVarPsi <- rep(1, ni)
errorVar[set[[1]]] <- totalVarPsi[set[[1]]]
iv <- NULL
ivCor <- corPsi[set[[1]], set[[1]]]
startVar <- totalVarPsi[set[[1]]]
ivCov <- suppressWarnings(cor2cov(as.matrix(ivCor), sqrt(startVar)))
for (i in seq_len(length(set) - 1)) {
iv <- c(iv, set[[i]])
dv <- set[[i + 1]]
tempBeta <- matrix(beta[dv, iv], nrow = length(dv), ncol = length(iv))
var.reg <- (tempBeta %*% ivCov %*% t(tempBeta))
tempPsi <- corPsi[dv, dv]
tempPsiSd <- rep(0, length(dv))
for (j in 1:length(dv)) {
errorVar[dv[j]] <- totalVarPsi[dv[j]] - var.reg[j, j]
if (is.na(errorVar[dv[j]]) || errorVar[dv[j]] < 0) {
tempPsiSd[j] <- NaN
} else {
tempPsiSd[j] <- sqrt(errorVar[dv[j]])
}
}
if (i < (length(set) - 1)) {
tempPsi <- suppressWarnings(cor2cov(as.matrix(tempPsi), tempPsiSd))
real.tempPsi <- matrix(0, length(iv) + length(dv), length(iv) + length(dv))
real.tempPsi[1:length(iv), 1:length(iv)] <- ivCov
real.tempPsi[(length(iv) + 1):(length(iv) + length(dv)), (length(iv) +
1):(length(iv) + length(dv))] <- tempPsi
agg <- c(iv, dv)
blank.path <- matrix(0, nrow = length(iv), ncol = length(agg))
temp.path2 <- beta[dv, agg]
temp.path2 <- rbind(blank.path, temp.path2)
ID <- matrix(0, length(agg), length(agg))
diag(ID) <- 1
ivCov <- solve(ID - temp.path2) %*% real.tempPsi %*% t(solve(ID - temp.path2))
}
}
if(!is.null(gamma)) {
errorVar <- errorVar[(nrow(covcov) + 1):length(errorVar)]
}
return(as.vector(errorVar))
}

# findFactorTotalVar: Find the factor total variance if regression
# coefficients, factor correlation, and factor residual variances are
# specified.

findFactorTotalVar <- function(beta, corPsi, residualVarPsi, gamma = NULL, covcov = NULL) {
if(!is.null(gamma)) {
beta <- parseGammaToBeta(beta, gamma)
corPsi <- parseCovCovToPsi(corPsi, cov2cor(covcov))
residualVarPsi <- c(diag(covcov), residualVarPsi)
}
ni <- nrow(beta)
set <- findRecursiveSet(beta)
real.psi <- suppressWarnings(cor2cov(as.matrix(corPsi), sqrt(residualVarPsi)))
ID <- matrix(0, ni, ni)
diag(ID) <- 1
iv.cov <- solve(ID - beta) %*% real.psi %*% t(solve(ID - beta))
factor.var <- diag(iv.cov)
if(!is.null(gamma)) {
factor.var <- factor.var[(nrow(covcov) + 1):length(factor.var)]
}
return(as.vector(factor.var))
}

# findFactorMean: Find the factor mean if regression coefficients and factor
# intercept are specified.

findFactorMean <- function(beta, alpha = NULL, gamma = NULL, covmean = NULL) {
if(!is.null(gamma)) {
beta <- parseGammaToBeta(beta, gamma)
alpha <- c(covmean, alpha)
}
ni <- nrow(beta)
set <- findRecursiveSet(beta)
factor.mean <- rep(0, ni)
if (is.null(alpha))
alpha <- rep(0, ni)
factor.mean[set[[1]]] <- alpha[set[[1]]]
iv <- NULL
iv.mean <- factor.mean[set[[1]]]
for (i in seq_len(length(set) - 1)) {
iv <- c(iv, set[[i]])
dv <- set[[i + 1]]
temp.path <- matrix(beta[dv, iv], nrow = length(dv), ncol = length(iv))
mean.reg <- (temp.path %*% iv.mean)
factor.mean[dv] <- alpha[dv] + mean.reg
if (i < (length(set) - 1)) {
agg <- c(iv, dv)
iv.mean <- factor.mean[agg]
}
}
if(!is.null(gamma)) {
factor.mean <- factor.mean[(length(covmean) + 1):length(factor.mean)]
}
return(as.vector(factor.mean))
}

# findFactorTotalCov: Find the factor total covariance if regression
# coefficients and factor covariances (which may be made from factor
# correlation, total factor variances, and error factor variances) are
# specified

findFactorTotalCov <- function(beta, psi = NULL, corPsi = NULL, totalVarPsi = NULL,
errorVarPsi = NULL, gamma = NULL, covcov = NULL) {
if (is.null(psi)) {
if (is.null(errorVarPsi))
errorVarPsi <- findFactorResidualVar(beta, corPsi, totalVarPsi)
psi <- suppressWarnings(cor2cov(as.matrix(corPsi), sqrt(errorVarPsi)))
}
iden <- diag(nrow(beta))
temp <- solve(iden - beta)
facTotalCov <- temp %*% psi %*% t(temp)
if(!is.null(gamma)) {
facTotalCov <- facTotalCov + (temp %*% gamma %*% covcov %*% t(gamma) %*% t(temp))
}
return(facTotalCov)
}

# total factor covariance, and measurement error variances.
findIndTotalVar <- function(lambda, totalFactorCov, residualVarTheta, kappa = NULL, covcov = NULL) {
factor.part <- lambda %*% totalFactorCov %*% t(lambda)
indicator.var <- diag(factor.part) + residualVarTheta

if(!is.null(kappa)) indicator.var <- indicator.var + diag(kappa %*% covcov %*% t(kappa))
return(as.vector(indicator.var))
}

# factor covariance, and total indicator variances are specified

findIndIntercept <- function(lambda, factorMean = NULL, indicatorMean = NULL, kappa = NULL, covmean = NULL) {
ni <- nrow(lambda)
nk <- ncol(lambda)
if (is.null(factorMean))
factorMean <- rep(0, nk)
if (is.null(indicatorMean))
indicatorMean <- rep(0, ni)
factor.part <- lambda %*% factorMean
intercept <- indicatorMean - factor.part
if(!is.null(kappa)) intercept <- intercept - (kappa %*% covmean)
return(as.vector(intercept))
}

# findIndResidualVar: Find the residual variances of indicators if factor

findIndResidualVar <- function(lambda, totalFactorCov, totalVarTheta = NULL, kappa = NULL, covcov = NULL) {
ni <- nrow(lambda)
if (is.null(totalVarTheta))
totalVarTheta <- rep(1, ni)
factor.part <- lambda %*% totalFactorCov %*% t(lambda)
error.var <- totalVarTheta - diag(factor.part)
if(!is.null(kappa)) error.var <- error.var - diag(kappa %*% covcov %*% t(kappa))
error.var[(error.var < 0) & (totalVarTheta == 0)] <- 0
return(as.vector(error.var))
}

# findIndMean: Find indicator means based on loading matrix, factor means, and
# measurement intercept.

findIndMean <- function(lambda, factorMean = NULL, tau = NULL, kappa = NULL, covmean = NULL) {
ni <- nrow(lambda)
nk <- ncol(lambda)
if (is.null(factorMean))
factorMean <- rep(0, nk)
if (is.null(tau))
tau <- rep(0, ni)
factor.part <- lambda %*% factorMean
indicator.mean <- tau + factor.part
if(!is.null(kappa)) indicator.mean <- indicator.mean + (kappa %*% covmean)
return(as.vector(indicator.mean))
}

# findPossibleFactorCor: From the set of regression coefficients, this function
# will find the elements that is possible to free covariances or correlations

findPossibleFactorCor <- function(beta) {
ni <- nrow(beta)
set <- findRecursiveSet(beta)
psi <- matrix(0, ni, ni)
diag(psi) <- 1
for (i in 1:length(set)) {
temp.set <- set[[i]]
if (length(temp.set) > 1) {
for (j in 2:length(temp.set)) {
for (k in 1:(j - 1)) {
psi[temp.set[j], temp.set[k]] <- NA
psi[temp.set[k], temp.set[j]] <- NA
}
}
}
}
return(psi)
}

# findRecursiveSet: Group variables together regarding the position in the
# mediation chain

findRecursiveSet <- function(beta) {
result <- list()
ni <- nrow(beta)
fix.variable <- rep(FALSE, ni)
ni.sofar <- 0
i <- 1
while (ni.sofar < ni) {
temp <- findRowZero(beta, fix.variable)
if (is.null(temp))
stop("The matrix is not recursive.")
fix.variable[temp] <- TRUE
result[[i]] <- temp
i <- i + 1
ni.sofar <- ni.sofar + length(temp)
}
return(result)
}

# \title{
# Find rows in a matrix that all elements are zero in non-fixed subset rows and columns.
# }
# \description{
# Find rows in a matrix that all elements are zero in non-fixed subset rows and columns. This function will be used in the \code{\link{findRecursiveSet}} function
# }
# \usage{
# findRowZero(square.matrix, is.row.fixed = FALSE)
# }
# \arguments{
# \item{square.matrix}{
# Any square matrix
# }
# \item{is.row.fixed}{
# A logical vector with the length equal to the dimension of the \code{square.matrix}. If \code{TRUE}, the function will skip examining this row.
# }
# }
# \value{
# A vector of positions that contain rows of all zeros
# }

findRowZero <- function(square.matrix, is.row.fixed = FALSE) {
ni <- nrow(square.matrix)
if (length(is.row.fixed) == 1) {
if (is.row.fixed == FALSE)
is.row.fixed <- rep(FALSE, ni)
}
result <- NULL
desired.zero <- sum(!is.row.fixed)
for (i in 1:ni) {
if (is.row.fixed[i] == FALSE) {
temp <- sum(square.matrix[i, !is.row.fixed] == 0, na.rm = TRUE)
if (temp == desired.zero)
result <- c(result, i)
}
}
return(result)
}


## Try the simsem package in your browser

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

simsem documentation built on March 29, 2021, 1:07 a.m.