# R/cumres.lvmfit.R In gof: Model Diagnostics Based on Cumulative Residuals

#### Documented in cumres.lvmfit

cumresMaxSet <- function(m,var,...) {
A <- t(lava::index(m)$A) Afix <- A; Afix[t(lava::index(m)$M0)==1] <- 0
A[A!=0] <- 1
P <- lava::index(m)$P k <- nrow(A) I <- diag(k) B <- rbind(I,solve(I-A)) VV <- B%*%P%*%t(B) u.var <- lava::index(m)$vars
V0 <- VV[seq(length(lava::vars(m))),seq(length(lava::vars(m)))+length(lava::vars(m))][,lava::endogenous(m)]
rownames(V0) <- lava::vars(m)
names(which(V0[var,]==0))
}

##' Cumulative residual processes for structural equation models
##'
##' Calculates GoF statistics based on cumulative residual processes for
##' structural equation models fitted with the \code{lava} package.
##'
##' With \code{y} and \code{x} given as functions the user can decide which
##' variables to use in the prediction of the outcome and predictor (use the
##' \code{predict} method as below).
##'
##' @param model \code{lvm} object
##' @param y A formula specifying the association to be checked. Alternatively
##' the outcome specified as a function or a string with the name of the outcome
##' in the model.
##' @param x Predictor. A function, vector or character
##' @param full If FALSE the prediction, Pr, of the variable that are ordered
##' after is only calculated based on the conditional distribution given
##' covariates. If TRUE the conditional expectation is based on the largest set
##' of covariates and endogenous variables such that the residual and Pr are
##' uncorrelated.
##' @param data data.frame (default is the model.frame of the model)
##' @param p Optional parameter vector
##' @param R Number of processes to simulate
##' @param b Moving average parameter
##' @param plots Number of processes to save for use with the plot method
##' @param \dots Additional arguments parsed on to lower-level functions
##' @method cumres lvmfit
##' @export
##' @return Returns a \code{cumres} object with associated
##' \code{plot},\code{print},\code{coef} methods
##' @author Klaus K. Holst
##' @references B.N. Sanchez and E. A. Houseman and L. M. Ryan (2009)
##' \emph{Residual-Based Diagnostics for Structural Equation Models}. Biometrics
##' Volume 65 (1), pp 104-115.
##' @keywords models regression
##' @examples
##'
##' \donttest{
##' library(lava)
##' m <- lvm(list(c(y1,y2,y3)~eta,eta~x)); latent(m) <- ~eta
##' ## simulate some data with non-linear covariate effect
##' functional(m,eta~x) <- function(x) 0.3*x^2
##' d <- sim(m,100)
##'
##' e <- estimate(m,d)
##' ## Checking the functional form of eta on x
##' g <- cumres(e,eta~x,R=1000)
##' plot(g)
##'
##' x <- function(p) predict(e,x=~y2+y3,p=p)[,"eta"]
##' ## Checking the functional form of y1 on eta
##' cumres(e,y1~eta,R=1000)
##' g <- cumres(e,"y1",x=x,R=1000)
##' plot(g)
##'
##' }
cumres.lvmfit <- function(model,y,x,full=FALSE,
data=model.frame(model),p,
R=1000, b=0, plots=min(R,50),
...) {
if (!requireNamespace("numDeriv")) stop("package 'numDeriv' not available")

cl <- match.call()

if (class(y)[1]=="formula") {
y <- lava::getoutcome(y)
x <- attributes(y)$x y <- lava::decomp.specials(y) } if (is.list(y) && class(y[[1]])[1]=="formula") { yy <- y y <- list(); x <- list() res <- c() cl[[1]] <- as.name("cumres") for (i in 1:length(yy)) { resp <- lava::getoutcome(yy[[i]]) y <- lava::decomp.specials(resp) x <- attributes(resp)$x
cl$y <- y; cl$x <- x
res <- c(res, list(eval.parent(cl)))
}
return(res)
}

cl$y <- y; cl$x <- x
cl[[1]] <- as.name("cumres")
if ((is.character(y) | is.list(y)) & length(y)>1) {
res <- c()
iy <- 0;
for (i in y) {
iy <- iy+1
cl$y <- i yname <- ifelse(is.character(i),i,paste("y",iy,sep="")) if (is.character(x) | is.list(x)) { ix <- 0 for (j in x) { ix <- ix+1 cl$x <- j
xname <- ifelse(is.character(j),j,paste("x",ix,sep=""))
newres <- list(eval.parent(cl))
names(newres) <- paste(yname,xname,sep="<-")
res <- c(res,newres)
}
} else {
xname <- ifelse(is.character(j),j,"x1")
newres <- list(eval.parent(cl))
names(newres) <- paste(yname,xname,sep="<-")
res <- c(res,newres)
}
}
return(res)
}

if ((is.character(x) | is.list(x)) & length(x)>1) {
res <- c()
ix <- 0
yname <- ifelse(is.character(y),y,"y1")
for (j in x) {
ix <- ix+1
cl$x <- j xname <- ifelse(is.character(j),j,paste("x",ix,sep="")) newres <- list(eval.parent(cl)) names(newres) <- paste(yname,xname,sep="<-") res <- c(res,newres) } return(res) } if (missing(p)) p <- lava::pars(model) x0 <- x if (is.function(x)) x0 <- x(p) if (is.character(x)) { if (full) { predictby <- cumresMaxSet(model,y,...) x0 <- predict(model,predictby)[,x] } else { if (x %in% lava::latent(model)) { x0 <- predict(model,~1)[,x] } else { x0 <- data[,x] } } } n <- length(x0) if (is.function(y)) { r <- y(p) grad <- attributes(r)$grad
grad <- numDeriv::jacobian(y,p,method=lava::lava.options()$Dmethod) } else { myres <- function(p) { rr <- predict(model,vars(model),residual=TRUE,p=p) if (is.matrix(rr)) return(rr[,y]) else return(rr) } r <- myres(p) grad <- numDeriv::jacobian(myres,p,method=lava::lava.options()$Dmethod)
}

ord <- order(x0)
variable <- data.frame(x0); names(variable) <- ifelse(is.character(x),x,"x")
res <- cumres.default(NULL, variable=variable, r=r, dr=grad, ic=iid(model),
R=R, plots=plots, b=b, ...)
res\$type <- "sem"
return(res)
}

################################################################################


## Try the gof package in your browser

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

gof documentation built on March 16, 2020, 1:07 a.m.