#' Non-linear structural equation models
#'
#' Fits non-linear SEM via 1. order Laplace Approximation
#'
#' The \code{model} arguments must be a named list: list( measure0, measure1,
#' measure2, pred0, pred1, pred2, model)
#'
#' where 'model' currently can be either 'nsem2' (2-factor model) or 'nsem3'
#' (3-factor model). Each of the remaining elements can be either character
#' vectors or formulas.
#'
#' Alternatively, a list of lists can be used as the \code{model} argument in
#' which case a list of data.frames also much be given as \code{data} argument.
#' A multigroup model is then fitted.
#'
#' @aliases nsem coef.lava.nlin logLik.lava.nlin score.lava.nlin vcov.lava.nlin
#' print.lava.nlin
#' @param model model list
#' @param data data.frame
#' @param laplace.control Options to the Laplace Approximation algorithm
#' @param control Options parsed on to optimizer (nlminb)
#' @param method Optimization method
#' @param vcov Boolean indicating whether calculation of asymptotic standard
#' @param p (Optional) evaluate log-likelihood in p
#' @param fast if TRUE do not spend time on returning lava object/graph
#' errors should be disabled
#' @param \dots Additional parameters parsed on to lower level functions
#' @return \code{lava.nlin} object (available methods: \code{coef},
#' \code{logLik}, \code{score}, \code{vcov}, \code{print}, ...)
#' @author Klaus K. Holst
#' @keywords models regression
#' @examples
#'
#' \dontrun{
#' model1 <- list(measure1=~parc+supfc+antcin+postcin,measure2=~thS+putS+midbrainS,pred0=~age+bmi,model="nsem3")
#' b <- nsem(model1,data=dtrp)
#' }
#'
#' @export nsem
nsem <- function(model,
data,
laplace.control=list(),
control=list(trace=1),
method=c("nlminb","BFGS","ucminf"),
vcov=TRUE,
p,fast=FALSE,
...) {
procmod <- function(M,data,...) {
if (class(M$measure0)[1]=="formula") M$measure0 <- all.vars(M$measure0)
if (class(M$measure1)[1]=="formula") M$measure1 <- all.vars(M$measure1)
if (class(M$measure2)[1]=="formula") M$measure2 <- all.vars(M$measure2)
if (class(M$pred0)[1]=="formula") M$pred0 <- all.vars(M$pred0)
if (class(M$pred1)[1]=="formula") M$pred1 <- all.vars(M$pred1)
if (class(M$pred2)[1]=="formula") M$pred2 <- all.vars(M$pred2)
mydata <- data[,c(M$measure0,M$measure1,M$measure2)]
if (!is.null(M$pred0)) mydata <- cbind(mydata,data[,M$pred0])
if (!is.null(M$pred1)) mydata <- cbind(mydata,data[,M$pred1])
if (!is.null(M$pred2)) mydata <- cbind(mydata,data[,M$pred2])
mydata <- na.omit(as.matrix(mydata))
nn <- c(M$measure0,M$measure1,M$measure2)
if (length(M$measure0)>0) nn <- c(nn,paste(M$measure0,"eta0",sep="~"))
if (length(M$measure1)>1) nn <- c(nn,paste(M$measure1[-1],"eta1",sep="~"))
if (length(M$measure2)>1) nn <- c(nn,paste(M$measure2[-1],"eta2",sep="~"))
if (length(M$pred0)>0) nn <- c(nn,paste("eta0",M$pred0,sep="~"))
if (length(M$pred1)>0) nn <- c(nn,paste("eta1",M$pred1,sep="~"))
if (length(M$pred2)>0) nn <- c(nn,paste("eta2",M$pred2,sep="~"))
nn <- c(nn,
switch(M$model,
nsem3b=c("eta2~eta0","eta2~eta0^2",
"eta1~eta0","eta1~eta0^2",
"eta0,eta0",
"eta1,eta1","eta2,eta2"),
nsem3=c("eta2~eta0","eta2~eta0^2",
"eta0,eta0","eta1,eta1","eta2,eta2"),
nsem2=c("eta2~eta1","eta2~eta1^2",
"eta1,eta1","eta2,eta2")))
omit <- switch(M$model,
nsem3b="eta0,eta0",
NULL)
nlatent <- switch(M$model,
nsem3b=3,
nsem3=3,
nsem2=2)
if (length(M$measure0)>0) nn <- c(nn,paste(M$measure0,M$measure0,sep=","))
if (length(M$measure1)>0) nn <- c(nn,paste(M$measure1,M$measure1,sep=","))
if (length(M$measure2)>0) nn <- c(nn,paste(M$measure2,M$measure2,sep=","))
if (M$model=="nsem2") {
nn[1] <- "eta1"
nn[1+length(M$measure1)] <- "eta2"
}
mm <- list(nlatent=nlatent, nvar0=length(M$measure0), nvar1=length(M$measure1), nvar2=length(M$measure2), npred0=length(M$pred0), npred1=length(M$pred1), npred2=length(M$pred2))
## npar <- c()
## theta0 <- rep(0,with(mm,2*nvar0 + 2*(nvar1+nvar2-1) + npred0+npred1+npred2+
## 2+nvar0+nvar1+nvar2+nlatent))
res <- list(data=mydata,names=nn,modelpar=mm,model=M$model,omit=omit,
measure0=M$measure0,measure1=M$measure1,measure2=M$measure2,
pred0=M$pred0,pred1=M$pred1,pred2=M$pred2)
return(res)
}
models <- c()
if (is.list(model[[1]])) { ## multigroup
for (i in 1:length(model)) {
models <- c(models, list(procmod(model[[i]],data[[i]])))
}
} else {
models <- list(procmod(model,data))
}
allnames <- unique(unlist(lapply(models,function(x) x$names)))
for (i in 1:length(models)) {
models[[i]]$idx <- match(models[[i]]$names,allnames)
}
allomit <- unique(unlist(lapply(models,function(x) x$omit)))
pidx <- which(!(allnames%in%allomit))
laplace.control0 <- list(lambda=0.5,niter=50,Dtol=1e-4,nq=0,eb0=FALSE)
laplace.control0[names(laplace.control)] <- laplace.control
laplace.control <- laplace.control0
eta0 <- NULL
if (is.matrix(laplace.control$eb0)) eta0 <- laplace.control$eb0
f <- function(p,eta.update=laplace.control$eb0,indiv=FALSE,...) { ## -log-likelihood
p0 <- rep(0,length(allnames))
p0[pidx] <- p
many <- (indiv||eta.update)&&(length(models)==1)
val <- with(models[[1]],Lapl(data,p0[idx],modelpar,model=model,control=laplace.control,indiv=many,eb0=eta0,...))
if (length(models)==1) {
if (eta.update) {
eta0 <<- val[,-1,drop=FALSE] ## Write eta0 in outer scope
val <- val[,1]
}
## val <- with(models[[1]],Lapl(data,p0[idx],modelpar,model=model,control=laplace.control,indiv=many,eb0=eta0,...))[,1]
return(-sum(val))
}
eta0 <- NULL
for (i in 2:length(models))
val <- -val - with(models[[i]],Lapl(data,p0[idx],modelpar,model=model,control=laplace.control,indiv=FALSE,...))
return(val)
}
if (!missing(p)) {
logL <- f(p,...)
return(logL)
}
theta0 <- rep(-0.2,length(setdiff(allnames,allomit))) ## Starting values
if (!is.null(control$start)) {
if (length(control$start)==length(theta0)) {
theta0 <- control$start
} else {
theta0[match(names(control$start),allnames)] <-
control$start[which(control$start%in%allnames)]
}
}
control$start <- NULL
control0 <- list(trace=1,kkt=FALSE)
control0[names(control)] <- control
control <- control0
## res.Laplace <- tryCatch(nlminb(theta0,f,control=control),error=function(e) NULL)
## Nelder-Mead',
## 'BFGS', 'CG', 'L-BFGS-B', 'nlm', 'nlminb', 'spg', 'ucminf',
## 'newuoa', 'bobyqa', 'nmkb', 'hjkb', 'Rcgmin', or 'Rvmmin'.
op <- tryCatch(optimx::optimx(theta0,f,method=method,control=control),error=function(e) NULL)
fop <- apply(coef(op),1,f)
iop <- which.min(fop)
res.Laplace <- list(par=as.vector(coef(op)[iop,]), value=op[iop,"value"], opt=op)
if (is.null(res.Laplace)) stop("Optimization error")
if (vcov) {
laplace.control$eb0 <- FALSE
S0 <- numDeriv::grad(f,res.Laplace$par)
res.Laplace$grad <- S0
H0 <- numDeriv::hessian(f,res.Laplace$par)
vcov0 <- tryCatch(solve(H0),error=function(e) matrix(NA,ncol=ncol(H0),nrow=nrow(H0)))
} else {
S0 <- NULL
vcov0 <- matrix(NA,nrow=length(theta0),ncol=length(theta0))
}
colnames(vcov0) <- rownames(vcov0) <- setdiff(allnames,allomit)
names(res.Laplace$par) <- setdiff(allnames,allomit)
mycoef <- cbind(res.Laplace$par,diag(vcov0)^0.5)
mycoef <- cbind(mycoef,mycoef[,1]/mycoef[,2],2*(1-pnorm(abs(mycoef[,1]/mycoef[,2]))))
rownames(mycoef) <- setdiff(allnames,allomit); colnames(mycoef) <- c("Estimate","Std.Err","Z-value","P(>|z|)")
res <- list(par=res.Laplace$par, coefmat=mycoef, vcov=vcov0, opt=res.Laplace, score=S0, data=data, models=models)
class(res) <- "lava.nlin"
if (model$model=="nsem2" && length(models)==1 && !fast) {
M <- models[[1]]
m <- lvm()
m <- regression(m,from="eta1",to=M$measure1)
m <- regression(m,from="eta2",to=M$measure2)
m <- regression(m,from=c("eta1","eta1^2"),to="eta2")
if (length(M$pred1)>0)
m <- regression(m,from=M$pred1,to="eta1")
if (length(M$pred2)>0)
m <- regression(m,from=M$pred2,to="eta2")
intercept(m,c("eta1^2",M$measure1[1],M$measure2[1])) <- 0
covariance(m,"eta1^2") <- 0
regression(m,from="eta1",to=M$measure1[1]) <- 1
regression(m,from="eta2",to=M$measure2[1]) <- 1
latent(m) <- c("eta1","eta2","eta1^2")
transform(m, y="eta1^2",x="eta1") <- function(x) x^2
labels(m) <- c("eta1"=expression(eta[1]),"eta2"=expression(eta[2]),"eta1^2"=expression(eta[1]^2))
pn <- names(pars(res))
## pn <- gsub("\\^","",pn)
pm <- coef(m)
m$order <- match(pn,pm)
ee <- structure(c(res, list(model=m,data=data)),class=c("lava.nlin","lvmfit"))
ee$coef <- ee$coefmat[match(pm,pn),]
p <- length(M$measure1)+length(M$measure2)
rownames(ee$coef)[seq(p+1,length(pm))] <- paste0("p",seq(length(pm)-p))
rownames(ee$coef)[seq(p)] <- paste0("m",seq(p))
return(ee)
}
return(res)
}
##' @export
plot.lava.nlin <- function(x,diag=FALSE,f,...) {
if (inherits(x,"lvmfit") & missing(f)) {
NextMethod("plot",diag=diag,...)
} else {
return(plot(estimate(NULL,coef=pars(x),vcov=vcov(x)),f=f,...))
}
}
##' @export
print.lava.nlin <- function(x,...) {
if (inherits(x,"lvmfit")) {
return(NextMethod("print",...))
}
printCoefmat(x$coefmat)
return(invisible(x))
}
##' @export
vcov.lava.nlin <- function(object,...) object$vcov
##' @export
score.lava.nlin <- function(x,...) x$score
##' @export
pars.lava.nlin <- function(x,...) x$par
##' @export
logLik.lava.nlin <- function(object,...) {
loglik <- -object$opt[["value"]]
if (is.null(attr(loglik, "nall")))
attr(loglik, "nall") <- nrow(object$data)
if (is.null(attr(loglik, "nobs")))
attr(loglik, "nobs") <- nrow(object$data) - length(pars(object))
if (is.null(attr(loglik, "df")))
attr(loglik, "df") <- length(pars(object))
class(loglik) <- "logLik"
return(loglik)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.