Nothing
#require(MASS)
#require(lavaan)
#require(numDeriv)
## Mirror the lower or upper diagonal of a square matrix
mirror.tri <- function(mat, upper.on.lower = FALSE) {
if (upper.on.lower) {
mat[lower.tri(mat)] <- t(mat)[lower.tri(mat)]
} else {
mat[upper.tri(mat)] <- t(mat)[upper.tri(mat)]
}
mat
}
## Factor analysis using principal components as the extraction method
factanal.prcomp <- function (covmat, nfactors = NULL) {
prc <- princomp(covmat = covmat, cor = TRUE)
eig <- prc$sdev^2
if (is.null(nfactors)) nfactors <- sum(eig >= 1)
loads <- prc$loadings[, 1:nfactors]
loads[, 1:nfactors] %*% diag(prc$sdev[1:nfactors])
}
## Get a covariance matrix from a correlation matrix and std deviations
cor2cov <- function(cormat, sds) {
outer(sds, sds, "*") * cormat
}
## Check if a square matrix is Hermitian
is.hermitian <- function(covmat) {
t.covmat <- t(covmat)
if (all(dim(covmat) == dim(t.covmat)) && all(covmat == t(covmat))) {
return(TRUE)
} else {
return(FALSE)
}
}
## Check if all covariances are within bounds
is.within.bounds <- function(covmat) {
if (!is.hermitian(covmat)) stop("Matrix is not Hermitian")
sds <- sqrt(diag(covmat))
max.cov <- outer(sds, sds, "*")
differences <- round(max.cov - covmat, 3)
res <- all(differences >= 0)
if (!res) {
outside.bounds <- which(differences < 0, arr.ind = TRUE)
varnames <- rownames(differences)
cat("Out of bounds covariances:\n\n")
if (!is.null(varnames)) {
for (i in 1:(NROW(outside.bounds)/2)) {
rn <- varnames[outside.bounds[i, 1]]
cn <- varnames[outside.bounds[i, 2]]
cat(sprintf("%s and %s: max - %f, obs - %f", rn,
cn, max.cov[rn, cn], covmat[rn, cn]), "\n")
}
} else {
for (i in 1:(NROW(outside.bounds)/2)) {
rn <- outside.bounds[i, 1]
cn <- outside.bounds[i, 2]
cat(sprintf("row %s and column %s: max - %f, obs - %f", rn,
cn, max.cov[rn, cn], covmat[rn, cn]), "\n")
}
}
cat("\n")
}
res
}
## Check if a covariance matrix is positive-definite
## We could try to invert it, but instead we opted to do a bunch of
## testing for learning purposes.
is.pd <- function(covmat) {
## Check if matrix is Hermitian (lower and upper diag are equal)
is.ht <- is.hermitian(covmat)
if (!is.ht) {
warning("Matrix isn't Hermitian")
return(FALSE)
}
## Check if all covariances are within bounds
is.wd <- is.within.bounds(covmat)
if (!is.wd) {
warning("Covariances aren't within bounds")
return(FALSE)
}
## Check if the matrix is invertible (this whole funcion could be
## replaced by a call to chol
is.invertible <- try(solve(covmat), silent = TRUE)
if (class(is.invertible) == "try-error") {
warning("Matrix isn't invertible")
return(FALSE)
}
## Check if all eigenvalues are positive
eigenvalues <- eigen(covmat)$values
if (!all(eigenvalues > 0)) {
warning("Matrix has negative eigenvalues")
return(FALSE)
}
## Check if the determinant is higher than 0
determ <- prod(eigenvalues)
if (determ <= 0) {
warning("The determinant is equal or lower than 0")
return(FALSE)
}
TRUE
}
## Creates a multivariate normal distribution from a covariance matrix.
rmnorm <- function(covmat, means = 0, sds = NULL, n) {
if (is.null(sds))
sds <- sqrt(diag(covmat))
pattern <- factanal.prcomp(covmat = covmat, nfactors = NROW(covmat))
nvar <- NROW(covmat)
X <- sapply(seq_len(nvar), function(x) {rnorm(n)})
X <- t(X)
Z <- pattern %*% X
Z <- t(Z)
Z <- sweep(Z, 2, sds, "*")
Z <- sweep(Z, 2, means, "+")
Z
}
## Identify the type of a lavaanified model
model.type <- function(model) {
if(all(model$op != "=~"))
type <- "pa"
if(any(model$op == "=~") & any(model$op == "~"))
type <- "sem"
if(any(model$op == "=~") & all(model$op != "~"))
type <- "cfa"
type
}
## Build but dont fit a model using lavaan to use as a starting point.
lvnfy <- function(model, data, n = NULL, ...) {
if (is.hermitian(data)) {
out <- sem(model, sample.cov = data, sample.nobs = n, do.fit
= FALSE, ...)
} else {
out <- sem(model, data = data, do.fit = FALSE, ...)
}
lv.list <- inspect(out, "list")
lv.list$ustart <-out@Fit@start
lv.list <- as.list(lv.list[c(1:4,7:8)])
vars <- variables(lv.list)
var.names <- rownames(vars)
latents <- var.names[!vars[,"lat"]]
regressed <- !vars[,"lat"] & (vars[,"endo"] | vars[,"exo"])
vars <- cbind(vars, "latenterise" = regressed)
attr(lv.list, "vars") <- vars
lv.list
}
## Get latent variable names
latents <- function(model) {
vars <- attr(model, "vars")
rownames(vars)[vars[,"lat"]]
}
## Get observed variables (even if they are not indicators)
observed <- function(model) {
vars <- attr(model, "vars")
rownames(vars)[!vars[,"lat"]]
}
## Get indicators and separate them by type
manifests <- function(model) {
vars <- attr(model, "vars")
rownames(vars)[vars[,"ind"]]
}
## Get endogenous variables
endogenous <- function(model) {
vars <- attr(model, "vars")
rownames(vars)[vars[,"endo"]]
}
## Get endogenous variables
exogenous <- function(model) {
vars <- attr(model, "vars")
rownames(vars)[vars[,"exo"]]
}
## Identity (help done)
II <- function(model) {
II <- BETA(model)
II[] <- 0
diag(II) <- 1
II
}
## Beta matrix (regression coefficients between endogenous variables)
BETA <- function(model) {
vars <- attr(model, "vars")
var.names <- rownames(vars)
## latents
latents <- vars[,"lat"]
obs.not.indicators <- !vars[,"lat"] & (vars[,"exo"] | vars[,"endo"])
sides <- var.names[latents | obs.not.indicators]
regressions <- lapply(model, `[`, model$op == "~" | (model$op ==
"=~" &
model$rhs %in%
var.names[latents]))
sides <- unique(var.names[latents | obs.not.indicators])
beta.mat <- matrix(0, length(sides), length(sides))
dimnames(beta.mat) <- list(sides, sides)
free.pars.pos <- NULL
index.matrix <- beta.mat
index.matrix[] <- seq_along(index.matrix)
free.pars <- regressions$free[regressions$free > 0L]
if (length(regressions) > 0) {
for (i in seq_len(length(regressions$lhs))) {
lhs <- regressions$lhs[i]
rhs <- regressions$rhs[i]
op <- regressions$op[i]
start.value <- regressions$ustart[i]
is.free <- regressions$free[i] > 0L
if (op == "~") {
beta.mat[lhs, rhs] <- start.value
if (is.free) {
free.pars.pos <- c(free.pars.pos, index.matrix[lhs,
rhs])
}
} else{
beta.mat[rhs, lhs] <- start.value
if (is.free) {
free.pars.pos <- c(free.pars.pos, index.matrix[rhs,
lhs])
}
}
}
}
attr(beta.mat, "free.pars") <- free.pars
attr(beta.mat, "free.pars.pos") <- free.pars.pos
beta.mat
}
## Psi matrix (variances and covariances of observed variables).
PSI <- function(model) {
vars <- attr(model, "vars")
var.names <- rownames(vars)
latents <- vars[,"lat"]
obs.not.indicators <- !vars[,"lat"] & (vars[,"exo"] | vars[,"endo"])
sides <- var.names[latents | obs.not.indicators]
psi.mat <- matrix(0, length(sides), length(sides))
dimnames(psi.mat) <- list(sides, sides)
var.covar <- lapply(model, `[`, model$op == "~~" & model$lhs %in%
sides & model$rhs %in% sides)
#free.pars.pos <- numeric(length(var.covar$lhs) * 2 - NCOL(psi.mat))
#free.pars <- numeric(length(free.pars.pos))
free.pars.pos <- free.pars <- NULL
index.matrix <- psi.mat
index.matrix[] <- seq_along(index.matrix)
if (length(var.covar) == 0) stop("Looks like you forgot to explicitly
assign at least the disturbances in the model.")
for (i in seq_len(length(var.covar$lhs))) {
lhs <- var.covar$lhs[i]
rhs <- var.covar$rhs[i]
start.value <- var.covar$ustart[i]
pos <- var.covar$free[i]
is.free <- pos > 0L
if (lhs == rhs) {
psi.mat[lhs, rhs] <- start.value
if (is.free) {
free.pars.pos <- c(free.pars.pos, index.matrix[lhs, rhs])
free.pars <- c(free.pars, pos)
}
} else {
psi.mat[lhs, rhs] <- psi.mat[rhs, lhs] <- start.value
if (is.free) {
free.pars.pos <- c(free.pars.pos, index.matrix[lhs, rhs],
index.matrix[rhs, lhs])
free.pars <- c(free.pars, rep(pos, 2))
}
}
}
if (length(psi.mat) == 0L) {
return(0)
}
attr(psi.mat, "free.pars") <- free.pars
attr(psi.mat, "free.pars.pos") <- free.pars.pos
psi.mat
}
## Coefficients between latent endogenous and observed
LAMBDA.Y <- function(model) {
vars <- attr(model, "vars")
var.names <- rownames(vars)
observed <- !vars[,"lat"]
latents <- vars[,"lat"] | (!vars[,"lat"] & (vars[,"exo"] |
vars[,"endo"]))
loadings.y <- lapply(model, `[`, model$op == "=~" & !model$rhs %in%
var.names[vars[,"lat"]])
lambda.y.mat <- matrix(0, ncol = sum(latents), nrow =
sum(observed))
colnames(lambda.y.mat) <- var.names[latents]
rownames(lambda.y.mat) <- var.names[observed]
free.pars.pos <- NULL
free.pars <- loadings.y$free[loadings.y$free > 0L]
index.matrix <- lambda.y.mat
index.matrix[] <- seq_along(index.matrix)
latenterised <- var.names[!vars[,"lat"] & (vars[,"exo"] |
vars[,"endo"])]
for (i in latenterised) {
lambda.y.mat[i, i] <- 1
}
for (i in seq_len(length(loadings.y$lhs))) {
lhs <- loadings.y$lhs[i]
rhs <- loadings.y$rhs[i]
start.value <- loadings.y$ustart[i]
is.free <- loadings.y$free[i] > 0L
lambda.y.mat[rhs, lhs] <-
loadings.y$ustart[i]
if (is.free) {
free.pars.pos <- c(free.pars.pos, index.matrix[rhs, lhs])
}
}
if (length(lambda.y.mat) == 0L) {
return(II(model))
free.pars.pos <- NULL
free.pars <- NULL
}
attr(lambda.y.mat, "free.pars") <- free.pars
attr(lambda.y.mat, "free.pars.pos") <- free.pars.pos
lambda.y.mat
}
## Variance/covariance of measurement errors of observed endogenous
THETA.EPSILON <- function(model) {
vars <- attr(model, "vars")
var.names <- rownames(vars)
observed <- !vars[,"lat"]
indicators <- !vars[,"lat"] & (!vars[,"exo"] &
!vars[,"endo"])
te <- matrix(0, sum(observed), sum(observed))
obs.names <- var.names[observed]
ind.names <- var.names[indicators]
dimnames(te) <- list(obs.names, obs.names)
var.covar <- lapply(model, `[`, model$op == "~~" & model$lhs
%in% ind.names)
free.pars.pos <- free.pars <- NULL
index.matrix <- te
index.matrix[] <- seq_along(index.matrix)
if (length(te) == 0L) {
return(matrix(0))
}
for (i in seq_len(length(var.covar$lhs))) {
lhs <- var.covar$lhs[i]
rhs <- var.covar$rhs[i]
start.value <- var.covar$ustart[i]
pos <- var.covar$free[i]
is.free <- pos > 0L
if (lhs == rhs) {
te[lhs, rhs] <- start.value
if (is.free) {
free.pars.pos <- c(free.pars.pos, index.matrix[lhs,
rhs])
free.pars <- c(free.pars, pos)
}
} else {
te[lhs, rhs] <- te[rhs, lhs] <- start.value
if (is.free) {
free.pars.pos <- c(free.pars.pos, index.matrix[lhs,
rhs], index.matrix[rhs, lhs])
free.pars <- c(free.pars, rep(pos, 2))
}
}
}
attr(te, "free.pars") <- free.pars
attr(te, "free.pars.pos") <- free.pars.pos
te
}
## Get the number of degrees of freedom of a given model
model.df <- function(model) {
## t-rule
## consider equality constraints eventually
nexo <- length(exogenous(model))
nvar <- length(observed(model))
npar <- max(model$free)
model <- lapply(model, `[`, model$op == "~~" & model$free == 0 & model$lhs == model$rhs)
nfix <- length(unique(c(model$lhs, model$rhs)))
ndf <- nvar*(nvar + 1)/2 - npar - nfix*(nfix + 1)/2
list(free = npar,
fixed = nfix,
total.df = ndf + npar,
model.df = ndf)
}
## Function to extract the starting values of the free parameters of a
## model
jumpstart <- function(model) {
model$ustart[model$free > 0L]
}
## Function to be passed to the optimizer to minimize the log likelihood
iterator <- function(pars, model, sigma, matrep, fun = "fml") {
model$ustart[model$free > 0L] <- pars
sigma.hat <- sem.matrices(model, matrep)
do.call(fun, list(sigma, sigma.hat))
}
## This function takes the 4 basic matrices and returns a covariance
## matrix based on the model parameters. It can be used to generate a
## 'population' cov matrix for testing purposes or to generate sigmas
## for optimization purposes for instance.
sem.matrices <- function(model, matrep) {
if (!is.null(model$est)) {
model$ustart <- model$est
}
if (any(is.na(model$ustart))) stop("Please give starting values for
all coefficients in the model.")
matrep$BE[attr(matrep$BE, "free.pars.pos")] <-
model$ustart[match(attr(matrep$BE, "free.pars"), model$free)]
matrep$PS[attr(matrep$PS, "free.pars.pos")] <-
model$ustart[match(attr(matrep$PS, "free.pars"), model$free)]
matrep$LY[attr(matrep$LY, "free.pars.pos")] <-
model$ustart[match(attr(matrep$LY, "free.pars"), model$free)]
matrep$TE[attr(matrep$TE, "free.pars.pos")] <-
model$ustart[match(attr(matrep$TE, "free.pars"), model$free)]
matrep$LY %*% solve(matrep$Inv - matrep$BE) %*% matrep$PS %*%
solve(t(matrep$Inv - matrep$BE)) %*% t(matrep$LY) + matrep$TE
}
## Maximum-likelihood fitting function.
fml <- function(sigma, sigma.hat) {
sigma.hat.inv <- try(chol2inv(chol(sigma.hat)), silent = TRUE)
if (class(sigma.hat.inv) == "try-error") return(Inf)
res <- (log(det(sigma.hat)) +
sum(sigma * sigma.hat.inv) -
log(det(sigma)) -
NROW(sigma.hat))
res
}
## Generalised least squares.
fgls <- function(sigma, sigma.hat) {
0.5 * sum(diag((solve(sigma) %*% (sigma - sigma.hat))^2))
}
## Function that calls the optimiser and builds a response
#fit.nlm <- function(model, sigma, n, fun = "fml") {
#model <- lvnfy(model, data = sigma, n = n)
#pars <- model$ustart[model$free > 0L]
#SCALE <- ifelse(abs(pars) > 1.00, abs(1/pars), 1)
#matrep <- matrix.representation(model)
#rnames <- rownames(sigma.hat(matrep))
#sigma <- sigma[rnames, rnames]
#res <- nlm(p = pars,
#f = iterator,
#hessian = TRUE,
#check.analyticals = TRUE,
#model = model,
#fun = fun,
#fscale = SCALE,
#sigma = sigma,
#matrep = matrep)
#hess <- res$hessian
#std.errors <- se.sem(hess, n)
#if (!std.errors$trust) warning("Could not compute standard errors.")
#model$est <- model$ustart
#model$est[model$free > 0L] <- res$estimate
#model$se <- rep(0, length(model$id))
#model$se[model$free > 0L] <- std.errors$values
#model$z <- model$est / model$se
#model$z <- ifelse(model$z == Inf, NA, model$z)
#model$p <- round(do.call("pval", list(model$z)), 3)
#res$model <- data.frame(model, stringsAsFactors = FALSE)
#res$chisquared <- n * res$minimum
#res$df <- model.df(model)$model.df
#res$n <- n
#res$pval <- 1 - pchisq(res$chisquared, res$df)
#res$estimator <- fun
#res$sigma.hat <- sem.matrices(model, matrep)
#res$sigma <- sigma
#class(res) <- c("sem", "list")
#invisible(res)
#}
#fit.optim <- function(model, sigma, n, fun = "fml") {
#model <- lvnfy(model, data = sigma, n = n)
#pars <- model$ustart[model$free > 0L]
#SCALE <- ifelse(abs(pars) > 1.00, abs(1.00/pars), 1.00)
#matrep <- matrix.representation(model)
#rnames <- rownames(sigma.hat(matrep))
#sigma <- sigma[rnames, rnames]
#res <- optim(par = pars,
#fn = iterator,
#gr = gradit,
#hessian = TRUE,
#model = model,
#fun = fun,
#sigma = sigma,
#matrep = matrep,
#method = "L-BFGS-B")
#hess <- res$hessian
#std.errors <- se.sem(hess, n)
#model$est <- model$ustart
#model$est[model$free > 0L] <- res$par
#model$se <- rep(0, length(model$id))
#model$se[model$free > 0L] <- std.errors
#model$z <- model$est / model$se
#model$z <- ifelse(model$z == Inf, NA, model$z)
#model$p <- round(do.call("pval", list(model$z)), 3)
#res$model <- data.frame(model, stringsAsFactors = FALSE)
#res$chisquared <- n * res$value
#res$df <- model.df(model)$model.df
#res$n <- n
#res$pval <- 1 - pchisq(res$chisquared, res$df)
#res$estimator <- fun
#res$sigma.hat <- sem.matrices(model, matrep)
#res$sigma <- sigma
#class(res) <- c("sem", "list")
#res
#}
#fit.optimx <- function(model, sigma, n, fun = "fml") {
#model <- lvnfy(model, data = sigma, n = n)
#pars <- model$ustart[model$free > 0L]
#SCALE <- ifelse(pars > 1.00, abs(pars/1.00), 1.00)
#matrep <- matrix.representation(model)
#rnames <- rownames(sigma.hat(matrep))
#sigma <- sigma[rnames, rnames]
#res <- optimx(par = pars,
#fn = iterator,
#gr = gradit,
#hessian = TRUE,
#model = model,
#fun = fun,
#sigma = sigma,
#matrep = matrep,
#control = list(factr = 1e-14, all.methods = TRUE),
#method = "L-BFGS-B")
#hess <- res$hessian
#std.errors <- se.sem(hess, n)
#model$est <- model$ustart
#model$est[model$free > 0L] <- res$par
#model$se <- rep(0, length(model$id))
#model$se[model$free > 0L] <- std.errors
#model$z <- model$est / model$se
#model$z <- ifelse(model$z == Inf, NA, model$z)
#model$p <- round(do.call("pval", list(model$z)), 3)
#res$model <- data.frame(model, stringsAsFactors = FALSE)
#res$chisquared <- n * res$value
#res$df <- model.df(model)$model.df
#res$n <- n
#res$pval <- 1 - pchisq(res$chisquared, res$df)
#res$estimator <- fun
#res$sigma.hat <- sem.matrices(model, matrep)
#res$sigma <- sigma
#class(res) <- c("sem", "list")
#res
#}
## Using stats::nlminb to try and find the best set of parameter
## estimates for a given model. We tried a bunch of other functions, but
## nlminb seems to be the most reliable.
fit.nlminb <- function(model, sigma, n, fun = "fml") {
## Setting things up
model <- lvnfy(model, data = sigma, n = n)
pars <- model$ustart[model$free > 0L]
## Setting the scale (got the idea from lavaan)
SCALE <- ifelse(abs(pars) > 1.00, abs(1.00/pars), 1)
matrep <- matrix.representation(model)
rnames <- rownames(sigma.hat(matrep))
## Making sure sigma follows sigma.hat
sigma <- sigma[rnames, rnames]
## Optimizer (finding a set of parameter estimates)
res <- nlminb(start = pars,
objective = iterator,
gradient = first.derivative.param,
## Using numDeriv::grad works as well
#gradient = gradit,
model = model,
fun = fun,
scale = SCALE,
sigma = sigma,
matrep = matrep,
control = list(iter.max = 10000, eval.max = 20000))
## Calculating the Hessian (standard errors come from this since
## this is the inverted of the vcov). There is no reason to expect
## this to be pd if there is multiple minima or if the best set of
## parameters wasnt found.
hess <- hessian(func = hessmin,
x = res$par,
sigma = sigma,
matrep = matrep,
model = model,
fun1 = fun)
## Gradients (indicate if the optimizer reach the minima). Here we
## used numDeriv::grad, which is a bit slower than
## lavaan::first.derivative.param.
gradd <- grad(func = gradmin,
x = res$par,
sigma = sigma,
matrep = matrep,
model = model,
fun = fun)
## Checking if any of the eigenvalues are positive.
## vcov and standard errors come from the inverse of this matrix. If
## it is not invertible, there is no way to calculate the confidence
## intervals.
if (any(eigen(hess)$values < 0L)) {
stop("Hessian matrix is not positive-definite. \
The optimizer did not converge to a unique \
solution. The model is probably empirically \
underidentified.")
}
## Any gradient not close to 0 might indicate that the system did
## not reach the minima.
#if (any(gradd > 1)) {
#warning("Some first order derivative are not \
#close to 0. The system probably did \
#not reach the minima. Results might \
#not be thurstworthy.")
#}
## More housecleaning and setting up
freeid <- model$id[model$free > 0L]
names(gradd) <- paste(model$lhs[freeid], model$op[freeid],
model$rhs[freeid], sep = "")
colnames(hess) <- rownames(hess) <- names(gradd)
## Calculating standard errors
std.errors <- se.sem(hess, n)
if (!std.errors$trust) warning("Could not compute standard errors.")
## Building the response
model$est <- model$ustart
model$est[model$free > 0L] <- res$par
model$se <- rep(0, length(model$id))
model$se[model$free > 0L] <- std.errors$values
## Calculating z-values
model$z <- model$est / model$se
model$z <- ifelse(model$z == Inf, NA, model$z)
## Significance of the parameter estimates
model$p <- round(do.call("pval", list(model$z)), 3)
out <- model[-1]
swap <- which(out$op == "~")
var1 <- out$lhs
var2 <- out$rhs
var1[swap] <- out$rhs[swap]
var2[swap] <- out$lhs[swap]
out$lhs <- var1
out$rhs <- var2
out$op[out$op == "~"] <- "-->"
out$op[out$op == "~~"] <- "<->"
out$op[out$op == "=~"] <- "-->"
out$free <- ifelse(out$free > 0, "free", "fixed")
out$sig <- rep("*", length(out$lhs))
out$sig[out$p >= 0.05] <- "ns"
out$sig[is.na(out$p)] <- NA
res$model <- data.frame(out, stringsAsFactors = FALSE)
colnames(res$model) <- c("var1", "op", "var2", "type", "start", "est", "se",
"z", "p", "sig")
## Calculating MLX2
res$chisquared <- n * res$objective
## Number of degrees of freedom
res$df <- model.df(model)$model.df
## Number of observations
res$n <- n
## Significance of the MLX2
res$pval <- 1 - pchisq(res$chisquared, res$df)
## Which function was minimised?
res$estimator <- fun
## Calculating the sigma hat with the best set of estimates
res$sigma.hat <- sem.matrices(model, matrep)
## Setting up stuff
res$sigma <- sigma
res$grad <- res$gradd
res$hessian <- hess
class(res) <- c("sem", "list")
invisible(res)
}
gradit <- function(pars, model, matrep = matrep, sigma, fun) {
grad(func = gradmin, x = pars, sigma = sigma,
model = model, matrep = matrep, fun = fun)
}
gradmin <- function(pars, model, matrep = matrep, sigma, fun) {
model$ustart[model$free > 0L] <- pars
sigma.hat <- sem.matrices(model, matrep)
do.call(fun, list(sigma, sigma.hat))
}
hessmin <- function(pars, model, matrep = matrep, sigma, fun1) {
model$ustart[model$free > 0L] <- pars
sigma.hat <- sem.matrices(model, matrep)
do.call(fun1, list(sigma, sigma.hat))
}
hessit <- function(pars, model, sigma, matrep, fun1) {
hessian(func = hessmin, x = pars, sigma = sigma,
model = model, matrep = matrep, fun1 = fun1)
}
se.sem <- function(hess, n) {
## From SAS:
## http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug
## _tcalis_sect072.htm
hess.inv <- try(chol2inv(chol(hess)))
if (class(hess.inv) == "try-error") {
return(list(trust = FALSE, values = rep(NA, NROW(hess))))
} else {
return(list(trust = TRUE, values =
sqrt(diag(2/(n-1)*hess.inv))))
}
}
## P-values of parameter estimates
pval <- function(z) {
2 * (1 - pnorm(abs(z)))
}
## Residuals
residuals.sem <- function(object, ...) {
round(object$sigma - object$sigma.hat, 4)
}
## Most of this function was taken from lavaan.
RMSEA <- function(fitted.model) {
rmsea <- sqrt(max(c((fitted.model$chisquared/fitted.model$n)/fitted.model$df - 1/fitted.model$n, 0)))
lower.lambda <- function(lambda) {
pchisq(fitted.model$chisquared, df=fitted.model$df, ncp=lambda) - 0.95
}
if (fitted.model$df < 1 || lower.lambda(0) < 0) {
ci.lower <- 0
} else {
lambda.l <- try(uniroot(f=lower.lambda, lower=0, upper=fitted.model$chisquared)$root)
if(inherits(lambda.l, "try-error")) { lambda.l <- NA }
ci.lower <- sqrt( lambda.l/(fitted.model$n*fitted.model$df) )
}
upper.lambda <- function(lambda) {
pchisq(fitted.model$chisquared, df=fitted.model$df, ncp=lambda) - 0.05
}
N.RMSEA <- max(fitted.model$n, 10000L)
if (fitted.model$df < 1 || upper.lambda(N.RMSEA) > 0) {
ci.upper <- 0
} else {
lambda.u <- try(uniroot(f=upper.lambda, lower=0, upper=N.RMSEA)$root)
if(inherits(lambda.u, "try-error")) { lambda.u <- NA }
ci.upper <- sqrt( lambda.u/(fitted.model$n*fitted.model$df) )
}
p <- 1 - pchisq(fitted.model$chisquared, df=fitted.model$df, ncp=(fitted.model$n*fitted.model$df*0.05^2))
list(RMSEA = rmsea, ci = c(lower = ci.lower, upper = ci.upper), p = p)
}
## Borrowed from lavaan
first.derivative.param <- function(x, model, matrep, sigma, fun) {
npar <- length(x)
h <- 10e-6
dx <- numeric(npar)
for(i in 1:npar) {
x.left <- x.left2 <- x.right <- x.right2 <- x
x.left[i] <- x[i] - h; x.left2[i] <- x[i] - 2*h
x.right[i] <- x[i] + h; x.right2[i] <- x[i] + 2*h
fx.left <- iterator(x.left, model = model, matrep = matrep,
sigma = sigma, fun = fun)
fx.left2 <- iterator(x.left2, model = model, matrep = matrep,
sigma = sigma, fun = fun)
fx.right <- iterator(x.right, model = model, matrep = matrep,
sigma = sigma, fun = fun)
fx.right2 <- iterator(x.right2, model = model, matrep = matrep,
sigma = sigma, fun = fun)
dx[i] <- (fx.left2 - 8*fx.left + 8*fx.right - fx.right2)/(12*h)
}
dx
}
summary.sem <- function(object, ...) {
space <- "\n\n"
if (object$convergence != 0) {
txt.con <- "The optimizer did not converge after %s iterations. Results are not trustworthy."
cat(space, sprintf(txt.con, object$iterations), space)
} else {
txt.con <- "The optimizer successfully converged after %s iterations"
cat(sprintf(txt.con, object$iterations), space)
cat("MODEL FIT", space)
cat(sprintf("Method %30s", toupper(object$estimator)), "\n")
cat(sprintf("Model chi-square %20.3f", object$chisquared), "\n")
cat(sprintf("P-value %29.3f", object$pval), "\n")
cat(sprintf("Degrees of freedom %18i", object$df), "\n")
cat(sprintf("Number of observations %14i", object$n), space)
cat("PARAMETER ESTIMATES\n\n")
#cat(sprintf("%30s %10s %10s %10s", "est", "std.error", "z",
#"p-value"), "\n")
#for (i in object$model$id) {
#cat(sprintf("%-5s %2s %5s %15.3f %10.3f %10.3f %10.3f", object$model$lhs[i],
#object$model$op[i], object$model$rhs[i], object$model$est[i],
#object$model$se[i], object$model$z[i], object$model$p[i]),
#"\n")
#ops <- unique(object$model$op[i:(i+1)])
#if (length(ops) > 1) {
#cat("\n")
#}
#}
print(object$model, digits = 3)
}
}
variables <- function(model) {
allvars <- unique(c(model$lhs, model$rhs))
latent <- unique(model$lhs[model$op == "=~"])
indicators <- unique(c(model$rhs[model$op == "=~" & !model$rhs %in%
latent]))
regressions <- lapply(model, `[`, model$op == "~" | (model$op == "=~"
& model$rhs %in%
latent))
exogenous <- unique(regressions$rhs[!regressions$rhs %in%
regressions$lhs])
endogenous <- unique(regressions$lhs)
latent.allvars <- allvars %in% latent
indicators.allvars <- allvars %in% indicators
exogenous.allvars <- allvars %in% exogenous
endogenous.allvars <- allvars %in% endogenous
#if (!all(exogenous.allvars) & ) exogenous.allvars <- latent.allvars
#latent.allvars
out <- cbind(latent.allvars, indicators.allvars, exogenous.allvars,
endogenous.allvars)
rownames(out) <- allvars
colnames(out) <- c("lat", "ind", "exo", "endo")
out
}
matrix.representation <- function(model) {
BE <- BETA(model)
PS <- PSI(model)
LY <- LAMBDA.Y(model)
TE <- THETA.EPSILON(model)
Inv <- II(model)
list(BE = BE, PS = PS, LY = LY, TE = TE, Inv = Inv)
}
sigma.hat <- function(matrep) {
matrep$LY %*% solve(matrep$Inv - matrep$BE) %*% matrep$PS %*%
solve(t(matrep$Inv - matrep$BE)) %*% t(matrep$LY) + matrep$TE
}
boot.lavaan <- function(fitted.model, n) {
sigma.hat <- fitted.values(fitted.model)$cov
nobs <- fitted.model@Sample@nobs[[1]]
mat <- mvrnorm(Sigma = sigma.hat, n = nobs, mu = rep(0,
NCOL(sigma.hat)), empirical = TRUE)
est.coefs <- coef(fitted.model)
X2 <- numeric(n)
coefs <- matrix(0, nrow = length(est.coefs), ncol = n)
for (i in 1:n) {
success <- FALSE
while (!success) {
Z <- mat[sample(1:NROW(mat), replace = TRUE), ]
new.fit <- update(fitted.model, sample.cov =
cov(Z), warn = FALSE)
if (sum(new.fit@Fit@se > 0) == length(est.coefs)) {
success <- TRUE
X2[i] <- new.fit@Fit@test[[1]]$stat
coefs[, i] <- coef(new.fit)
}
}
}
p.X2 <- sum(X2 > fitted.model@Fit@test[[1]]$stat)/n
se <- apply(coefs, 1, sd)
z <- est.coefs/se
p.est <- pval(z)
invisible(list(X2 = X2, p.X2 = p.X2, est = est.coefs, se = se, z = z, p.est =
p.est, coefs = coefs))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.