From 9c62111e666efa141454b90f18bb490c30a68f30 Mon Sep 17 00:00:00 2001 From: Daniel Date: Sun, 30 Jun 2024 17:30:36 +0200 Subject: [PATCH] Support for svy2lme (#985) * Support for svy2lme * fix * fix * news, desc * add tests * some fixes * fix tests * fix partial matching * wordlist * fix --- DESCRIPTION | 4 +- NAMESPACE | 4 + NEWS.md | 6 ++ R/extract_random_variances.R | 36 ++++++++ R/methods_BayesFactor.R | 9 +- R/methods_aov.R | 8 -- R/methods_htest.R | 7 -- R/methods_other.R | 7 -- R/methods_svy2lme.R | 113 ++++++++++++++++++++++++++ R/n_clusters_easystats.R | 4 +- R/utils_format.R | 8 +- inst/WORDLIST | 1 + man/model_parameters.BFBayesFactor.Rd | 3 - man/model_parameters.aov.Rd | 3 - man/model_parameters.cgam.Rd | 3 - man/model_parameters.htest.Rd | 3 - tests/testthat/_snaps/svylme.md | 26 ++++++ tests/testthat/test-rstanarm.R | 4 +- tests/testthat/test-svylme.R | 27 ++++++ 19 files changed, 225 insertions(+), 51 deletions(-) create mode 100644 R/methods_svy2lme.R create mode 100644 tests/testthat/_snaps/svylme.md create mode 100644 tests/testthat/test-svylme.R diff --git a/DESCRIPTION b/DESCRIPTION index ed19045b2..bf4f3bfc1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: parameters Title: Processing of Model Parameters -Version: 0.22.0 +Version: 0.22.0.1 Authors@R: c(person(given = "Daniel", family = "Lüdecke", @@ -197,6 +197,7 @@ Suggests: sparsepca, survey, survival, + svylme, testthat (>= 3.2.1), tidyselect, tinytable (>= 0.1.0), @@ -215,3 +216,4 @@ Config/testthat/edition: 3 Config/testthat/parallel: true Config/Needs/website: easystats/easystatstemplate Config/rcmdcheck/ignore-inconsequential-notes: true +Remotes: easystats/insight, easystats/bayestestR diff --git a/NAMESPACE b/NAMESPACE index e2cf309e1..1bfa0c66b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -170,6 +170,7 @@ S3method(degrees_of_freedom,rqss) S3method(degrees_of_freedom,selection) S3method(degrees_of_freedom,serp) S3method(degrees_of_freedom,summary.lm) +S3method(degrees_of_freedom,svy2lme) S3method(degrees_of_freedom,systemfit) S3method(degrees_of_freedom,truncreg) S3method(degrees_of_freedom,vgam) @@ -382,6 +383,7 @@ S3method(model_parameters,stanfit) S3method(model_parameters,stanmvreg) S3method(model_parameters,stanreg) S3method(model_parameters,summary_emm) +S3method(model_parameters,svy2lme) S3method(model_parameters,svyglm) S3method(model_parameters,svytable) S3method(model_parameters,systemfit) @@ -536,6 +538,7 @@ S3method(p_value,speedlm) S3method(p_value,stanreg) S3method(p_value,summary.lm) S3method(p_value,survreg) +S3method(p_value,svy2lme) S3method(p_value,svyglm) S3method(p_value,svyglm.nb) S3method(p_value,svyglm.zip) @@ -872,6 +875,7 @@ S3method(standard_error,sem) S3method(standard_error,stanreg) S3method(standard_error,summary.lm) S3method(standard_error,survreg) +S3method(standard_error,svy2lme) S3method(standard_error,svyglm) S3method(standard_error,svyglm.nb) S3method(standard_error,svyglm.zip) diff --git a/NEWS.md b/NEWS.md index 6dbe22731..782e4b541 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +# parameters 0.22.1 + +## New supported models + +* Support for `svy2lme` models from package *svylme*. + # parameters 0.22.0 ## Breaking changes diff --git a/R/extract_random_variances.R b/R/extract_random_variances.R index 620aabfa0..9b042a79e 100644 --- a/R/extract_random_variances.R +++ b/R/extract_random_variances.R @@ -108,6 +108,42 @@ .extract_random_variances.MixMod <- .extract_random_variances.glmmTMB +# svy2lme ------------------------ + +.extract_random_variances.svy2lme <- function(model, ci = 0.95, effects = "random", ...) { + s <- sqrt(as.vector(model$s2)) + stdev <- matrix(s * sqrt(diag(model$L)), ncol = 1) + vcnames <- c(paste0("SD (", model$znames, ")"), "SD (Observations)") + grp_names <- names(model$znames) + if (is.null(grp_names)) { + grp_names <- model$znames + } + + out <- data.frame( + Parameter = vcnames, + Level = NA, + Coefficient = c(as.vector(stdev), s), + SE = NA, + CI_low = NA, + CI_high = NA, + t = NA, + df_error = NA, + p = NA, + Effects = "random", + Group = c(grp_names, "Residual"), + stringsAsFactors = FALSE + ) + + # fix intercept names + out$Parameter <- gsub("(Intercept)", "Intercept", out$Parameter, fixed = TRUE) + + if (effects == "random") { + out[c("t", "df_error", "p")] <- NULL + } + + rownames(out) <- NULL + out +} diff --git a/R/methods_BayesFactor.R b/R/methods_BayesFactor.R index a5fd14625..cb8b85cd8 100644 --- a/R/methods_BayesFactor.R +++ b/R/methods_BayesFactor.R @@ -63,7 +63,6 @@ model_parameters.BFBayesFactor <- function(model, es_type = NULL, include_proportions = FALSE, verbose = TRUE, - effectsize_type = NULL, ...) { insight::check_if_installed("BayesFactor") @@ -84,12 +83,6 @@ model_parameters.BFBayesFactor <- function(model, return(NULL) } - ## TODO: remove deprecation warning later - if (!is.null(effectsize_type)) { - insight::format_warning("Argument `effectsize_type` is deprecated. Use `es_type` instead.") - es_type <- effectsize_type - } - out <- bayestestR::describe_posterior( model, centrality = centrality, @@ -143,12 +136,12 @@ model_parameters.BFBayesFactor <- function(model, # needs {effectsize} to be installed insight::check_if_installed("effectsize") - ## TODO: add back ci-argument, once effectsize >= 0.7.1 is on CRAN. tryCatch( { effsize <- effectsize::effectsize(model, centrality = centrality, dispersion = dispersion, + ci = ci, ci_method = ci_method, rope_ci = rope_ci, type = es_type, diff --git a/R/methods_aov.R b/R/methods_aov.R index c74b3c2a0..0d13590c0 100644 --- a/R/methods_aov.R +++ b/R/methods_aov.R @@ -35,7 +35,6 @@ #' (e.g., `"g"`, `"l"`, `"two"`...). See section *One-Sided CIs* in #' the [effectsize_CIs vignette](https://easystats.github.io/effectsize/). #' @inheritParams model_parameters.default -#' @param effectsize_type Deprecated. Use `es_type` instead. #' @param ... Arguments passed to [`effectsize::effectsize()`]. For example, #' to calculate _partial_ effect sizes types, use `partial = TRUE`. For objects #' of class `htest` or `BFBayesFactor`, `adjust = TRUE` can be used to return @@ -110,18 +109,11 @@ model_parameters.aov <- function(model, drop = NULL, table_wide = FALSE, verbose = TRUE, - effectsize_type = NULL, ...) { # save model object, for later checks original_model <- model object_name <- insight::safe_deparse_symbol(substitute(model)) - ## TODO: remove deprecation warning later - if (!is.null(effectsize_type)) { - insight::format_warning("Argument `effectsize_type` is deprecated. Use `es_type` instead.") - es_type <- effectsize_type - } - if (inherits(model, "aov") && !is.null(type) && type > 1) { if (requireNamespace("car", quietly = TRUE)) { model <- car::Anova(model, type = type) diff --git a/R/methods_htest.R b/R/methods_htest.R index 4e070c61d..4cf38ba3f 100644 --- a/R/methods_htest.R +++ b/R/methods_htest.R @@ -46,14 +46,7 @@ model_parameters.htest <- function(model, bootstrap = FALSE, es_type = NULL, verbose = TRUE, - effectsize_type = NULL, ...) { - ## TODO: remove deprecation warning later - if (!is.null(effectsize_type)) { - insight::format_warning("Argument `effectsize_type` is deprecated. Use `es_type` instead.") - es_type <- effectsize_type - } - if (bootstrap) { insight::format_error("Bootstrapped h-tests are not yet implemented.") } else { diff --git a/R/methods_other.R b/R/methods_other.R index f5110bf59..402f29ab8 100644 --- a/R/methods_other.R +++ b/R/methods_other.R @@ -41,14 +41,7 @@ model_parameters.Gam <- function(model, type = NULL, table_wide = FALSE, verbose = TRUE, - effectsize_type = NULL, ...) { - ## TODO: remove deprecation warning later - if (!is.null(effectsize_type)) { - insight::format_warning("Argument `effectsize_type` is deprecated. Use `es_type` instead.") - es_type <- effectsize_type - } - model_parameters( summary(model)$parametric.anova, es_type = es_type, diff --git a/R/methods_svy2lme.R b/R/methods_svy2lme.R new file mode 100644 index 000000000..90899b91a --- /dev/null +++ b/R/methods_svy2lme.R @@ -0,0 +1,113 @@ +#' @export +model_parameters.svy2lme <- function(model, + ci = 0.95, + effects = "all", + keep = NULL, + drop = NULL, + verbose = TRUE, + include_sigma = FALSE, + ...) { + dots <- list(...) + # which component to return? + effects <- match.arg(effects, choices = c("fixed", "random", "all")) + params <- params_variance <- NULL + + if (effects %in% c("fixed", "all")) { + # Processing + fun_args <- list( + model, + ci = ci, + ci_method = "wald", + standardize = NULL, + p_adjust = NULL, + wb_component = FALSE, + keep_parameters = keep, + drop_parameters = drop, + verbose = verbose, + include_sigma = include_sigma, + summary = FALSE, + vcov = NULL, + vcov_args = NULL + ) + fun_args <- c(fun_args, dots) + params <- do.call(".extract_parameters_mixed", fun_args) + + params$Effects <- "fixed" + } + + att <- attributes(params) + + if (effects %in% c("random", "all")) { + params_variance <- .extract_random_variances( + model, + ci = ci, + effects = effects + ) + } + + # merge random and fixed effects, if necessary + if (!is.null(params) && !is.null(params_variance)) { + params$Level <- NA + params$Group <- "" + params <- params[match(colnames(params_variance), colnames(params))] + } + + params <- rbind(params, params_variance) + # remove empty column + if (!is.null(params$Level) && all(is.na(params$Level))) { + params$Level <- NULL + } + + # due to rbind(), we lose attributes from "extract_parameters()", + # so we add those attributes back here... + if (!is.null(att)) { + attributes(params) <- utils::modifyList(att, attributes(params)) + } + + params <- .add_model_parameters_attributes( + params, + model, + ci = ci, + exponentiate = FALSE, + bootstrap = FALSE, + iterations = 1000, + ci_method = "wald", + p_adjust = NULL, + verbose = verbose, + summary = FALSE, + group_level = FALSE, + wb_component = FALSE, + ... + ) + + attr(params, "object_name") <- insight::safe_deparse_symbol(substitute(model)) + class(params) <- c("parameters_model", "see_parameters_model", class(params)) + + params +} + + +#' @export +standard_error.svy2lme <- function(model, ...) { + .data_frame( + Parameter = .remove_backticks_from_string(colnames(model$Vbeta)), + SE = as.vector(sqrt(diag(model$Vbeta))) + ) +} + + +#' @export +p_value.svy2lme <- function(model, ...) { + stat <- insight::get_statistic(model) + p <- 2 * stats::pnorm(abs(stat$Statistic), lower.tail = FALSE) + .data_frame( + Parameter = stat$Parameter, + p = as.vector(p) + ) +} + + +#' @export +degrees_of_freedom.svy2lme <- function(model, ...) { + Inf +} diff --git a/R/n_clusters_easystats.R b/R/n_clusters_easystats.R index 6a9205a82..486d1c4c9 100644 --- a/R/n_clusters_easystats.R +++ b/R/n_clusters_easystats.R @@ -527,6 +527,6 @@ plot.n_clusters_dbscan <- plot.n_clusters_elbow #' @export plot.n_clusters_hclust <- function(x, ...) { insight::check_if_installed("pvclust") - graphics::plot(attributes(x)$model) - pvclust::pvrect(attributes(x)$model, alpha = attributes(x)$ci, pv = "si") + graphics::plot(attributes(x)[["model"]]) + pvclust::pvrect(attributes(x)[["model"]], alpha = attributes(x)$ci, pv = "si") } diff --git a/R/utils_format.R b/R/utils_format.R index 6c7a98295..7a4956657 100644 --- a/R/utils_format.R +++ b/R/utils_format.R @@ -954,16 +954,16 @@ } # fix column output - if (inherits(attributes(x)$model, c("lavaan", "blavaan")) && "Label" %in% colnames(x)) { + if (inherits(attributes(x)[["model"]], c("lavaan", "blavaan")) && "Label" %in% colnames(x)) { x$From <- ifelse(!nzchar(as.character(x$Label), keepNA = TRUE) | x$Label == x$To, x$From, paste0(x$From, " (", x$Label, ")")) # nolint x$Label <- NULL } - if (inherits(attributes(x)$model, c("lavaan", "blavaan")) && !"Parameter" %in% colnames(x)) { + if (inherits(attributes(x)[["model"]], c("lavaan", "blavaan")) && !"Parameter" %in% colnames(x)) { parameter_column <- colnames(x)[1] } - if (inherits(attributes(x)$model, c("lavaan", "blavaan")) && "Defined" %in% x$Component) { + if (inherits(attributes(x)[["model"]], c("lavaan", "blavaan")) && "Defined" %in% x$Component) { x$From[x$Component == "Defined"] <- "" x$Operator[x$Component == "Defined"] <- "" x$To <- ifelse(x$Component == "Defined", paste0("(", x$To, ")"), x$To) @@ -1175,7 +1175,7 @@ # fix non-equal length of columns final_table <- .fix_nonmatching_columns( final_table, - is_lavaan = inherits(attributes(x)$model, c("lavaan", "blavaan")) + is_lavaan = inherits(attributes(x)[["model"]], c("lavaan", "blavaan")) ) do.call(rbind, final_table) } else { diff --git a/inst/WORDLIST b/inst/WORDLIST index 879650eb9..b48027070 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -318,6 +318,7 @@ strengejacke subclusters subscale subscales +svylme systemfit th tidymodels diff --git a/man/model_parameters.BFBayesFactor.Rd b/man/model_parameters.BFBayesFactor.Rd index b88312181..5bb4ecf8e 100644 --- a/man/model_parameters.BFBayesFactor.Rd +++ b/man/model_parameters.BFBayesFactor.Rd @@ -17,7 +17,6 @@ es_type = NULL, include_proportions = FALSE, verbose = TRUE, - effectsize_type = NULL, ... ) } @@ -70,8 +69,6 @@ information is often redundant.} \item{verbose}{Toggle off warnings.} -\item{effectsize_type}{Deprecated. Use \code{es_type} instead.} - \item{...}{Additional arguments to be passed to or from methods.} } \value{ diff --git a/man/model_parameters.aov.Rd b/man/model_parameters.aov.Rd index 4033ecc41..bd24bb816 100644 --- a/man/model_parameters.aov.Rd +++ b/man/model_parameters.aov.Rd @@ -18,7 +18,6 @@ drop = NULL, table_wide = FALSE, verbose = TRUE, - effectsize_type = NULL, ... ) @@ -97,8 +96,6 @@ be in the same row. Default: \code{FALSE}.} \item{verbose}{Toggle warnings and messages.} -\item{effectsize_type}{Deprecated. Use \code{es_type} instead.} - \item{...}{Arguments passed to \code{\link[effectsize:effectsize]{effectsize::effectsize()}}. For example, to calculate \emph{partial} effect sizes types, use \code{partial = TRUE}. For objects of class \code{htest} or \code{BFBayesFactor}, \code{adjust = TRUE} can be used to return diff --git a/man/model_parameters.cgam.Rd b/man/model_parameters.cgam.Rd index 7b684a368..2cfcad6c6 100644 --- a/man/model_parameters.cgam.Rd +++ b/man/model_parameters.cgam.Rd @@ -39,7 +39,6 @@ type = NULL, table_wide = FALSE, verbose = TRUE, - effectsize_type = NULL, ... ) @@ -160,8 +159,6 @@ ANOVA-tables using \code{car::Anova()} will be returned. (Ignored for \item{table_wide}{Logical that decides whether the ANOVA table should be in wide format, i.e. should the numerator and denominator degrees of freedom be in the same row. Default: \code{FALSE}.} - -\item{effectsize_type}{Deprecated. Use \code{es_type} instead.} } \value{ A data frame of indices related to the model's parameters. diff --git a/man/model_parameters.htest.Rd b/man/model_parameters.htest.Rd index c988acf4d..b0a069c49 100644 --- a/man/model_parameters.htest.Rd +++ b/man/model_parameters.htest.Rd @@ -12,7 +12,6 @@ bootstrap = FALSE, es_type = NULL, verbose = TRUE, - effectsize_type = NULL, ... ) @@ -46,8 +45,6 @@ models, can also be a character vector with multiple effect size names.} \item{verbose}{Toggle warnings and messages.} -\item{effectsize_type}{Deprecated. Use \code{es_type} instead.} - \item{...}{Arguments passed to or from other methods. For instance, when \code{bootstrap = TRUE}, arguments like \code{type} or \code{parallel} are passed down to \code{bootstrap_model()}.} diff --git a/tests/testthat/_snaps/svylme.md b/tests/testthat/_snaps/svylme.md new file mode 100644 index 000000000..b00a05ce8 --- /dev/null +++ b/tests/testthat/_snaps/svylme.md @@ -0,0 +1,26 @@ +# model_parameters svylme + + Code + print(mp) + Output + # Fixed Effects + + Parameter | Coefficient | SE | 95% CI | t | p + -------------------------------------------------------------------- + (Intercept) | -60.98 | 34.48 | [-128.57, 6.61] | -1.77 | 0.077 + ell | 0.92 | 0.26 | [ 0.41, 1.42] | 3.56 | < .001 + mobility | -0.38 | 0.24 | [ -0.85, 0.08] | -1.60 | 0.109 + api99 | 1.10 | 0.03 | [ 1.03, 1.17] | 31.44 | < .001 + + # Random Effects + + Parameter | Coefficient + ----------------------------------- + SD (Intercept: dnum1) | 1.19 + SD (api99: dnum2) | 1.39e-03 + SD (Residual) | 20.00 + Message + + Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed + using a Wald t-distribution approximation. + diff --git a/tests/testthat/test-rstanarm.R b/tests/testthat/test-rstanarm.R index 7cece5df7..706f6eb57 100644 --- a/tests/testthat/test-rstanarm.R +++ b/tests/testthat/test-rstanarm.R @@ -25,7 +25,7 @@ test_that("mp2", { data(pbcLong, package = "rstanarm") pbcLong$ybern <- as.integer(pbcLong$logBili >= mean(pbcLong$logBili)) set.seed(123) - invisible(capture.output( + invisible(capture.output({ model <- rstanarm::stan_mvmer( formula = list( ybern ~ year + (1 | id), @@ -35,7 +35,7 @@ test_that("mp2", { refresh = 0, seed = 123 ) - )) + })) mp <- suppressWarnings(model_parameters(model, centrality = "mean")) s <- summary(model) diff --git a/tests/testthat/test-svylme.R b/tests/testthat/test-svylme.R new file mode 100644 index 000000000..f8a6bc470 --- /dev/null +++ b/tests/testthat/test-svylme.R @@ -0,0 +1,27 @@ +skip_on_cran() +skip_on_os(c("mac", "linux", "solaris")) + +skip_if_not_installed("withr") +skip_if_not_installed("survey") +skip_if_not_installed("lme4") +skip_if_not_installed("svylme") + +withr::with_environment( + new.env(), + test_that("model_parameters svylme", { + data(api, package = "survey") + # two-stage cluster sample + dclus2 <- survey::svydesign( + id = ~ dnum + snum, + fpc = ~ fpc1 + fpc2, + data = apiclus2 + ) + m <- svylme::svy2lme( + api00 ~ ell + mobility + api99 + (1 + api99 | dnum), + design = dclus2, + method = "nested" + ) + mp <- model_parameters(m) + expect_snapshot(print(mp)) + }) +)