Nothing
# Adjusting formula for lmList
#
# When dealing with nested data, sometimes a factor will only take on one level
# within a group. This causes an error, since a factor must have 2 or more levels
# to be included as a predictor within a formula for a linear model. When a
# factor takes on only one level, the effect of that factor can then be included
# in the intercept of the model. This function takes a formula for use with \code{lmList}
# and will make the necessary adjustments to avoid errors while fitting the separate LS models.
#
# @param formula a linear formula that is used by \code{lmList()} e.g. y ~ x1 + x2 + ... + xn | g
# @param data the model frame from the model fit by \code{lmer()}
# @return a list containing the adjusted formulas
# @author Adam Loy \email{loyad01@@gmail.com}
# @keywords models regression
problem_factor_groups <- function(formula, data){
form <- formula(formula)
model_frame <- data
g <- deparse(form[[3]][[3]])
ngroups <- length(unique(data[,g]))
y <- deparse(form[[2]])
xs <- setdiff(names(model_frame), c(deparse(form[[2]]), g))
xs_frame <- model_frame[, xs]
are_factors <- rapply(xs_frame, class)
are_factors <- names(are_factors)[which(are_factors == "factor")]
id_vars <- c(g, are_factors)
molten_model_frame <- melt(model_frame, id = id_vars)
one_level_factor <- matrix(NA, nrow = ngroups, ncol = 1 + length(are_factors))
colnames(one_level_factor) <- c(g, are_factors)
for(i in 1:length(are_factors)){
cast_model_frame <- dcast(molten_model_frame, paste(g, "~", are_factors[i]), length)
if(i == 1) one_level_factor[,1] <- cast_model_frame[,1]
one_level_factor[,(i + 1)] <- apply(cast_model_frame[,-1], 1, function(x){
empty_factors <- sum(x == 0)
empty_factors > 0
})
}
return(one_level_factor)
}
subbars <- function(term)
#@ Substitute the '+' function for the '|' function (from lmer)
{
if (is.name(term) || !is.language(term)) return(term)
if (length(term) == 2) {
term[[2]] <- subbars(term[[2]])
return(term)
}
stopifnot(length(term) >= 3)
if (is.call(term) && term[[1]] == as.name('|'))
term[[1]] <- as.name('+')
for (j in 2:length(term)) term[[j]] <- subbars(term[[j]])
term
}
adjust_formula_lmList <- function(formula, data){
lm_form <- formula(formula)
lm_form[[3]] <- lm_form[[3]][[2]]
mf <- match.call()
m <- match(c("data", "subset", "weights", "na.action", "offset"), names(mf), 0)
mf <- mf[c(1, m)]
mf$formula <- subbars(formula)
mf$drop.unused.levels <- TRUE
mf[[1]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
problem_mat <- data.frame(problem_factor_groups(formula, mf))
problem_list <- split(problem_mat, problem_mat[,1])
formula_list <- lapply(problem_list, function(x){
#formula_list <- vector("list", length(problem_list))
#for(i in 1:length(problem_list)){
# x <- problem_list[[i]]
changes <- x[,-1]
if(dim(x)[2] == 2){
if(changes == 1){changes <- colnames(x)[2]
changes <- paste(changes, ":.", sep = "")
}
else changes <- NULL
}
else{
nchanges <- sum(changes)
if(nchanges == 0) changes <- NULL
else{
if(nchanges == 1){
changes <- colnames(changes)[which(changes == 1)]
changes <- paste(changes, ":.", sep = "")
}
else{
#if(nchanges > 1){
changes <- colnames(changes)[which(changes == 1)]
changes <- paste(changes, ":.", sep = "")
changes <- paste(changes, collapse = "-")
}
}
}
if(is.null(changes) != TRUE) {change_form <- update(lm_form, paste(". ~ .", changes, sep = "-"))}
else{ change_form <- lm_form }
#formula_list[[i]] <- change_form
})
return(formula_list)
}
pooledSD <- function(x, ...)
{
stopifnot(is(x, "adjust_lmList"))
sumsqr <- apply(sapply(x,
function(el) {
if (is.null(el)) {
c(0,0)
} else {
res <- resid(el)
c(sum(res^2), length(res) - length(coef(el)))
}
}), 1, sum)
if (sumsqr[2] == 0) {
stop("No degrees of freedom for estimating std. dev.")
}
val <- sqrt(sumsqr[1]/sumsqr[2])
attr(val, "df") <- sumsqr[2]
val
}
#' @export
adjust_lmList <- function(object, data, pool){
UseMethod("adjust_lmList", object)
}
#'Fitting Common Models via \code{lm}
#'
#'Separate linear models are fit via \code{lm} similar to \code{lmList},
#'however, \code{adjust_lmList} can handle models where a factor takes only one
#'level within a group. In this case, the \code{formula} is updated eliminating
#'the offending factors from the formula for that group as the effect is
#'absorbed into the intercept.
#'
#'
#'@aliases adjust_lmList adjust_lmList,formula,data.frame-method
#'@keywords models regression
#'@param object a linear formula such as that used by \code{lmList}, e.g.
#'\code{y ~ x1 + ... + xn | g}, where \code{g} is a grouping factor.
#'@param data a data frame containing the variables in the model.
#'@param pool a logical value that indicates whether the pooled standard
#'deviation/error should be used.
#'@seealso \code{\link[lme4]{lmList}, \link[stats]{lm}}
#'@references Douglas Bates, Martin Maechler and Ben Bolker (2012). lme4:
#'Linear mixed-effects models using S4 classes. R package version 0.999999-0.
#'@export
#'@examples
#'
#' data(Exam, package = 'mlmRev')
#' sepLM <- adjust_lmList(normexam ~ standLRT + sex + schgend | school, data = Exam)
#' confint(sepLM)
#'
adjust_lmList.formula <- function(object, data, pool){
# options(show.error.messages = FALSE)
# lmList_result <- try(lmList(formula = formula, data = data), silent = TRUE)
lmList_result <- suppressWarnings(lme4::lmList(formula = object, data = data))
# options(show.error.messages = TRUE)
orig_names <- names(lmList_result)
check_results <- unlist(lapply(lmList_result, is.null))
form <- formula(object)
g <- deparse(form[[3]][[3]])
ngroups <- length(unique(data[,g]))
if(sum(check_results) != 0){ #return(lmList_result)
#else{
new_formulas <- adjust_formula_lmList(object, data)
problem_cases <- as.numeric(which(check_results == TRUE))
#print(problem_cases)
split_data <- split(data, data[,g])
for(i in problem_cases){
lmList_result[[i]] <- lm(formula = new_formulas[[i]], data = split_data[[i]])
}
}
class(lmList_result) <- "adjust_lmList"
if(missing(pool)) pool <- TRUE
#lmList_result <- new("adjust_lmList", lmList_result, call = match.call(), pool = pool)
attr(lmList_result, "names") <- orig_names
attr(lmList_result,"call") <- match.call()
lmList_result
}
# setClass("adjust_lmList", representation(call = "call", pool = "logical"), contains = "list")
#
# setClass("adjust_lmList.confint", contains = "array")
# #' @export
# setMethod("adjust_lmList", signature(formula = "formula", data = "data.frame"),
# function(formula, data, pool){
# # options(show.error.messages = FALSE)
# # lmList_result <- try(lmList(formula = formula, data = data), silent = TRUE)
# lmList_result <- lmList(formula = formula, data = data)
# # options(show.error.messages = TRUE)
# orig_names <- names(lmList_result)
# check_results <- unlist(lapply(lmList_result, is.null))
# form <- formula(formula)
# g <- deparse(form[[3]][[3]])
# ngroups <- length(unique(data[,g]))
#
# if(sum(check_results) != 0){ #return(lmList_result)
#
# #else{
# new_formulas <- adjust_formula_lmList(formula, data)
# problem_cases <- as.numeric(which(check_results == TRUE))
# #print(problem_cases)
# split_data <- split(data, data[,g])
#
# for(i in problem_cases){
# lmList_result[[i]] <- lm(formula = new_formulas[[i]], data = split_data[[i]])
# }
#
# }
# #class(lmList_result) <- "adjust_lmList"
# if(missing(pool)) pool <- TRUE
# lmList_result <- new("adjust_lmList", lmList_result, call = match.call(), pool = pool)
# attr(lmList_result, "names") <- orig_names
# lmList_result
# })
#' @export
#' @method coef adjust_lmList
#' @importFrom plyr rbind.fill
# setMethod("coef", signature(object = "adjust_lmList"),
coef.adjust_lmList <- function(object, ...){
coefs <- lapply(object, coef)
non.null <- !unlist(lapply(coefs, is.null))
if(sum(non.null) > 0){
coefs <- lapply(coefs, function(x) as.data.frame(t(x)))
coefs <- do.call('rbind.fill', coefs)
rownames(coefs) <- names(object)
#dimnames(coefs) <- list(names(object), names(coefs))
#coefs <- as.data.frame(coefs)
#effectNames <- names(coefs)
}
coefs
}
# )
#' @export
#' @method print adjust_lmList
# setMethod("show", signature(object = "adjust_lmList"),
print.adjust_lmList <- function(x, ...){
cat("Call:", deparse(attr(x, "call")), "\n")
cat("Coefficients:\n")
invisible(print(coef(x)))
}
# )
#' @export
#' @method confint adjust_lmList
# setMethod("confint", signature(object = "adjust_lmList"),
confint.adjust_lmList <- function(object, parm, level = 0.95, pool = NULL, ...){
#if(length(object < 1))
# return(new("adjust_lmList.confint", array(numeric(0), c(0,0,0))))
coefs <- coef(object)
val <- array(NA, c(dim(coefs)[2], 2, length(object)) )
dimnames(val) <- list(names(coefs),
c(paste((1-level)/2 * 100, "%"), paste((1 - (1-level)/2) * 100, "%")),
names(object))
#val <- aperm(val, c(3, 2, 1))l
if(is.null(pool)) pool <- object@pool
if(length(pool) > 0 && pool[1]){
sd <- pooledSD(object)
a <- (1 - level)/2
fac <- sd * qt(c(a, 1 - a), attr(sd, "df"))
parm <- names(coefs)
for (i in seq_along(object)){
cis <- coef(object[[i]]) +
sqrt(diag(summary(object[[i]], corr = FALSE)$cov.unscaled)) %o% fac
for(j in names(coef(object[[i]]))){
val[j, , i] <- cis[j,]
}
}
}
else{
cis <- lapply(object, confint)
for (i in seq_along(object)) {
for(j in names(coef(object[[i]]))){
val[j, , i] <- cis[[i]][j,]
}
}
}
class(val) <- "adjust_lmList.confint"
val
} #, valueClass = "adjust_lmList.confint"
# )
#' @export
#' @method plot adjust_lmList.confint
# setMethod("plot", signature(x = "adjust_lmList.confint"),
plot.adjust_lmList.confint <- function(x, y, ...){
# stopifnot(require("ggplot2"))
group <- intervals <- NULL # Make codetools happy
cis <- as(x, "array")
df <- adply(cis, c(1, 2, 3), identity)
colnames(df) <- c("what", "end", "group", "intervals")
p <- ggplot(df, aes(x = group, y = intervals))
p + geom_line() +
geom_errorbar(aes(ymin = intervals, ymax = intervals)) +
facet_wrap( ~ what, nrow = 1, scales = "free")
# coord_flip() + ylab(NULL)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.