#' Remove random effects from all.vars
#'
#' @keywords internal
#'
all_vars_merMod <- function(formula.) {
if(!any(class(formula.) %in% c("formula", "formula.cerror"))) formula. <- formula(formula.)
if(inherits(formula., "formula.cerror"))
gsub(" " , "", unlist(strsplit(formula., "~~"))) else {
n <- rownames(attr(terms(formula.), "factors"))
if(any(grepl("\\|", n)))
all.vars(lme4::nobars(formula.)) else
all.vars(formula.)
}
}
#' Get vector of untransformed variables
#'
#' @keywords internal
#'
all_vars_notrans <- function(formula.) {
if(!all(class(formula.) %in% c("formula", "formula.cerror"))) formula. <- formula(formula.)
if(inherits(formula., "formula")) {
if(any(grepl("\\|", formula.))) formula. <- lme4::nobars(formula.)
formula. <- all_vars_trans(formula.)
if(any(grepl(":", formula.))) {
idx <- which(grepl(":", formula.))
for(i in idx) formula.[i] <- paste(sapply(strsplit(formula.[i], ":"), stripTransformations), collapse = ":")
for(j in (1:length(formula.))[-idx]) formula.[j] <- stripTransformations(formula.[j])
} else {
formula. <- sapply(formula., stripTransformations)
}
} else formula. <- unlist(strsplit(formula., " ~~ "))
ret <- gsub("(,.*)", "", formula.)
return(ret)
}
#' Get vector of transformed variables
#'
#' @keywords internal
#'
all_vars_trans <- function(formula., smoothed = FALSE) {
if(!all(class(formula.) %in% c("formula", "formula.cerror"))) formula. <- formula(formula.)
if(inherits(formula., "formula")) {
if(formula.[[3]] == 1) ret <- deparse(formula.[[2]]) else {
if(any(grepl("\\|", formula.))) formula. <- lme4::nobars(formula.)
ret <- c(rownames(attr(terms(formula.), "factors"))[1], labels(terms(formula.)))
if(smoothed == FALSE) ret <- gsub("(.*)\\,.*", "\\1", gsub("s\\((.*)\\).*", "\\1", ret))
# else {
# ret <- gsub("(s\\(.*),.*", "\\1", ret)
#
# if(any(grepl("s\\(", ret))) ret <- sapply(ret, function(x)
# ifelse(grepl("s\\(", x) & !grepl("\\)", x), paste0(x, ")"), x))
# }
# ret <- gsub("(,.*)", "", ret)
}
return(ret)
} else unlist(strsplit(formula., " ~~ "))
}
#' Captures output table
#'
#' @keywords internal
#'
captureTable <- function(g, row.names = FALSE) {
g1 <- capture.output(print(g, row.names = row.names))
if(all(g1 == "data frame with 0 columns and 0 rows"))
g1 <- "No independence claims present. Tests of directed separation not possible."
g1 <- paste0(g1, "\n")
return(g1)
}
#' Bind data.frames of differing dimensions
#'
#' From: https://stackoverflow.com/a/31678079
#'
#' @param ... data.frames to be bound, separated by commas
#'
#' @keywords internal
#'
cbind_fill <- function(...) {
nm <- list(...)
dfdetect <- grepl("data.frame|matrix", unlist(lapply(nm, function(cl) paste(class(cl), collapse = " ") )))
vec <- data.frame(nm[!dfdetect])
n <- max(sapply(nm[dfdetect], nrow))
vec <- data.frame(lapply(vec, function(x) rep(x, n)))
if (nrow(vec) > 0) nm <- c(nm[dfdetect], list(vec))
nm <- lapply(nm, as.data.frame)
do.call(cbind, lapply(nm, function (df1)
rbind(df1, as.data.frame(matrix(NA, ncol = ncol(df1), nrow = n-nrow(df1), dimnames = list(NULL, names(df1))))) ))
}
#' Transform variables based on model formula and store in new data frame
#'
#' @keywords internal
#'
dataTrans <- function(formula., data) {
notrans <- all_vars_merMod(formula.)
if(inherits(formula., "formula.cerror")) notrans <- gsub(".*\\((.*)\\)", "\\1", notrans)
trans <- all_vars_trans(formula.)
trans <- unlist(strsplit(trans, "\\:"))
trans <- trans[!duplicated(trans)]
if(any(grepl("scale\\(.*\\)", trans))) {
#
# trans[which(grepl("scale(.*)", trans))] <- notrans[which(grepl("scale(.*)", trans))]
#
warning("`scale` applied directly to variable. Use argument `standardize = TRUE` instead.", call. = FALSE)
}
if(any(!notrans %in% trans)) {
for(k in 1:length(notrans)) {
if(is.factor(data[, notrans[k]])) next else
if(grepl("scale(.*)", trans[k])) data[, notrans[k]] <- scale(data[, notrans[k]]) else
data[, notrans[k]] <-
sapply(data[, notrans[k]], function(x) eval(parse(text = gsub(notrans[k], x, trans[k]))))
}
}
colnames(data) <- notrans
return(data)
}
#' Get ANOVA results from `merMod`
#'
#' @keywords internal
#'
getAnova <- function(model, test.statistic = "F", test.type = "III") {
if(inherits(model, "glmmTMB")) test.statistic = "Chisq"
krp <- as.data.frame(car::Anova(model, test.statistic = test.statistic, type = test.type))
ct <- summary(model)$coefficients
colnames(ct)[2] <- "Std.Error"
ret <- do.call(rbind, lapply(1:nrow(krp), function(i) {
if(rownames(krp)[i] %in% rownames(ct)) {
cbind.data.frame(
ct[i, 1:2, drop = FALSE],
DF = krp[i, 3],
Crit.Value = krp[i, 1],
P = krp[i, ncol(krp)],
row.names = NULL
)
} else {
data.frame(
Estimate = NA,
Std.Error = NA,
DF = krp[i, 3],
Crit.Value = krp[i, 1],
P = krp[i, ncol(krp)]
)
}
} ) )
# ret <- cbind.data.frame(
# ct[, 1:2],
# DF = krp[, 3],
# Crit.Value = krp[, 1],
# P = krp[, ncol(krp)]
# )
names(ret)[ncol(ret)] <- "Pr(>|t|)"
rownames(ret) <- rownames(krp)
return(ret)
}
#' Get random effects from lme
#'
#' @keywords internal
#'
findbars.lme <- function(model) {
rand <- model$call$random
sapply(rand, function(i) {
i = gsub(".*\\|(.*)", "\\1", as.character(i)[2])
strsplit(gsub(" " , "", i), "\\/")[[1]]
} )
}
#' Get data from model list
#'
#' @keywords internal
#'
GetData <- function(modelList) {
if(!all(class(modelList) %in% c("psem", "list"))) modelList <- list(modelList)
modelList <- removeData(modelList, formulas = 1)
data.list <- lapply(modelList, GetSingleData)
data.list <- data.list[!sapply(data.list, is.null)]
data.list <- unname(data.list)
if(all(sapply(data.list, class) == "comparative.data"))
data <- data.list[[1]] else
data <- do.call(cbind_fill, data.list)
data <- data[, !duplicated(colnames(data), fromLast = TRUE)]
# colnames(data) <- gsub(".*\\((.*)\\).*", "\\1", colnames(data))
data <- as.data.frame(data)
rownames(data) <- 1:nrow(data)
return(data)
}
#' Get data from one model
#'
#' @keywords internal
#'
GetSingleData <- function(model) {
dat <- data.frame()
switch(class(model)[1],
"phylolm" = {
stop("Please provide `data =` argument to `psem`.", call. = FALSE)
},
"phyloglm" = {
stop("Please provide `data =` argument to `psem`.", call. = FALSE)
},
"lm" ={
dat <- eval(getCall(model)$data, environment(formula(model)))
},
"negbin" = {
dat <- eval(getCall(model)$data, environment(formula(model)))
},
"Sarlm" = {
dat <- eval(getCall(model)$data, environment(formula(model)))
},
"glm" = {
dat <- model$data
},
"glmmPQL" = {
dat <- model$data
},
"pgls" = {
dat <- model$data
},
"lmerMod" = {
dat <- lme4::getData(model) #model@frame
},
"glmerMod" = {
dat <- lme4::getData(model) #model@frame
},
"lmerModLmerTest" = {
dat <- lme4::getData(model) #model@frame
},
"glmmTMB" = {
dat <- model$frame
},
"gls" = {
dat <- nlme::getData(model)
},
"lme" = {
dat <- nlme::getData(model)
},
"gam" = {
dat <- eval(getCall(model)$data, environment(formula(model)))
}
)
return(dat)
}
#' Obtain (observation-level) random effects from a generalized linear mixed model
#'
#' RE = "all" all random effects are reported
#' RE = "RE" just group effects are reported
#' RE = "OLRE" just observation-level effects are reported
#'
#' @keywords internal
#'
GetOLRE <- function(sigma, model, X, data, RE = c("all", "RE", "OLRE")) {
if(class(model) %in% c("lmerMod", "glmerMod")) {
if(is.null(X)) X <- model.matrix(model)
rand <- sapply(lme4::findbars(formula(model)), function(x) as.character(x)[3])
rand <- rand[!duplicated(rand)]
}
# else if(class(model) %in% c("lme", "glmmPQL")) { }
idx <- sapply(sapply(strsplit(rand, "\\:"), function(x) gsub("\\(|\\)", "", x)), function(x) {
length(unique(data[, x])) == nrow(data)
} )
sigma.names <- unlist(names(sigma)) # unlist(strsplit(names(sigma), "\\."))
idx. <- sapply(sigma.names, function(x) !any(x %in% rand[idx]))
if(RE == "RE")
out <- sapply(sigma[idx.], function(i) {
if(all(rownames(i) %in% colnames(X))) X. <- X else
X. <- do.call(cbind, model.matrix(model, type = "randomListRaw"))
Z <- as.matrix(X.[, rownames(i), drop = FALSE])
sum(rowSums(Z %*% i) * Z) / nrow(X.)
} ) else if(RE == "OLRE") {
if(all(idx == FALSE)) out <- 0 else {
out <- sapply(sigma[idx], function(i) {
Z <- as.matrix(X[, rownames(i), drop = FALSE])
sum(rowSums(Z %*% i) * Z) / nrow(X)
} ) } } else if(RE == "all")
out <- sapply(sigma, function(i) {
Z <- as.matrix(X[, rownames(i), drop = FALSE])
sum(rowSums(Z %*% i) * Z) / nrow(X)
} )
if(length(out) == 0) out <- 0
return(out)
}
#' Get random effects variance-covariance from lme
#'
#' @keywords internal
#'
GetVarCov <- function(model) {
vc <- try(getVarCov(model), silent = TRUE)
if(any(class(vc) == "try-error")) {
vc <- nlme::VarCorr(model)
v <- suppressWarnings(as.numeric(vc[, 1]))
names(v) <- gsub(" =", "", rownames(vc))
vm <- as.list(na.omit(v[-length(v)]))
vl <- lapply(1:length(vm), function(i) matrix(vm[[i]], dimnames = list(names(vm)[i], names(vm)[i])))
names(vl) <- names(which(is.na(v)))
vl
} else list(vc)
}
#' Assess significance
#'
#' @keywords internal
#'
isSig <- function(p) {
ifelse(p > 0.01 & p < 0.05, "*",
ifelse(p > 0.001 & p <= 0.01, "**",
ifelse(p <= 0.001, "***", "")))
}
#' Recompute P-values using Kenward-Rogers approximation
#'
#' @keywords internal
#'
# KRp <- function(model, vars, data, intercepts = FALSE) {
#
# # if(any(grepl("\\*", all_vars_notrans(formula(model)))) & !all(grepl("\\*", vars))) {
#
# f <- all_vars_trans(formula(model))
#
# model <- update(model, as.formula(paste(f[1], " ~ ", paste(f[-1], collapse = " + "), " + ", paste(onlyBars(formula(model)), collapse = " + "))))
#
# # }
#
# out <- data.frame()
#
# for(x in vars) { #sapply(vars, function(x) {
#
# reduceModel <- update(model, as.formula(paste(". ~ . -", x)))
#
# if(nobs(model) != nobs(reduceModel)) stop("Different sample sizes for `KRmodcomp`. Remove all NAs and re-run")
#
# kr <- try(pbkrtest::KRmodcomp(model, reduceModel), silent = TRUE)
#
# if(class(kr) == "try-error")
#
# stop("Cannot obtain P-values from `lmerMod` using `pbkrtest::KRmodcopm`. Consider fitting using `nlme::lme`") else {
#
# d <- round(kr$stats$ddf, 2)
#
# p <- kr$stats$p.valueU
#
# out <- rbind(out, data.frame(d, p))
#
# }
#
# } # )
#
# if(intercepts == TRUE) {
#
# reduceModelI <- update(model, as.formula(paste("~ . - 1")), data = data)
#
# krI <- try(pbkrtest::KRmodcomp(model, reduceModelI), silent = TRUE)
#
# if(class(krI) == "try-error")
#
# stop("Cannot obtain P-values from `lmerMod` using `pbkrtest::KRmodcomp`. Consider re-fitting using `nlme::lme`")else {
#
# dI <- krI$stats$ddf
#
# pI <- krI$stats$p.valueU
#
# out <- rbind(data.frame(d = dI, p = pI), out)
#
# }
#
# }
#
# return(out)
#
# }
#' Get list of formula from a `psem` object
#'
#' @keywords internal
#'
listFormula <- function(modelList, formulas = 0) {
modelList <- removeData(modelList, formulas)
if(!all(class(modelList) %in% c("psem", "list"))) modelList <- list(modelList)
fList <- lapply(modelList, function(i) if(any(class(i) %in% c("formula.cerror"))) i else formula(i) )
fList <- lapply(fList, lme4::nobars)
return(fList)
}
#' Get number of observations from a model
#'
#' @keywords internal
#'
nObs <- function(object, ...) if(any(class(object) %in% c("phylolm", "phyloglm", "Sarlm"))) length(fitted(object)) else nobs(object, ...)
#' Get random effects from merMod
#'
#' @keywords internal
#'
onlyBars <- function(formula., slopes = TRUE) {
f <- lme4::findbars(formula.)
if(slopes == TRUE) paste(sapply(f, function(x) paste0("(", deparse(x), ")")), collapse = " + ") else {
# paste(sapply(f, function(x) paste0("(1 ", gsub(".*(\\|.*)", "\\1", f), ")")), collapse = "+")
f <- f[sapply(f, function(x) grepl("1\\||1 \\|", deparse(x)))]
paste(sapply(f, function(x) paste0("(", deparse(x), ")")), collapse = " + ")
}
}
#' Do not print attributes with custom functions
#'
#' @keywords internal
#'
#' @method print attr
#'
print.attr <- function(x, ...) {
attributes(x) <- NULL
noquote(x)
}
#' Remove data from the model list
#'
#' formulas = 0, keep everything
#' formulas = 1, remove all formulas including correlated errors
#' formulas = 2, remove only formula but keep correlated errors
#' formulas = 3, remove correlated errors but keep formula
#'
#' @keywords internal
#'
removeData <- function(modelList, formulas = 0) {
remove <- c("character", "matrix", "data.frame", "SpatialPointsDataFrame", "comparative.data")
if(formulas == 1) remove <- c(remove, "formula", "formula.cerror")
if(formulas == 2) remove <- c(remove, "formula")
if(formulas == 3) remove <- c(remove, "formula.cerror")
modelList[!sapply(modelList, function(x) any(class(x) %in% remove))]
}
#' Strip transformations
#'
#' @keywords internal
#'
stripTransformations <- function(x) {
x <- gsub(".*\\((.*)\\).*", "\\1", x)
gsub(" ", "", gsub("(.*)\\+.*", "\\1", x))
}
#' Get Response Name as a Character
#'
#' @keywords internal
#'
get_response <- function(mod) {
mod <- removeData(mod)
f <- lapply(mod, formula)
r <- lapply(f, function(x) x[[2]])
return(as.character(r))
}
#' Get Left-hand side of formulae
#'
#' @keywords internal
#'
getLHS <- function(formulaList){
sapply(formulaList, function(x) as.character(x[[2]]))
}
#' Get Right-hand side of formulae
#'
#' @keywords internal
#'
getRHS <- function(formulaList){
rhs <- sapply(formulaList, function(x) all.vars(x)[-1])
unique(do.call(c, rhs))
}
#' Operator for non-overlap in sets
#'
#' @keywords internal
#'
"%not_in%" <- function(x, y) x[!x %in% y]
#' Get a sorted psem object in DAG order
#'
#' @description Takes a [psem] object, pulls out the
#' DAG, and then sorts the psem object into the order
#' of the DAG (from exogenous to terminal endogenous
#' variable) for use by other functions. Note: removes
#' correlated errors.
#'
#' @param object A fit [psem] object
#' @param keepdata Defaults to TRUE. Should the
#' data with the psem be included in the returned
#' object?
#'
#' @return A new [psem] object, without the data.
#' @export
getSortedPsem <- function(object, keepdata = TRUE){
#first, remove data
dat <- object$data
object <- removeData(object, formulas = 1)
#Now, get formulae
formulaList <- listFormula(object)
lhs <- getLHS(formulaList)
names(object)<- lhs
#sort some dags so we do things in the right order
object_dag <- getDAG(formulaList)
sorted_dag <- sortDag(object_dag, formulaList)
lhs_sorted <- colnames(sorted_dag)
lhs_sorted <- lhs_sorted[which(lhs_sorted %in% lhs)]
#Sort the object
object <- object[lhs_sorted]
#should we include the data?
if(keepdata) object$data <- dat
#return
return(object)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.