#' Model predicted residual outliers
#'
#' Compute model predicted residuals for each variable using regression
#' estimated factor scores.
#'
#' @aliases obs.resid
#' @param data matrix or data.frame
#' @param model if a single numeric number declares number of factors to extract in
#' exploratory factor analysis. If \code{class(model)} is a sem (semmod), or lavaan (character),
#' then a confirmatory approach is performed instead
#' @author Phil Chalmers \email{rphilip.chalmers@@gmail.com}
#' @seealso
#' \code{\link{gCD}}, \code{\link{LD}}, \code{\link{robustMD}}
#' @references
#'
#' Chalmers, R. P. & Flora, D. B. (2015). faoutlier: An R Package for Detecting
#' Influential Cases in Exploratory and Confirmatory Factor Analysis.
#' \emph{Applied Psychological Measurement, 39}, 573-574. \doi{10.1177/0146621615597894}
#'
#' Flora, D. B., LaBrish, C. & Chalmers, R. P. (2012). Old and new ideas for data
#' screening and assumption testing for exploratory and confirmatory factor analysis.
#' \emph{Frontiers in Psychology, 3}, 1-21. \doi{10.3389/fpsyg.2012.00055}
#' @keywords covariance
#' @export obs.resid
#' @examples
#'
#' \dontrun{
#' data(holzinger)
#' data(holzinger.outlier)
#'
#' #Exploratory
#' nfact <- 3
#' (ORresult <- obs.resid(holzinger, nfact))
#' (ORresult.outlier <- obs.resid(holzinger.outlier, nfact))
#' plot(ORresult)
#' plot(ORresult.outlier)
#'
#' #-------------------------------------------------------------------
#' #Confirmatory with sem
#' model <- sem::specifyModel()
#' F1 -> Remndrs, lam11
#' F1 -> SntComp, lam21
#' F1 -> WrdMean, lam31
#' F2 -> MissNum, lam41
#' F2 -> MxdArit, lam52
#' F2 -> OddWrds, lam62
#' F3 -> Boots, lam73
#' F3 -> Gloves, lam83
#' F3 -> Hatchts, lam93
#' F1 <-> F1, NA, 1
#' F2 <-> F2, NA, 1
#' F3 <-> F3, NA, 1
#'
#' (ORresult <- obs.resid(holzinger, model))
#' (ORresult.outlier <- obs.resid(holzinger.outlier, model))
#' plot(ORresult)
#' plot(ORresult.outlier)
#'
#' #-------------------------------------------------------------------
#' #Confirmatory with lavaan
#' model <- 'F1 =~ Remndrs + SntComp + WrdMean
#' F2 =~ MissNum + MxdArit + OddWrds
#' F3 =~ Boots + Gloves + Hatchts'
#'
#' (obs.resid2 <- obs.resid(holzinger, model, orthogonal=TRUE))
#' (obs.resid2.outlier <- obs.resid(holzinger.outlier, model, orthogonal=TRUE))
#' plot(obs.resid2)
#' plot(obs.resid2.outlier)
#'
#' }
obs.resid <- function(data, model, ...)
{
ret <- list()
rownames(data) <- 1:nrow(data)
N <- nrow(data)
if(any(is.na(data)))
stop('All routines require complete datasets (no NA\'s)')
if(!all(sapply(data, function(x) is.numeric(x) || is.integer(x))))
stop('data must conist of numeric elements')
if(is.numeric(model)){
R <- cor(data)
mod <- factanal(data, model, rotation='none', scores = 'regression', ...)
scores <- mod$scores
ret$fascores <- scores
Lambda <- unclass(mod$loadings)
Theta <- diag(mod$uniquenesses)
e <- data - scores %*% t(Lambda)
VAR <- Theta %*% solve(R) %*% Theta
eji <- t(solve(diag(sqrt(diag(VAR)))) %*% t(scale(e, scale = FALSE)))
colnames(eji) <- colnames(e) <- colnames(data)
ret$res <- e
ret$std_res <- eji
} else if(class(model) == "semmod"){
C <- cov(data)
vnames <- colnames(C)
mod <- sem::sem(model, C, N, ...)
scores <- sem::fscores(mod, data)
ret$fascores <- scores
lnames <- setdiff(colnames(mod$P), vnames)
Lambda <- mod$A[1:length(vnames), (length(vnames)+1):ncol(mod$A)]
Theta <- mod$P[1:length(vnames),1:length(vnames)]
e <- data - scores %*% t(Lambda)
VAR <- Theta %*% solve(C) %*% Theta
eji <- t(solve(diag(sqrt(diag(VAR)))) %*% t(scale(e, scale = FALSE)))
colnames(eji) <- colnames(e) <- colnames(data)
ret$res <- e
ret$std_res <- eji
} else if(class(model) == "character"){
dots <- list(...)
if(!is.null(dots$ordered))
stop('observed residuals only defined for continuous data')
mod <- lavaan::sem(model, data=data, ...)
C <- mod@SampleStats@cov[[1L]]
scores <- lavaan::predict(mod)
ret$fascores <- scores
Lambda <- mod@Model@GLIST$lambda
Theta <- mod@Model@GLIST$theta
Psi <- mod@Model@GLIST$psi
dat <- data[,mod@Data@ov.names[[1L]]]
e <- dat - scores %*% t(Lambda)
VAR <- Theta %*% solve(C) %*% Theta
eji <- t(solve(diag(sqrt(diag(VAR)))) %*% t(scale(e, scale = FALSE)))
colnames(eji) <- colnames(e) <- colnames(dat)
ret$res <- e
ret$std_res <- eji
} else {
stop('model class not supported')
}
ret$id <- rownames(data)
class(ret) <- 'obs.resid'
ret
}
#' @rdname obs.resid
#' @param x an object of class \code{obs.resid}
#' @param restype type of residual used, either \code{'obs'} for observation value
#' (inner product), \code{'res'} or \code{'std_res'} for unstandardized and standardized
#' for each variable, respectively
#' @param ... additional parameters to be passed
#' @export
print.obs.resid <- function(x, restype = 'obs', ...)
{
if(restype == 'res') return(print(x$res))
if(restype == 'std_res') return(print(x$std_res))
stat <- c()
for(i in 1:length(x$id))
stat[i] <- x$std_res[i, ] %*% x$std_res[i, ]
if(restype == 'obs') return(stat)
}
#' @rdname obs.resid
#' @param y a \code{NULL} value ignored by the plotting function
#' @param main the main title of the plot
#' @param type type of plot to use, default displays points and lines
#' @export
plot.obs.resid <- function(x, y = NULL, main = 'Observed Residuals',
type = c('p','h'), restype = 'obs', ...)
{
ylab <- switch(restype,
obs = 'Observed residuals',
res = 'Observed variable residuals',
std_res = 'Standardized variable residuals')
ID <- as.numeric(x$id)
stat <- c()
for(i in 1:length(x$id))
stat[i] <- x$std_res[i, ] %*% x$std_res[i, ]
if(restype == 'obs'){
dat <- data.frame(stat=stat,ID=ID)
ret <- lattice::xyplot(stat~ID, dat, type = type, main = main, ylab = ylab, ...)
}
if(restype == 'res' || restype == 'std_res'){
if(restype == 'res') dat <- data.frame(ID=ID,x$res)
else dat <- data.frame(ID=ID,x$std_res)
rownames(dat) <- ID
dat2 <- reshape(dat, v.names='res', direction = 'long', varying=2:ncol(dat),
times = colnames(dat)[-1], timevar='variable', sep='')
ret <- lattice::xyplot(res~ID|as.factor(variable), dat2, type = type, main = main,
ylab = ylab, ...)
}
return(ret)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.