Nothing
#' R-Squared Measures for Multilevel and Linear Mixed Effects Models by Rights
#' and Sterba (2019), Manually Inputting Parameter Estimates
#'
#' This function computes R-squared measures by Rights and Sterba (2019) for
#' multilevel and linear mixed effects models by manually inputting parameter
#' estimates.
#'
#' @param data a matrix or data frame with the level-1 and level-2 predictors
#' and outcome variable used in the model.
#' @param within a character vector with the variable names in \code{data} or
#' numeric vector with numbers corresponding to the columns in
#' \code{data} of the level-1 predictors used in the model. If
#' none used, set to \code{NULL}.
#' @param between a character vector with the variable names in \code{data} or
#' numeric vector with numbers corresponding to the columns in
#' \code{data} of the level-2 predictors used in the model. If
#' none used, set to \code{NULL}.
#' @param random a character vector with the variable names in \code{data} or
#' numeric vector with numbers corresponding to the columns in
#' \code{data} of the level-1 predictors that have random slopes
#' in the model. If no random slopes specified, set to \code{NULL}.
#' @param gamma.w a numeric vector of fixed slope estimates for all level-1
#' predictors, to be entered in the order of the predictors
#' listed in the argument \code{within}.
#' @param gamma.b a numeric vector of the intercept and fixed slope estimates
#' for all level-2predictors, to be entered in the order of the
#' predictors listed in the argument \code{between}. Note that
#' the first element is the parameter estimate for the intercept
#' if \code{intercept = TRUE}.
#' @param tau a matrix indicating the random effects covariance matrix, the
#' first row/column denotes the intercept variance and covariances
#' (if intercept is fixed, set all to 0) and each subsequent
#' row/column denotes a given random slope's variance and covariances
#' (to be entered in the order listed in the argument \code{random}).
#' @param sigma2 a numeric value indicating the level-1 residual variance.
#' @param intercept logical: if \code{TRUE} (default), the first element in the
#' \code{gamma.b} is assumed to be the fixed intercept estimate;
#' if set to \code{FALSE}, the first element in the argument
#' \code{gamma.b} is assumed to be the first fixed level-2
#' predictor slope.
#' @param center logical: if \code{TRUE} (default), all level-1 predictors are
#' assumed to be cluster-mean-centered and the function will
#' output all decompositions; if set to \code{FALSE}, function
#' will output only the total decomposition.
#' @param digits an integer value indicating the number of decimal places to be
#' used.
#' @param plot logical: if \code{TRUE}, bar chart showing the decomposition
#' of scaled total, within-cluster, and between-cluster outcome
#' variance into five (total), three (within-cluster), and two
#' (between-cluster) proportions is drawn. Note that the \pkg{ggplot2}
#' package is required to draw the bar chart.
#' @param gray logical: if \code{TRUE}, graphical parameter to draw the bar
#' chart in gray scale.
#' @param start a numeric value between 0 and 1, graphical parameter to specify
#' the gray value at the low end of the palette.
#' @param end a numeric value between 0 and 1, graphical parameter to specify
#' the gray value at the high end of the palette.
#' @param color a character vector, graphical parameter indicating the color
#' of bars in the bar chart in the following order: Fixed slopes
#' (Within), Fixed slopes (Between), Slope variation (Within),
#' Intercept variation (Between), and Residual (Within). By default,
#' colors from the colorblind-friendly palettes are used.
#' @param filename a character string indicating the \code{filename}
#' argument including the file extension in the \code{ggsave}
#' function. Note that one of \code{".eps"}, \code{".ps"},
#' \code{".tex"}, \code{".pdf"} (default),
#' \code{".jpeg"}, \code{".tiff"}, \code{".png"},
#' \code{".bmp"}, \code{".svg"} or \code{".wmf"} needs
#' to be specified as file extension in the \code{file}
#' argument. Note that plots can only be saved when
#' \code{plot = TRUE}.
#' @param width a numeric value indicating the \code{width} argument
#' (default is the size of the current graphics device)
#' in the \code{ggsave} function.
#' @param height a numeric value indicating the \code{height} argument
#' (default is the size of the current graphics device)
#' in the \code{ggsave} function.
#' @param units a character string indicating the \code{units} argument
#' (default is \code{in}) in the \code{ggsave} function.
#' @param dpi a numeric value indicating the \code{dpi} argument
#' (default is \code{600}) in the \code{ggsave} function.
#' @param write a character string naming a text file with file extension
#' \code{".txt"} (e.g., \code{"Output.txt"}) for writing the
#' output into a text file.
#' @param write a character string naming a text file with file extension
#' \code{".txt"} (e.g., \code{"Output.txt"}) for writing the
#' output into a text file.
#' @param append logical: if \code{TRUE} (default), output will be appended
#' to an existing text file with extension \code{.txt} specified
#' in \code{write}, if \code{FALSE} existing text file will be
#' overwritten.
#' @param check logical: if \code{TRUE} (default), argument specification is
#' checked.
#' @param output logical: if \code{TRUE}, (default) output is shown on the console.
#'
#' @details
#' A number of R-squared measures for multilevel and linear mixed effects models
#' have been developed in the methodological literature (see Rights & Sterba, 2018).
#' R-squared measures by Rights and Sterba (2019) provide an integrative framework
#' of R-squared measures for multilevel and linear mixed effects models with random
#' intercepts and/or slopes. Their measures are based on partitioning model implied
#' variance from a single fitted model, but they provide a full partitioning of
#' the total outcome variance to one of five specific sources. See the help page
#' of the \code{\link{multilevel.r2}} function for more details.
#'
#' @author
#' Jason D. Rights, Sonya K. Sterba, Jessica K. Flake, and Takuya Yanagida
#'
#' @seealso
#' \code{\link{multilevel.r2}}, \code{\link{multilevel.cor}},
#' \code{\link{multilevel.descript}}, \code{\link{multilevel.icc}},
#' \code{\link{multilevel.indirect}}
#'
#' @references
#' Rights, J. D., & Cole, D. A. (2018). Effect size measures for multilevel models
#' in clinical child and adolescent research: New r-squared methods and recommendations.
#' \emph{Journal of Clinical Child and Adolescent Psychology, 47}, 863-873.
#' https://doi.org/10.1080/15374416.2018.1528550
#'
#' Rights, J. D., & Sterba, S. K. (2019). Quantifying explained variance in multilevel
#' models: An integrative framework for defining R-squared measures. \emph{Psychological Methods, 24},
#' 309-338. https://doi.org/10.1037/met0000184
#'
#' @return
#' Returns an object of class \code{misty.object}, which is a list with following
#' entries:
#' \tabular{ll}{
#' \code{call} \tab function call \cr
#' \code{type} \tab type of analysis \cr
#' \code{data} \tab matrix or data frame specified in \code{data} \cr
#' \code{plot} \tab ggplot2 object for plotting the results \cr
#' \code{args} \tab specification of function arguments \cr
#' \code{result} \tab list with result tables \cr
#' }
#'
#' @note
#' This function is based on a copy of the function \code{r2mlm_manual()} in the
#' \pkg{r2mlm} package by Mairead Shaw, Jason Rights, Sonya Sterba, and Jessica
#' Flake.
#'
#' @export
#'
#' @examples
#' \dontrun{
#' # Load misty and lme4 package
#' misty::libraries(misty, lme4)
#'
#' # Load data set "Demo.twolevel" in the lavaan package
#' data("Demo.twolevel", package = "lavaan")
#'
#' #----------------------------------------------------------------------------
#'
#' # Cluster mean centering, center() from the misty package
#' Demo.twolevel <- center(Demo.twolevel, x2, type = "CWC", cluster = "cluster")
#'
#' # Compute group means, cluster.scores() from the misty package
#' Demo.twolevel <- cluster.scores(Demo.twolevel, x2, cluster = "cluster", name = "x2.b")
#'
#' # Estimate random intercept model using the lme4 package
#' mod1 <- lmer(y1 ~ x2.c + x2.b + w1 + (1| cluster), data = Demo.twolevel,
#' REML = FALSE, control = lmerControl(optimizer = "bobyqa"))
#'
#' # Estimate random intercept and slope model using the lme4 package
#' mod2 <- lmer(y1 ~ x2.c + x2.b + w1 + (1 + x2.c | cluster), data = Demo.twolevel,
#' REML = FALSE, control = lmerControl(optimizer = "bobyqa"))
#'
#' #----------------------------------------------------------------------------
#' # Example 1: Random intercept model
#'
#' # Fixed slope estimates
#' fixef(mod1)
#'
#' # Random effects variance-covariance matrix
#' as.data.frame(VarCorr(mod1))
#'
#' # R-squared measures according to Rights and Sterba (2019)
#' multilevel.r2.manual(data = Demo.twolevel,
#' within = "x2.c", between = c("x2.b", "w1"),
#' gamma.w = 0.41127956,
#' gamma.b = c(0.01123245, -0.08269374, 0.17688507),
#' tau = 0.9297401,
#' sigma2 = 1.813245794)
#'
#' #----------------------------------------------------------------------------
#' # Example 2: Random intercept and slope model
#'
#' # Fixed slope estimates
#' fixef(mod2)
#'
#' # Random effects variance-covariance matrix
#' as.data.frame(VarCorr(mod2))
#'
#' # R-squared measures according to Rights and Sterba (2019)
#' multilevel.r2.manual(data = Demo.twolevel,
#' within = "x2.c", between = c("x2.b", "w1"), random = "x2.c",
#' gamma.w = 0.41127956,
#' gamma.b = c(0.01123245, -0.08269374, 0.17688507),
#' tau = matrix(c(0.931008649, 0.004110479, 0.004110479, 0.017068857), ncol = 2),
#' sigma2 = 1.813245794)
#' }
multilevel.r2.manual <- function(data, within = NULL, between = NULL, random = NULL,
gamma.w = NULL, gamma.b = NULL, tau, sigma2,
intercept = TRUE, center = TRUE, digits = 3,
plot = FALSE, gray = FALSE, start = 0.15, end = 0.85,
color = c("#D55E00", "#0072B2", "#CC79A7", "#009E73", "#E69F00"),
filename = NULL, width = NA, height = NA,
units = c("in", "cm", "mm", "px"), dpi = 600,
write = NULL, append = TRUE, check = TRUE, output = TRUE) {
#_____________________________________________________________________________
#
# Initial Check --------------------------------------------------------------
# Check if input 'data' is missing
if (isTRUE(missing(data))) { stop("Please specify a matrix or data frame for the argument 'data'.", call. = FALSE) }
# Check if input 'data' is NULL
if (isTRUE(is.null(data))) { stop("Input specified for the argument 'data' is NULL.", call. = FALSE) }
# Check if input 'tau' is missing
if (isTRUE(missing(tau))) { stop("Please specify a matrix for the argument 'tau'.", call. = FALSE) }
# Check if input 'tau' is NULL
if (isTRUE(is.null(tau))) { stop("Input specified for the argument 'tau' is NULL.", call. = FALSE) }
# Check if input 'sigma2' is missing
if (isTRUE(missing(sigma2))) { stop("Please specify a numeric value for the argument 'sigma2'.", call. = FALSE) }
# Check if input 'sigma2' is NULL
if (isTRUE(is.null(sigma2))) { stop("Input specified for the argument 'sigma2' is NULL.", call. = FALSE) }
#_____________________________________________________________________________
#
# Input Check ----------------------------------------------------------------
# Check inputs
.check.input(logical = c("intercept", "center", "plot", "gray", "append", "output"),
numeric = list(start = 1L, end = 1L),
args = c("digits", "write1"),
envir = environment(), input.check = check)
# Additional checks
if (isTRUE(check)) {
# Package ggplot2
if (isTRUE(plot)) { if (isTRUE(!nzchar(system.file(package = "ggplot2")))) { stop("Package \"ggplot2\" is needed to draw a plot, please install the package.", call. = FALSE) } }
#......
# Check input 'within'
# Character vector
if (isTRUE(is.character(within))) {
# Check if level-1 predictors are in the data
(!within %in% colnames(data)) |> (\(y) if (isTRUE(any(y))) { stop(paste0("Predictors specified in the argument 'within' were not found in 'data': ", paste(within[which(y)], collapse = ", ")), call. = FALSE) })()
} else {
if (isTRUE(ncol(data) < max(within))) { stop("Colummn numbers specified in the argument 'within' were not found in 'data'", call. = FALSE) }
}
#......
# Check input 'between'
# Character vector
if (isTRUE(is.character(between))) {
# Check if level-2 predictors are in the data
(!between %in% colnames(data)) |> (\(y) if (isTRUE(any(y))) { stop(paste0("Predictors specified in the argument 'between' were not found in 'data': ", paste(between[which(y)], collapse = ", ")), call. = FALSE) })()
} else {
if (isTRUE(ncol(data) < max(between))) { stop("Colummn numbers specified in the argument 'between' were not found in 'data'", call. = FALSE) }
}
# Check input 'random'
# Character vector
if (isTRUE(is.character(between))) {
# Check if level-1 predictors with random slopes are in 'within'
(!random %in% within) |> (\(y) if (isTRUE(any(y))) { stop(paste0("Predictors specified in the argument 'random' were not found in 'within': ", paste(between[which(y)], collapse = ", ")), call. = FALSE) })()
} else {
if (isTRUE(any(!random %in% within))) { stop("Colummn numbers specified in the argument 'random' were not all found in 'within'", call. = FALSE) }
}
# Check input 'gamma.w'
if (isTRUE(length(gamma.w) != length(within))) { stop("The number of parameter estimates specified in 'gamma.w' does not match with 'within'", call. = FALSE) }
# Check input 'gamma.b'
if (isTRUE(intercept)) {
if (isTRUE(length(gamma.b) != (length(between) + 1))) { stop("The number of parameter estimates specified in 'gamma.b' does not match with the length of 'between' plus one for the intercept.", call. = FALSE) }
} else {
if (isTRUE(length(gamma.b) != length(between))) { stop("The number of parameter estimates specified in 'gamma.b' does not match with 'between'.", call. = FALSE) }
}
# Check input 'tau'
if (isTRUE(!is.null(random))) {
if (isTRUE(intercept)) {
if (ncol(tau) != (length(random) + 1L)) { stop("The matrix specified in 'tau' does not match the number of level-1 predictors with random slopes plus one for the intercept.", call. = FALSE) }
} else {
if (ncol(tau) != (length(random))) { stop("The matrix specified in 'tau' does not match the number of level-1 predictors with random slopes.", call. = FALSE) }
}
} else {
if (isTRUE(intercept)) {
if (isTRUE(length(tau) != 1L)) { stop("Please specify a numeric value for the random intercept variance.", call. = FALSE) }
}
}
# Check input 'sigma2'
if (isTRUE(length(sigma2) != 1L || !is.numeric(sigma2))) { stop("Please specify a numeric value for the argument 'sigma2'.", call. = FALSE) }
# Check input 'start'
if (isTRUE(start < 0L || start > 1L)) { stop("Please specify a numeric value between 0 and 1 for the argument 'start'", call. = FALSE) }
# Check input 'end'
if (isTRUE(end < 0L || end > 1L)) { stop("Please specify a numeric value between 0 and 1 for the argument 'end'", call. = FALSE) }
}
#_____________________________________________________________________________
#
# Data and Variables ---------------------------------------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## 'within', 'between', and 'random' Arguments ####
if (isTRUE(is.character(within))) { within_covs <- match(within, colnames(data)) } else { within_covs <- within }
if (isTRUE(is.character(between))) { between_covs <- match(between, colnames(data)) } else { between_covs <- between }
if (isTRUE(is.character(random))) { random_covs <- match(random, colnames(data)) } else { random_covs <- random }
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## 'units' Argument ####
# Default setting
if (isTRUE(all(c("in", "cm", "mm", "px") %in% units))) { units <- "in" }
#_____________________________________________________________________________
#
# Main Function --------------------------------------------------------------
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Internal functions from the r2mlm package ####
#### r2mlm_manual() Function ####
r2mlm_manual <- function(data, within_covs, between_covs, random_covs, gamma_w, gamma_b, Tau, sigma2, has_intercept = TRUE, clustermeancentered = TRUE) {
if (isTRUE(has_intercept)) {
if (isTRUE(length(gamma_b) > 1L)) gamma <- c(1L, gamma_w, gamma_b[2:length(gamma_b)])
if (isTRUE(length(gamma_b) == 1L)) gamma <- c(1L, gamma_w)
if (isTRUE(is.null(within_covs))) gamma_w <- 0L
}
if (isTRUE(!has_intercept)) {
gamma <- c(gamma_w, gamma_b)
if (isTRUE(is.null(within_covs))) gamma_w <- 0L
if (isTRUE(is.null(between_covs))) gamma_b <- 0L
}
if (isTRUE(is.null(gamma))) gamma <- 0L
# Compute phi
phi <- var(cbind(1L, data[, c(within_covs)], data[, c(between_covs)]), na.rm = TRUE)
if (isTRUE(!has_intercept)) phi <- var(cbind(data[, c(within_covs)], data[, c(between_covs)]), na.rm = TRUE)
if (isTRUE(is.null(within_covs) && is.null(within_covs) && !has_intercept)) phi <- 0L
phi_w <- var(data[, within_covs], na.rm = TRUE)
if (isTRUE(is.null(within_covs))) phi_w <- 0L
phi_b <- var(cbind(1L, data[, between_covs]), na.rm = TRUE)
if (isTRUE(is.null(between_covs))) phi_b <- 0L
# Compute psi and kappa
var_randomcovs <- var(cbind(1, data[, c(random_covs)]), na.rm = TRUE)
if (isTRUE(length(Tau) > 1L)) psi <- matrix(c(diag(Tau)), ncol = 1L)
if (isTRUE(length(Tau) == 1L)) psi <- Tau
if (isTRUE(length(Tau) > 1L)) kappa <- matrix(c(Tau[lower.tri(Tau) == TRUE]), ncol = 1)
if (isTRUE(length(Tau) == 1L)) kappa <- 0L
v <- matrix(c(diag(var_randomcovs)), ncol = 1L)
r <- matrix(c(var_randomcovs[lower.tri(var_randomcovs) == TRUE]), ncol = 1L)
if (isTRUE(is.null(random_covs))) {
v <- 0L
r <- 0L
m <- matrix(1L, ncol = 1L)
}
if (isTRUE(length(random_covs) > 0L)) m <- matrix(c(colMeans(cbind(1, data[, c(random_covs)]), na.rm = TRUE)), ncol = 1L)
# Total variance
totalvar_notdecomp <- t(v)%*%psi + 2L*(t(r)%*%kappa) + t(gamma)%*%phi%*%gamma + t(m)%*%Tau%*%m + sigma2
totalwithinvar <- (t(gamma_w)%*%phi_w%*%gamma_w) + (t(v)%*%psi + 2L*(t(r)%*%kappa)) + sigma2
totalbetweenvar <- (t(gamma_b)%*%phi_b%*%gamma_b) + Tau[1L]
totalvar <- totalwithinvar + totalbetweenvar
# Total decomp
decomp_fixed_notdecomp <- (t(gamma)%*%phi%*%gamma) / totalvar_notdecomp
decomp_varslopes_notdecomp <- (t(v)%*%psi + 2L*(t(r)%*%kappa)) / totalvar_notdecomp
decomp_varmeans_notdecomp <- (t(m)%*%Tau%*%m) / totalvar_notdecomp
decomp_sigma_notdecomp <- sigma2/totalvar_notdecomp
decomp_fixed_within <- (t(gamma_w)%*%phi_w%*%gamma_w) / totalvar
decomp_fixed_between <- (t(gamma_b)%*%phi_b%*%gamma_b) / totalvar
decomp_fixed <- decomp_fixed_within + decomp_fixed_between
decomp_varslopes <- (t(v)%*%psi + 2L*(t(r)%*%kappa)) / totalvar
decomp_varmeans <- (t(m)%*%Tau%*%m) / totalvar
decomp_sigma <- sigma2 / totalvar
# Within decomp
decomp_fixed_within_w <- (t(gamma_w)%*%phi_w%*%gamma_w) / totalwithinvar
decomp_varslopes_w <- (t(v)%*%psi + 2L*(t(r)%*%kappa)) / totalwithinvar
decomp_sigma_w <- sigma2 / totalwithinvar
# Between decomp
decomp_fixed_between_b <- (t(gamma_b)%*%phi_b%*%gamma_b) / totalbetweenvar
decomp_varmeans_b <- Tau[1L] / totalbetweenvar
# New measures
if (isTRUE(clustermeancentered)) {
R2_f <- decomp_fixed
R2_f1 <- decomp_fixed_within
R2_f2 <- decomp_fixed_between
R2_fv <- decomp_fixed + decomp_varslopes
R2_fvm <- decomp_fixed + decomp_varslopes + decomp_varmeans
R2_v <- decomp_varslopes
R2_m <- decomp_varmeans
R2_f_w <- decomp_fixed_within_w
R2_f_b <- decomp_fixed_between_b
R2_fv_w <- decomp_fixed_within_w + decomp_varslopes_w
R2_v_w <- decomp_varslopes_w
R2_m_b <- decomp_varmeans_b
}
if (isTRUE(!clustermeancentered)) {
R2_f <- decomp_fixed_notdecomp
R2_fv <- decomp_fixed_notdecomp + decomp_varslopes_notdecomp
R2_fvm <- decomp_fixed_notdecomp + decomp_varslopes_notdecomp + decomp_varmeans_notdecomp
R2_v <- decomp_varslopes_notdecomp
R2_m <- decomp_varmeans_notdecomp
}
if (isTRUE(clustermeancentered)) {
decomp_table <- matrix(c(decomp_fixed_within, decomp_fixed_between, decomp_varslopes, decomp_varmeans, decomp_sigma,
decomp_fixed_within_w, "NA", decomp_varslopes_w, "NA", decomp_sigma_w,
"NA", decomp_fixed_between_b, "NA", decomp_varmeans_b, "NA"), ncol = 3L)
decomp_table <- suppressWarnings(apply(decomp_table, 2, as.numeric))
rownames(decomp_table) <- c("fixed, within", "fixed, between", "slope variation", "mean variation", "sigma2")
colnames(decomp_table) <- c("total", "within", "between")
R2_table <- matrix(c(R2_f1, R2_f2, R2_v, R2_m, R2_f, R2_fv, R2_fvm,
R2_f_w, "NA", R2_v_w, "NA", "NA", R2_fv_w, "NA",
"NA", R2_f_b, "NA", R2_m_b, "NA", "NA", "NA"), ncol = 3L)
R2_table <- suppressWarnings(apply(R2_table, 2, as.numeric)) # make values numeric, not character
rownames(R2_table) <- c("f1", "f2", "v", "m", "f", "fv", "fvm")
colnames(R2_table) <- c("total", "within", "between")
}
if (isTRUE(!clustermeancentered)) {
decomp_table <- matrix(c(decomp_fixed_notdecomp, decomp_varslopes_notdecomp, decomp_varmeans_notdecomp, decomp_sigma_notdecomp), ncol = 1L)
decomp_table <- suppressWarnings(apply(decomp_table, 2L, as.numeric))
rownames(decomp_table) <- c("fixed", "slope variation", "mean variation", "sigma2")
colnames(decomp_table) <- c("total")
R2_table <- matrix(c(R2_f, R2_v, R2_m, R2_fv, R2_fvm), ncol = 1L)
R2_table <- suppressWarnings(apply(R2_table, 2L, as.numeric))
rownames(R2_table) <- c("f", "v", "m", "fv", "fvm")
colnames(R2_table) <- c("total")
}
Output <- list(noquote(decomp_table), noquote(R2_table))
names(Output) <- c("Decompositions", "R2s")
return(Output)
}
#### prepare_data() Function ####
prepare_data <- function(model, calling_function, cluster_variable, second_model = NULL) {
# Step 1a: Pull dataframe associated with model
if (isTRUE(calling_function == "lme4")) {
data <- model@frame
} else {
data <- model[["data"]]
}
# Step 2a) pull interaction terms into list
interaction_vars <- get_interaction_vars(model)
if (isTRUE(!is.null(second_model))) {
interaction_vars_2 <- get_interaction_vars(second_model)
interaction_vars <- unique(append(interaction_vars, interaction_vars_2))
}
# Step 2b) split interaction terms into halves, multiply halves to create new columns in dataframe
data <- add_interaction_vars_to_data(data, interaction_vars)
return(data)
}
#### add_interaction_vars_to_data() Function ####
add_interaction_vars_to_data <- function(data, interaction_vars) {
for (whole in interaction_vars) {
half1 <- unlist(strsplit(whole, ":"))[1L]
half2 <- unlist(strsplit(whole, ":"))[2L]
newcol <- data[[half1]] * data[[half2]]
data <- within(data, assign(whole, newcol))
}
return(data)
}
#### get_covs() Function ####
get_covs <- function(variable_list, data) {
cov_list <- c()
i <- 1L
for (variable in variable_list) {
tmp <- match(variable, names(data))
cov_list[i] <- tmp
i <- i + 1L
}
return(cov_list)
}
#### get_random_slope_vars() Function ####
get_random_slope_vars <- function(model, has_intercept, calling_function) {
# Visible global function defnition
ranef <- NULL
if (isTRUE(calling_function == "lme4")) {
temp_cov_list <- ranef(model)[[1L]]
} else if (calling_function == "nlme") {
temp_cov_list <- nlme::ranef(model)
}
if (isTRUE(has_intercept == 1L)) {
running_count <- 2L
} else {
running_count <- 1L
}
random_slope_vars <- c()
x <- 1L
while (running_count <= length(temp_cov_list)) {
random_slope_vars[x] <- names(temp_cov_list[running_count])
x <- x + 1L
running_count <- running_count + 1L
}
return(random_slope_vars)
}
#### get_cwc() Function ####
get_cwc <- function(l1_vars, cluster_variable, data) {
for (variable in l1_vars) {
t <- tapply(data[, variable], data[, cluster_variable], sum, na.rm = TRUE)
temp_tracker <- 0L
# sum all of the sums
for (i in t) { temp_tracker <- temp_tracker + i }
if (abs(temp_tracker) < 0.0000001) {
centeredwithincluster <- TRUE
} else {
centeredwithincluster <- FALSE
break
}
}
return(centeredwithincluster)
}
#### get_interaction_vars() Function ####
get_interaction_vars <- function(model) {
interaction_vars <- c()
x <- 1L
for (term in attr(terms(model), "term.labels")) {
if (isTRUE(grepl(":", term))) {
interaction_vars[x] <- term
x <- x + 1L
}
}
return(interaction_vars)
}
#### sort_variables() Function ####
sort_variables <- function(data, predictors, cluster_variable) {
l1_vars <- c()
l2_vars <- c()
l1_counter <- 1L
l2_counter <- 1L
for (variable in predictors) {
t <- tapply(data[, variable], data[, cluster_variable], var, na.rm = TRUE)
counter <- 1L
while (counter <= length(t)) {
if (is.na(t[[counter]])) { t[[counter]] <- 0L }
counter <- counter + 1L
}
variance_tracker <- 0L
for (i in t) { variance_tracker <- variance_tracker + i }
if (isTRUE(variance_tracker == 0L)) {
l2_vars[l2_counter] <- variable
l2_counter <- l2_counter + 1L
} else {
l1_vars[l1_counter] <- variable
l1_counter <- l1_counter + 1L
}
}
return(list("l1_vars" = l1_vars, "l2_vars" = l2_vars))
}
#...................
### R-squared measures ####
r2mlm.out <- r2mlm_manual(data = data,
within_covs = within_covs, between_covs = between_covs,
random_covs = random_covs,
gamma_w = gamma.w, gamma_b = gamma.b,
Tau = tau, sigma2 = sigma2,
has_intercept = intercept, clustermeancentered = center)
rs <- suppressWarnings(list(decomp = matrix(apply(r2mlm.out$Decomposition, 2L, as.numeric),
ncol = ncol(r2mlm.out$Decomposition),
dimnames = list(rownames(r2mlm.out$Decompositions),
colnames(r2mlm.out$Decompositions))),
r2 = matrix(apply(r2mlm.out$R2s, 2L, as.numeric), ncol = ncol(r2mlm.out$R2s),
dimnames = list(rownames(r2mlm.out$R2s), colnames(r2mlm.out$R2s)))))
#_____________________________________________________________________________
#
# Return Object --------------------------------------------------------------
object <- list(call = match.call(),
type = "multilevel.r2.manual",
data = data,
args = list(within = within, between = between, random = random,
gamma.w = gamma.w, gamma.b = gamma.b, tau = tau,
sigma2 = sigma2, intercept = intercept, center = center,
digits = digits, plot = plot, gray = gray,
start = start, end = end, color = color,
width = width, height = height, units = units, dpi = dpi,
write = write, append = append, check = check, output = output),
plot = NULL,
result = list(decomp = rs$decomp,
total = data.frame(f1 = ifelse(ncol(rs$r2) > 1L, rs$r2[row.names(rs$r2) == "f1", "total"], NA),
f2 = ifelse(ncol(rs$r2) > 1L, rs$r2[row.names(rs$r2) == "f2", "total"], NA),
f = rs$r2[row.names(rs$r2) == "f", "total"],
v = rs$r2[row.names(rs$r2) == "v", "total"],
m = rs$r2[row.names(rs$r2) == "m", "total"],
fv = rs$r2[row.names(rs$r2) == "fv", "total"],
fvm = rs$r2[row.names(rs$r2) == "fvm", "total"]),
within = data.frame(f1 = ifelse(ncol(rs$r2) > 1L, rs$r2[row.names(rs$r2) == "f1", "within"], NA),
v = ifelse(ncol(rs$r2) > 1L, rs$r2[row.names(rs$r2) == "v", "within"], NA),
f1v = ifelse(ncol(rs$r2) > 1L, rs$r2[row.names(rs$r2) == "fv", "within"], NA)),
between = data.frame(f2 = ifelse(ncol(rs$r2) > 1L, rs$r2[row.names(rs$r2) == "f2", "between"], NA),
m = ifelse(ncol(rs$r2) > 1L, rs$r2[row.names(rs$r2) == "m", "between"], NA))))
class(object) <- "misty.object"
#_____________________________________________________________________________
#
# Plot and Save Results ------------------------------------------------------
if (isTRUE(plot)) { object$plot <- plot(object, filename = filename, width = width, height = height, units = units, dpi = dpi, check = FALSE) |> (\(y) suppressMessages(suppressWarnings(print(y))))() }
#_____________________________________________________________________________
#
# Write Results --------------------------------------------------------------
if (isTRUE(!is.null(write))) {
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## Text file ####
# Send R output to textfile
sink(file = write, append = ifelse(isTRUE(file.exists(write)), append, FALSE), type = "output", split = FALSE)
if (isTRUE(append && file.exists(write))) { write("", file = write, append = TRUE) }
# Print object
print(object, check = FALSE)
# Close file connection
sink()
}
#_____________________________________________________________________________
#
# Output ---------------------------------------------------------------------
if (isTRUE(output)) { print(object, check = FALSE) }
return(invisible(object))
}
#_______________________________________________________________________________
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.