R/lisrel.R

In kkholst/lava: Latent Variable Models

lisrel <- function(model,p,X=NULL,muX=NULL,varX=NULL,...) {
    mom <- moments(model,p)
    A <- t(index(model)$M)

    eta.idx <- match(latent(model),vars(model))
    exo.idx <- match(exogenous(model),vars(model))
    y <- setdiff(manifest(model), exogenous(model))
    y.idx <- match(y, vars(model))

    ##  Jy <- Jx <- Jeta <- I <- diag(length(vars(model)))
    ##  if (length(eta.idx)>0)
    ##    J[eta.idx,eta.idx] <- 0; J <- J[-eta.idx,]
    ##  Jeta[obs.idx,obs.idx] <- 0; Jeta <- J[-obs.idx,]

    A <- t(mom$A)
    Lambda <- A[y.idx,eta.idx,drop=FALSE]
    K <- A[y.idx,exo.idx,drop=FALSE]
    B <- A[eta.idx,eta.idx,drop=FALSE]
    I <- diag(nrow=nrow(B))
    Gamma <- A[eta.idx,exo.idx,drop=FALSE]
    V <- mom$P
    Psi <- V[eta.idx,eta.idx] ## Residual variance

    Theta <- V[y.idx,y.idx] ## -
    IBi <- if (ncol(I)>0) solve(I-B) else I
    LIBi <- Lambda%*%IBi
    Phi <- LIBi%*%Gamma + K

    Veta.x <- IBi%*%Psi%*%IBi  ## Variance of eta given x
    COVetay.x <- Veta.x%*%t(Lambda) ## Covariance of eta,y given x
    ##  Vy.x <- Lambda%*%COVetay.x + Theta ## Omega
    Vy.x <- LIBi%*%Psi%*%t(LIBi) + Theta

    if (!is.null(X)) {
        Ey.x <- t(apply(as.matrix(X)%*% t(LIBi%*%Gamma + K),1,function(x) x + mom$v[y.idx]))
    } else Ey.x <- NULL

    CV <- COVetay.x%*%Vy.x
    ##  Sigma <- Vy.x + Phi%*%varX%*%t(Phi)

    return(list(mu=mom$v,
                Lambda=Lambda, K=K, B=B, I=I, Gamma=Gamma, Psi=Psi, Theta=Theta, IBi=IBi, LIBi=LIBi, Phi=Phi,
                Vy.x=Vy.x, Veta.x=Veta.x, COVetay.x=COVetay.x, CV=CV, Ey.x=Ey.x))
}