diff --git a/DESCRIPTION b/DESCRIPTION index 1bb0eca..1c41858 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: esci Title: Estimation Statistics with Confidence Intervals -Version: 1.0.2 +Version: 1.0.3 Authors@R: person("Robert", "Calin-Jageman", email = "rcalinjageman@dom.edu", role = c("aut", "cre", "cph")) Description: A collection of functions and 'jamovi' module for the estimation approach to inferential statistics, the approach which emphasizes effect sizes, interval estimates, and meta-analysis. Nearly all functions are based on 'statpsych' and 'metafor'. This package is still under active development, and breaking changes are likely, especially with the plot and hypothesis test functions. Data sets are included for all examples from Cumming & Calin-Jageman (2024) . @@ -28,7 +28,7 @@ Imports: Rdpack, stringr, mathjaxr -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 Roxygen: list(markdown = TRUE) RdMacros: Rdpack, mathjaxr Suggests: diff --git a/NAMESPACE b/NAMESPACE index b9e9075..7d20240 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -71,17 +71,13 @@ importFrom(rlang,enquo) importFrom(statpsych,ci.lc.mean.bs) importFrom(statpsych,ci.lc.median.bs) importFrom(statpsych,ci.lc.prop.bs) -importFrom(statpsych,ci.mean1) -importFrom(statpsych,ci.median1) importFrom(statpsych,ci.oddsratio) importFrom(statpsych,ci.phi) -importFrom(statpsych,ci.prop1) importFrom(statpsych,ci.ratio.mean.ps) importFrom(statpsych,ci.ratio.mean2) importFrom(statpsych,ci.ratio.median.ps) importFrom(statpsych,ci.ratio.median2) importFrom(statpsych,ci.stdmean.ps) -importFrom(statpsych,ci.stdmean1) importFrom(stats,aggregate) importFrom(stats,anova) importFrom(stats,complete.cases) @@ -90,10 +86,12 @@ importFrom(stats,cor) importFrom(stats,line) importFrom(stats,median) importFrom(stats,na.omit) +importFrom(stats,pbinom) importFrom(stats,pnorm) importFrom(stats,predict) importFrom(stats,pt) importFrom(stats,qnorm) +importFrom(stats,qt) importFrom(stats,quantile) importFrom(stats,rnorm) importFrom(stats,sd) diff --git a/NEWS.md b/NEWS.md index 5c21e5b..bb67c91 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,12 @@ +esci version 1.0.3 (Release data: July 2024) +=========== + +Changes: + +* Tweaks to deal with changes in statpsych 1.6. No substantive updates. + + + esci version 1.0.2 (Release data: March 2024) =========== diff --git a/R/CI_smd_ind_contrast.R b/R/CI_smd_ind_contrast.R index 6520f05..cf5eac0 100644 --- a/R/CI_smd_ind_contrast.R +++ b/R/CI_smd_ind_contrast.R @@ -2,11 +2,11 @@ #' contrast #' #' -#' @description `CI_smd_ind_contrast` returns the point estimate +#' @description \loadmathjax `CI_smd_ind_contrast` returns the point estimate #' and confidence interval for a standardized mean difference (smd aka Cohen's #' *d* aka Hedges *g*). A standardized mean difference is a difference in means standardized #' to a standard deviation: -#' d = psi/s +#' \mjdeqn{d = \frac{ \psi }{s}}{d = psi/s} #' #' #' @param means A vector of 2 or more means @@ -88,16 +88,16 @@ #' When equal variance is assumed, the standardized mean difference is #' d_s, defined in Kline, p. 196: #' -#' d_s = psi / sd_s +#' \mjdeqn{ d_s = \frac{ \psi }{ sd_{pooled} } }{ d_s = psi / sd_s} #' #' #' where psi is defined in Kline, equation 7.8: #' -#' psi = sum(contrasts*means) +#' \mjdeqn{ \psi =\sum_{i=1}^{a}c_iM_i }{psi = sum(contrasts*means)} #' #' #' and where sd_pooled is defined in Kline, equation 3.11 -#' sqrt(sum(variances*dfs) / sum(dfs)) +#' \mjdeqn{sd_{pooled} = { \frac{ \sum_{i=1}^{a} (n_i -1) s_i^2 } { \sum_{i=1}^{a} (n_i-1) } }}{sqrt(sum(variances*dfs) / sum(dfs))} #' #' #' The CI for d_s is derived from lambda-prime transformation from Lecoutre, @@ -115,12 +115,12 @@ #' assumed, the standardized mean difference is d_avg, defined in Bonett, #' equation 6: #' -#' d_avg = psi / sd_avg +#' \mjdeqn{ d_{avg} = \frac{ \psi }{ sd_{avg} }}{ d_avg = psi / sd_avg} #' #' Where sd_avg is the square root of the average of the group variances, as #' given in Bonett, explanation of equation 6: #' -#' sqrt(mean(variances)) +#' \mjdeqn{sd_{avg} = \sqrt{ \frac{ \sum_{i=1}^{a} s_i^2 }{ a } }}{sqrt(mean(variances))} #' #' #' #### If only 2 groups #### diff --git a/R/esci-package.R b/R/esci-package.R index e6f8b97..61a6dfc 100644 --- a/R/esci-package.R +++ b/R/esci-package.R @@ -13,16 +13,12 @@ #' @importFrom statpsych ci.lc.mean.bs #' @importFrom statpsych ci.lc.median.bs #' @importFrom statpsych ci.lc.prop.bs -#' @importFrom statpsych ci.mean1 -#' @importFrom statpsych ci.median1 #' @importFrom statpsych ci.oddsratio #' @importFrom statpsych ci.phi -#' @importFrom statpsych ci.prop1 #' @importFrom statpsych ci.ratio.mean2 #' @importFrom statpsych ci.ratio.mean.ps #' @importFrom statpsych ci.ratio.median2 #' @importFrom statpsych ci.ratio.median.ps -#' @importFrom statpsych ci.stdmean1 #' @importFrom statpsych ci.stdmean.ps #' @importFrom stats aggregate #' @importFrom stats anova @@ -35,8 +31,10 @@ #' @importFrom stats qnorm #' @importFrom stats predict #' @importFrom stats pt +#' @importFrom stats pbinom #' @importFrom stats pnorm #' @importFrom stats qnorm +#' @importFrom stats qt #' @importFrom stats quantile #' @importFrom stats rnorm #' @importFrom stats sd diff --git a/R/estimate_magnitude.R b/R/estimate_magnitude.R index e59293b..d8daf4d 100644 --- a/R/estimate_magnitude.R +++ b/R/estimate_magnitude.R @@ -16,9 +16,11 @@ #' If you want to compare your sample to a known value or reference, then #' use [esci::estimate_mdiff_one()]. #' -#' The estimated mean is from [statpsych::ci.mean1()]. +#' The estimated mean is from [statpsych::ci.mean1()] (renamed ci.mean as of +#' statpsych 1.6). #' -#' The estimated median is from [statpsych::ci.median1()] +#' The estimated median is from [statpsych::ci.median1()] (renamed ci.median +#' as of statpsych 1.6) #' #' #' @param data For raw data - A data frame or tibble @@ -466,7 +468,7 @@ estimate_magnitude.vector <- function( estimate$properties$data_source <- NULL # 2 alpha CI for median - mdn_2a <- statpsych::ci.median1( + mdn_2a <- wrapper_ci.median( alpha = (1 - conf_level)*2, y =outcome_variable[!is.na(outcome_variable)] ) diff --git a/R/estimate_mdiff_one.R b/R/estimate_mdiff_one.R index 09832f3..96cc486 100644 --- a/R/estimate_mdiff_one.R +++ b/R/estimate_mdiff_one.R @@ -16,11 +16,13 @@ #' it with [esci::plot_mdiff()] and you can test hypotheses with #' [esci::test_mdiff()]. #' -#' The estimated mean differences are from [statpsych::ci.mean1()]. +#' The estimated mean differences are from [statpsych::ci.mean1()] (renamed +#' ci.mean as of statpsych 1.6). #' #' The estimated SMDs are from [esci::CI_smd_one()]. #' -#' The estimated median differences are from [statpsych::ci.median1()] +#' The estimated median differences are from [statpsych::ci.median1()] (renamed +#' ci.median as of statpsych 1.6) #' #' #' @param data For raw data - a data frame or tibble diff --git a/R/estimate_pdiff_one.R b/R/estimate_pdiff_one.R index a376b49..30b2026 100644 --- a/R/estimate_pdiff_one.R +++ b/R/estimate_pdiff_one.R @@ -15,7 +15,8 @@ #' it with [esci::plot_pdiff()] and you can test hypotheses with #' [esci::test_pdiff()]. #' -#' The estimated proportion differences are from [statpsych::ci.prop1()]. +#' The estimated proportion differences are from [statpsych::ci.prop1()] (renamed +#' ci.prop as of statpsych 1.6). #' #' #' @param data For raw data - a dataframe or tibble diff --git a/R/estimate_proportion.R b/R/estimate_proportion.R index 5afedd4..dab176a 100644 --- a/R/estimate_proportion.R +++ b/R/estimate_proportion.R @@ -14,7 +14,8 @@ #' If you want to compare your estimate to a known value or reference, then #' use [esci::estimate_pdiff_one()]. #' -#' The estimated proportions are from [statpsych::ci.prop1()]. +#' The estimated proportions are from [statpsych::ci.prop1()] (renamed +#' ci.prop as of statpsych 1.6). #' #' #' @param data For raw data - a data frame or tibble diff --git a/R/overview.R b/R/overview.R index 06e198a..789d75e 100644 --- a/R/overview.R +++ b/R/overview.R @@ -285,7 +285,7 @@ overview.base <- function( # Analysis ---------------------------------------- if (n_means == 1) { res <- as.data.frame( - statpsych::ci.mean1( + wrapper_ci.mean ( alpha = 1 - conf_level, m = overview_table$mean, sd = overview_table$sd, @@ -725,7 +725,7 @@ wrap_ci_median1 <- function( if(length(na.omit(x)) < 2) return(c(median(x), NA, NA, NA)) - res <- statpsych::ci.median1( + res <- wrapper_ci.median( alpha = 1 - conf_level, y = if(na.rm) x[!is.na(x)] else x ) @@ -743,7 +743,7 @@ wrap_ci_mean1 <- function( if(length(x) < 2) return(c(mean(x), NA, NA, NA)) - res <- statpsych::ci.mean1( + res <- wrapper_ci.mean( alpha = 1 - conf_level, m = mean(x, na.rm = na.rm), sd = sd(x, na.rm = na.rm), diff --git a/R/overview_nominal.R b/R/overview_nominal.R index 2c6d4e6..236c2a2 100644 --- a/R/overview_nominal.R +++ b/R/overview_nominal.R @@ -244,7 +244,7 @@ overview_nominal.base <- function( res_ta <- res } else { res <- as.data.frame( - statpsych::ci.prop1( + wrapper_ci.prop( alpha = 1 - conf_level, f = overview_table$cases[x], n = overview_table$n[x] @@ -252,7 +252,7 @@ overview_nominal.base <- function( ) res_ta <- as.data.frame( - statpsych::ci.prop1( + wrapper_ci.prop( alpha = (1 - conf_level)*2, f = overview_table$cases[x], n = overview_table$n[x] diff --git a/R/plot_aesthetics.R b/R/plot_aesthetics.R index 678ecb3..62b5e61 100644 --- a/R/plot_aesthetics.R +++ b/R/plot_aesthetics.R @@ -197,6 +197,11 @@ esci_plot_mdiff_aesthetics <- function( esci_plot_simple_aesthetics <- function(myplot, use_ggdist = TRUE) { # Customize plot ------------------------------- # No legend + + ggplot_version <- as.numeric(gsub("\\.", "", utils::packageVersion("ggplot2"))) + + + myplot <- myplot + ggplot2::theme(legend.position = "none") # Points @@ -214,24 +219,43 @@ esci_plot_simple_aesthetics <- function(myplot, use_ggdist = TRUE) { if (use_ggdist) { myplot <- myplot + ggplot2::discrete_scale( - c("size", "point_size"), - "point_size_d", - function(n) return(c("raw" = 1, "summary" = 3)) + aesthetics = c("size", "point_size"), + scale_name = if (ggplot_version >= 350) NULL else "point_size_d", + palette = function(n) return(c("raw" = 1, "summary" = 3)) ) } else { + + if (ggplot_version >= 350) { + myplot <- myplot + ggplot2::discrete_scale( + aesthetics = c("size"), + palette = function(n) return(c("raw" = 1, "summary" = 1)) + ) + + } else { + myplot <- myplot + ggplot2::discrete_scale( + aesthetics = c("size"), + scale_name = "point_size_d", + palette = function(n) return(c("raw" = 1, "summary" = 1)) + ) + + } + + } + + if (ggplot_version >= 350) { myplot <- myplot + ggplot2::discrete_scale( - c("size"), - "point_size_d", - function(n) return(c("raw" = 1, "summary" = 1)) + aesthetics = c("alpha", "point_alpha"), + palette = function(n) return(c("raw" = 0.8, "summary" = 1)) + ) + } else { + myplot <- myplot + ggplot2::discrete_scale( + c("alpha", "point_alpha"), + "point_alpha_d", + function(n) return(c("raw" = 0.8, "summary" = 1)) ) } - myplot <- myplot + ggplot2::discrete_scale( - c("alpha", "point_alpha"), - "point_alpha_d", - function(n) return(c("raw" = 0.8, "summary" = 1)) - ) # Error bars myplot <- myplot + ggplot2::scale_linetype_manual( @@ -241,16 +265,31 @@ esci_plot_simple_aesthetics <- function(myplot, use_ggdist = TRUE) { values = c("summary" = "black"), aesthetics = "interval_color" ) - myplot <- myplot + ggplot2::discrete_scale( - "interval_alpha", - "interval_alpha_d", - function(n) return(c("summary" = 1)) - ) - myplot <- myplot + ggplot2::discrete_scale( - "interval_size", - "interval_size_d", - function(n) return(c("summary" = 3)) - ) + + if (ggplot_version >= 350) { + myplot <- myplot + ggplot2::discrete_scale( + aesthetics = "interval_alpha", + palette = function(n) return(c("summary" = 1)) + ) + myplot <- myplot + ggplot2::discrete_scale( + aesthetics = "interval_size", + palette = function(n) return(c("summary" = 3)) + ) + + } else { + + myplot <- myplot + ggplot2::discrete_scale( + aesthetics = "interval_alpha", + scale_name = "interval_alpha_d", + palette = function(n) return(c("summary" = 1)) + ) + myplot <- myplot + ggplot2::discrete_scale( + aesthetics = "interval_size", + scale_name = "interval_size_d", + palette = function(n) return(c("summary" = 3)) + ) + + } # Slab @@ -258,11 +297,21 @@ esci_plot_simple_aesthetics <- function(myplot, use_ggdist = TRUE) { values = c("summary" = "gray"), aesthetics = "slab_fill" ) - myplot <- myplot + ggplot2::discrete_scale( - "slab_alpha", - "slab_alpha_d", - function(n) return(c("summary" = 1)) - ) + + if (ggplot_version >= 350) { + myplot <- myplot + ggplot2::discrete_scale( + aesthetics = "slab_alpha", + palette = function(n) return(c("summary" = 1)) + ) + + } else { + myplot <- myplot + ggplot2::discrete_scale( + aesthetics = "slab_alpha", + scale_name = "slab_alpha_d", + palette = function(n) return(c("summary" = 1)) + ) + + } return(myplot) } diff --git a/R/statpsych_apply.R b/R/statpsych_apply.R index 5057e9a..338f7ff 100644 --- a/R/statpsych_apply.R +++ b/R/statpsych_apply.R @@ -222,7 +222,7 @@ apply_ci_prop1 <- function( ) { res <- as.data.frame( - statpsych::ci.prop1( + wrapper_ci.prop( alpha = 1 - conf_level, f = myrow[["cases"]], n = myrow[["N"]] @@ -279,7 +279,7 @@ apply_ci_mean1 <- function( ) { res <- as.data.frame( - statpsych::ci.mean1( + wrapper_ci.mean( alpha = 1 - conf_level, m = myrow[["mean"]] - reference_mean, sd = myrow[["sd"]], diff --git a/R/statpsych_wrapper.R b/R/statpsych_wrapper.R index 45077d1..306e7fc 100644 --- a/R/statpsych_wrapper.R +++ b/R/statpsych_wrapper.R @@ -9,7 +9,7 @@ wrapper_ci.stdmean1 <- function( # Result from statpsych ------------------------------- res <- as.data.frame( - statpsych::ci.stdmean1( + stolen_ci.stdmean( alpha = 1 - conf_level, m = comparison_mean, sd = comparison_sd, @@ -405,7 +405,7 @@ wrapper_ci.median.ps <- function( ) { res <- as.data.frame( - statpsych::ci.median1( + wrapper_ci.median( alpha = 1 - conf_level, y = comparison_measure ) @@ -414,7 +414,7 @@ wrapper_ci.median.ps <- function( res <- rbind( res, as.data.frame( - statpsych::ci.median1( + wrapper_ci.median( alpha = 1 - conf_level, y = reference_measure ) @@ -628,3 +628,99 @@ wrapper_ci.cor2 <- function( return(res) } + + +wrapper_ci.mean <- function(alpha, m, sd, n) { + + # return(statpsych::ci.mean(alpha, m, sd, n)) + + df <- n - 1 + tcrit <- qt(1 - alpha/2, df) + se <- sd/sqrt(n) + ll <- m - tcrit*se + ul <- m + tcrit*se + out <- t(c(m, se, ll, ul)) + colnames(out) <- c("Estimate", "SE", "LL", "UL") + rownames(out) <- "" + return(out) + +} + +wrapper_ci.median <- function(alpha, y) { + # return(statpsych::ci.median(alpha, y)) + + n <- length(y) + y <- sort(y) + z <- qnorm(1 - alpha/2) + median <- median(y) + c1 <- round((n - z*sqrt(n))/2) + if (c1 < 1) {c1 = 1} + ll <- y[c1] + ul <- y[n - c1 + 1] + a <- round(n/2 - sqrt(n)) + if (a < 1) {a = 1} + ll1 <- y[a] + ul1 <- y[n - a + 1] + p <- pbinom(a - 1, size = n, prob = .5) + z0 <- qnorm(1 - p) + se <- (ul1 - ll1)/(2*z0) + out <- t(c(median, se, ll, ul)) + colnames(out) <- c("Estimate", "SE", "LL", "UL") + rownames(out) <- "" + return(out) + + +} + +wrapper_ci.prop <- function(alpha, f, n) { + + # return( + # statpsych::ci.prop( + # alpha = alpha, + # f = f, + # n = n + # ) + # ) + + if (f > n) {stop("f cannot be greater than n")} + z <- qnorm(1 - alpha/2) + p.mle <- f/n + se.mle <- sqrt(p.mle*(1 - p.mle)/n) + b1 <- 2*n*p.mle + z^2 + b2 <- 2*(n + z^2) + LL.wil <- (b1 - 1 - z*sqrt(z^2 - 2 - 1/n + 4*p.mle*(n*(1 - p.mle) + 1)))/b2 + UL.wil <- (b1 + 1 + z*sqrt(z^2 + 2 - 1/n + 4*p.mle*(n*(1 - p.mle) - 1)))/b2 + if (p.mle == 0) {LL.wil = 0} + if (p.mle == 1) {UL.wil = 1} + p.adj <- (f + 2)/(n + 4) + se.adj <- sqrt(p.adj*(1 - p.adj)/(n + 4)) + LL.adj <- p.adj - z*se.adj + UL.adj <- p.adj + z*se.adj + if (LL.adj < 0) {LL.adj = 0} + if (UL.adj > 1) {UL.adj = 1} + out1 <- t(c(p.adj, se.adj, LL.adj, UL.adj)) + out2 <- t(c(p.mle, se.mle, LL.wil, UL.wil)) + out <- rbind(out1, out2) + colnames(out) <- c("Estimate", "SE", "LL", "UL") + rownames(out) <- c("Adjusted Wald", "Wilson with cc") + return(out) + + +} + + +stolen_ci.stdmean <- function(alpha, m, sd, n, h) { + z <- qnorm(1 - alpha/2) + df <- n - 1 + adj <- 1 - 3/(4*df - 1) + est <- (m - h)/sd + estu <- adj*est + se <- sqrt(est^2/(2*df) + 1/df) + ll <- est - z*se + ul <- est + z*se + out <- t(c(est, estu, se, ll, ul)) + colnames(out) <- c("Estimate", "adj Estimate", "SE", "LL", "UL") + rownames(out) <- "" + return(out) +} + diff --git a/README.md b/README.md index 17b8d3f..dea126b 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# esci +# esci esci logo esci provides student-friendly tools for estimation statistics: @@ -90,4 +90,4 @@ plot_mdiff(estimate) ``` - +Example difference plot diff --git a/_pkgdown.yml b/_pkgdown.yml index 8829b59..c78dbdd 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -2,6 +2,7 @@ url: https://rcalinjageman.github.io/esci template: bootstrap: 5 + math-rendering: mathjax development: mode: auto diff --git a/docs/404.html b/docs/404.html index 949ad5a..bf59d05 100644 --- a/docs/404.html +++ b/docs/404.html @@ -8,59 +8,46 @@ Page not found (404) • esci - - - - + + + Skip to contents - -