Fixes easystats/report#245

# Load effectsize
devtools::load_all(".../forks/effectsize")

# "Right" method
y <- t.test(mtcars$mpg ~ mtcars$vs)
cohens_d(y)
#> Cohen's d |         95% CI
#> --------------------------
#> -1.70     | [-2.55, -0.82]
#> 
#> - Estimated using un-pooled SD.

# Approximation method
x <- t.test(mpg ~ vs, data = mtcars)
cohens_d(x)
#> Warning: Unable to retrieve data from htest object. Returning an approximate
#>   effect size using t_to_d().
#> d     |         95% CI
#> ----------------------
#> -1.96 | [-2.94, -0.95]

# New method with data specified
cohens_d(x, data = mtcars)
#> Cohen's d |         95% CI
#> --------------------------
#> -1.70     | [-2.55, -0.82]
#> 
#> - Estimated using un-pooled SD.

# Comparison
cohens_d(y) == cohens_d(x, data = mtcars)
#>                                                 Cohens_d   CI CI_low CI_high
#> difference in means between group 0 and group 1     TRUE TRUE   TRUE    TRUE

# Approximation method
report::report_table(x)
#> Warning: Unable to retrieve data from htest object. Returning an approximate
#>   effect size using t_to_d().
#> Welch Two Sample t-test
#> 
#> Parameter | Group | Mean_Group1 | Mean_Group2 | Difference |          95% CI | t(22.72) |      p |     d |          d  CI
#> -------------------------------------------------------------------------------------------------------------------------
#> mpg       |    vs |       16.62 |       24.56 |      -7.94 | [-11.46, -4.42] |    -4.67 | < .001 | -1.96 | [-2.94, -0.95]
#> 
#> Alternative hypothesis: two.sided

# New method with data specified
report::report_table(x, data = mtcars)
#> Welch Two Sample t-test
#> 
#> Parameter | Group | Mean_Group1 | Mean_Group2 | Difference |          95% CI | t(22.72) |      p | Cohen's d |   Cohens_d  CI
#> -----------------------------------------------------------------------------------------------------------------------------
#> mpg       |    vs |       16.62 |       24.56 |      -7.94 | [-11.46, -4.42] |    -4.67 | < .001 |     -1.70 | [-2.55, -0.82]
#> 
#> Alternative hypothesis: two.sided

^{Created on 2022-10-16 with reprex v2.0.2}

Oct 16 '22 23:10 rempsyc

Codecov Report

Attention: Patch coverage is 90.78947% with 7 lines in your changes are missing coverage. Please review.

Project coverage is 90.75%. Comparing base (16d832b) to head (1dc280f). Report is 1 commits behind head on main.

:exclamation: Current head 1dc280f differs from pull request most recent head 21c6ad3. Consider uploading reports for the commit 21c6ad3 to get more accurate results

Files	Patch %	Lines
R/effectsize.htest.R	90.14%	7 Missing :warning:

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #522      +/-   ##
==========================================
+ Coverage   90.71%   90.75%   +0.04%     
==========================================
  Files          57       56       -1     
  Lines        3563     3440     -123     
==========================================
- Hits         3232     3122     -110     
+ Misses        331      318      -13

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

Oct 16 '22 23:10 codecov-commenter

I haven't given a close look yet, but this seems like it does not cover the full scope of possible ~garbage~ htest data.name patterns...

See @strengejacke code here: https://github.com/easystats/insight/blob/main/R/utils_get_data.R

Oct 17 '22 18:10 mattansb

Ok. (When you have time) would it be possible to give me an example of what you mean by "htest data.name patterns"? Do you mean like to support the different htest objects (it wasn't clear to me from the link you provided)? I just added support for cohen's d and t-test for now, but am happy to add support for more if useful.

Oct 17 '22 21:10 rempsyc

[x] chisq.test (no data argument so formula interface not possible)
[x] fisher.test (no data argument so formula interface not possible)
[x] mcnemar.test (no data argument so formula interface not possible)
[x] cor.test (already works)
[x] oneway.test (already works)
[x] t.test (done)
[x] wilcox.test (done)
[x] friedman.test (done)
[x] kruskal.test (done)
[ ] chisq.test_gof (not sure what function to use to test this object type)

Dec 16 '22 23:12 rempsyc

devtools::load_all("D:/github/forks/effectsize")
#> ℹ Loading effectsize

# t.test
x <- t.test(mpg ~ vs, data = mtcars)
effectsize(x)
#> Warning: Unable to retrieve data from htest object. Returning an approximate
#>   effect size using t_to_d().
#> d     |         95% CI
#> ----------------------
#> -1.96 | [-2.94, -0.95]
effectsize(x, data = mtcars)
#> Cohen's d |         95% CI
#> --------------------------
#> -1.70     | [-2.55, -0.82]
#> 
#> - Estimated using un-pooled SD.

# cor.test
# no need to specify the data argument
x <- cor.test(~ qsec + drat, data = mtcars)
effectsize(x)
#> Warning: This 'htest' method is not (yet?) supported.
#>   Returning 'parameters::model_parameters(model)'.
#> Pearson's product-moment correlation
#> 
#> Parameter1 | Parameter2 |    r |        95% CI | t(30) |     p
#> --------------------------------------------------------------
#> qsec       |       drat | 0.09 | [-0.27, 0.43] |  0.50 | 0.620
#> 
#> Alternative hypothesis: true correlation is not equal to 0

# wilcox.test
x <- wilcox.test(mpg ~ vs, data = mtcars)
#> Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
#> compute exact p-value with ties
effectsize(x)
#> Error: Unable to retrieve data from htest object.
#>   Try using 'rb()' directly.
effectsize(x, data = mtcars)
#> Warning: 'y' is numeric but has only 2 unique values.
#>   If this is a grouping variable, convert it to a factor.
#> r (rank biserial) |       95% CI
#> --------------------------------
#> 1.00              | [1.00, 1.00]

# chisq.test
# No data argument, so it is not possible to use the formula interface.

# fisher.test
# No data argument, so it is not possible to use the formula interface.

# chisq.test_gof
# ???

# mcnemar.test
# No data argument, so it is not possible to use the formula interface.

# friedman.test
wb <- aggregate(warpbreaks$breaks, by = list(
  w = warpbreaks$wool, t = warpbreaks$tension), FUN = mean)
wb
#>   w t        x
#> 1 A L 44.55556
#> 2 B L 28.22222
#> 3 A M 24.00000
#> 4 B M 28.77778
#> 5 A H 24.55556
#> 6 B H 18.77778
x <- friedman.test(x ~ w | t, data = wb)
effectsize(x)
#> Error: Unable to retrieve data from htest object.
#>   Try using 'kendalls_w()' directly.
effectsize(x, data = wb)
#> Kendall's W |       95% CI
#> --------------------------
#> 0.11        | [0.11, 1.00]
#> 
#> - One-sided CIs: upper bound fixed at [1.00].

# chisq.test_gof
# no need to specify the data argument

# oneway.test
# no need to specify the data argument
x <- oneway.test(extra ~ group, data = sleep, var.equal = TRUE)
effectsize(x)
#> Eta2 |       95% CI
#> -------------------
#> 0.16 | [0.00, 1.00]
#> 
#> - One-sided CIs: upper bound fixed at [1.00].

# kruskal.test
# Cannot get it to work
airquality2 <- airquality
airquality2$Month <- as.factor(airquality2$Month)
x <- kruskal.test(Ozone ~ Month, data = airquality2)
effectsize(x)
#> Warning in grepl("Kruskal-Wallis", x$method, fixed = TRUE) &&
#> grepl("^list\\(", : 'length(x) = 2 > 1' in coercion to 'logical(1)'
#> Error: Unable to retrieve data from htest object.
#>   Try using 'rank_epsilon_squared()' directly.
effectsize(x, data = airquality2)
#> Warning in grepl("Kruskal-Wallis", x$method, fixed = TRUE) &&
#> grepl("^list\\(", : 'length(x) = 2 > 1' in coercion to 'logical(1)'
#> Warning: Missing values detected. NAs dropped.
#> Epsilon2 (rank) |       95% CI
#> ------------------------------
#> 0.25            | [0.15, 1.00]
#> 
#> - One-sided CIs: upper bound fixed at [1.00].

^{Created on 2022-12-16 with reprex v2.0.2}

Dec 16 '22 23:12 rempsyc

Hi! I'm performing a paired t-test and a signed-rank Wilcoxon, I still get the very same warning. What should I do?

Jan 31 '23 01:01 Larissa-Cury

Hey @Larissa-Cury ,

Here is a toy example of your error:

library(effectsize)

wt <- wilcox.test(mpg ~ am, data = mtcars)
effectsize(wt)
#> Error: Unable to retrieve data from htest object.

You have two options:

Call the {effectsize} effect size function directly.

rank_biserial(mpg ~ am, data = mtcars)
#> r (rank biserial) |         95% CI
#> ----------------------------------
#> -0.66             | [-0.84, -0.36]

Use the vector input (and not the formula input) to the *.test() function.

wt2 <- wilcox.test(mtcars$mpg[mtcars$am == 0], mtcars$mpg[mtcars$am == 1])
effectsize(wt2)
#> r (rank biserial) |         95% CI
#> ----------------------------------
#> -0.66             | [-0.84, -0.36]

Jan 31 '23 07:01 mattansb

@rempsyc chisq.test doesn't need any data - it is the only htest that actually contains the data 🥳

To merge this, we would need a lot of unit tests that show that this feature can be used as intended. Here are two examples that fail because of the unrecoverable data - can they be made to work? (given here for t.test but applicable to anything really)

Example 1: advance formula interface
Example 2: vector interface, but data is removed
Example 3: test run with an index (similar to https://github.com/easystats/report/issues/295)

library(effectsize)
library(testthat)

tt1 <- t.test(mpg ~ I(am + cyl == 4), data = mtcars)
dd1 <- cohens_d(mpg ~ I(am + cyl == 4), data = mtcars, pooled_sd = FALSE)

dat <- mtcars
tt2 <- t.test(dat$mpg[dat$am==1], dat$mpg[dat$am==0])
dd2 <- cohens_d(dat$mpg[dat$am==1], dat$mpg[dat$am==0], pooled_sd = FALSE)
rm("dat")

col.y <- "mpg"
tt3 <- t.test(mtcars[[col.y]])
dd3 <- cohens_d(mtcars[[col.y]])
rm("col.y")

expect_equal(effectsize(tt1, data = mtcars)[[1]], dd1[[1]])
#> Error in `[.data.frame`(dots$data, , vars) : undefined columns selected

expect_equal(effectsize(tt2, data = mtcars)[[1]], dd2[[1]])
#> Error in `[.data.frame`(dots$data, , vars) : undefined columns selected

expect_equal(effectsize(tt3, data = mtcars)[[1]], dd3[[1]])
#> Error in `[.data.frame`(dots$data, , vars) : undefined columns selected

Use internal functions to not repeat yourself in code
If some examples cannot be made to work:
- If this feature really worth the trouble just for a formula interface (that sucks)? If so...
- Add informative errors to catch these "extreme" cases.

Jul 03 '23 06:07 mattansb

Thanks, that's very useful and specific! I'll look into implementing these suggestions :)

Jul 03 '23 14:07 rempsyc

[x] unit tests
[x] internal functions to avoid repeating code
[x] informative messages

Doesn’t work well with advanced formula interfaces for now but I return NULL for data if it can’t find it, so it defaults to the same as not providing the argument (plus a message about it). Working on trying to fix the formula stuff with regex right now but not sure if that’s the right approach…

Also, to fix the lints, I was going to change all the .effectsize_t.test to .effectsize_t_test, etc. (name style should match snake_case). But do you even want that?

Jul 03 '23 18:07 rempsyc

Don't bother with the lints now. (Also see me other comments.)

Jul 04 '23 05:07 mattansb

hehe oops...

Perhaps it would be best to focus on just the simple formula interface?

Here is some code to reconstruct the input formula and obtain the data:

library(effectsize)

t-test

## Two 

model <- t.test(mpg ~ am, data = mtcars)
data.name <- model$data.name

data.name <- gsub("by", "~", data.name, fixed = TRUE)
(form <- as.formula(data.name))
#> mpg ~ am

mf <- model.frame(form, data = mtcars)
mf[[2]] <- factor(mf[[2]])
head(mf)
#>                    mpg am
#> Mazda RX4         21.0  1
#> Mazda RX4 Wag     21.0  1
#> Datsun 710        22.8  1
#> Hornet 4 Drive    21.4  0
#> Hornet Sportabout 18.7  0
#> Valiant           18.1  0

## One

model <- t.test(mpg ~ 1, data = mtcars)
data.name <- model$data.name

(form <- as.formula(paste0(data.name, "~1")))
#> mpg ~ 1

mf <- model.frame(form, data = mtcars)
head(mf)
#>                    mpg
#> Mazda RX4         21.0
#> Mazda RX4 Wag     21.0
#> Datsun 710        22.8
#> Hornet 4 Drive    21.4
#> Hornet Sportabout 18.7
#> Valiant           18.1


## Paired

model <- t.test(Pair(mpg, hp) ~ 1, data = mtcars)
data.name <- model$data.name

(form <- as.formula(paste0(data.name, "~1")))
#> Pair(mpg, hp) ~ 1

mf <- model.frame(form, data = mtcars)
head(mf)
#>                   Pair(mpg, hp).x Pair(mpg, hp).y
#> Mazda RX4                    21.0           110.0
#> Mazda RX4 Wag                21.0           110.0
#> Datsun 710                   22.8            93.0
#> Hornet 4 Drive               21.4           110.0
#> Hornet Sportabout            18.7           175.0
#> Valiant                      18.1           105.0

One Way

model <- oneway.test(mpg ~ gear, data = mtcars)
data.name <- model$data.name

data.name <- gsub("and", "~", data.name, fixed = TRUE)
(form <- as.formula(data.name))
#> mpg ~ gear

mf <- model.frame(form, data = mtcars)
mf[[2]] <- factor(mf[[2]])
head(mf)
#>                    mpg gear
#> Mazda RX4         21.0    4
#> Mazda RX4 Wag     21.0    4
#> Datsun 710        22.8    4
#> Hornet 4 Drive    21.4    3
#> Hornet Sportabout 18.7    3
#> Valiant           18.1    3

Wilcox test

## rank sum

model <- wilcox.test(mpg ~ am, data = mtcars)
#> Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
#> compute exact p-value with ties
data.name <- model$data.name

data.name <- gsub("by", "~", data.name, fixed = TRUE)
(form <- as.formula(data.name))
#> mpg ~ am

mf <- model.frame(form, data = mtcars)
mf[[2]] <- factor(mf[[2]])
head(mf)
#>                    mpg am
#> Mazda RX4         21.0  1
#> Mazda RX4 Wag     21.0  1
#> Datsun 710        22.8  1
#> Hornet 4 Drive    21.4  0
#> Hornet Sportabout 18.7  0
#> Valiant           18.1  0

## signed rank (1)

model <- wilcox.test(mpg ~ 1, data = mtcars)
#> Warning in wilcox.test.default(x = respVar, ...): cannot compute exact p-value
#> with ties
data.name <- model$data.name

(form <- as.formula(paste0(data.name, "~1")))
#> mpg ~ 1

mf <- model.frame(form, data = mtcars)
head(mf)
#>                    mpg
#> Mazda RX4         21.0
#> Mazda RX4 Wag     21.0
#> Datsun 710        22.8
#> Hornet 4 Drive    21.4
#> Hornet Sportabout 18.7
#> Valiant           18.1


## signed rank (2)

model <- wilcox.test(Pair(mpg, hp) ~ 1, data = mtcars)
#> Warning in wilcox.test.default(x = respVar[, 1L], y = respVar[, 2L], paired =
#> TRUE, : cannot compute exact p-value with ties
data.name <- model$data.name

(form <- as.formula(paste0(data.name, "~1")))
#> Pair(mpg, hp) ~ 1

mf <- model.frame(form, data = mtcars)
head(mf)
#>                   Pair(mpg, hp).x Pair(mpg, hp).y
#> Mazda RX4                    21.0           110.0
#> Mazda RX4 Wag                21.0           110.0
#> Datsun 710                   22.8            93.0
#> Hornet 4 Drive               21.4           110.0
#> Hornet Sportabout            18.7           175.0
#> Valiant                      18.1           105.0

Kruskal-Wallis Rank Sum Test

model <- kruskal.test(Ozone ~ Month, data = airquality)
data.name <- model$data.name

data.name <- gsub("by", "~", data.name, fixed = TRUE)
(form <- as.formula(data.name))
#> Ozone ~ Month

mf <- model.frame(form, data = airquality)
mf[[2]] <- factor(mf[[2]])
head(mf)
#>   Ozone Month
#> 1    41     5
#> 2    36     5
#> 3    12     5
#> 4    18     5
#> 6    28     5
#> 7    23     5

Friedman Rank Sum Test

wb <- aggregate(warpbreaks$breaks,
                by = list(w = warpbreaks$wool,
                          t = warpbreaks$tension),
                FUN = mean)

model <- friedman.test(x ~ w | t, data = wb)
data.name <- model$data.name

data.name <- sub("and", "~", data.name)
data.name <- sub("and", "+", data.name)
(form <- as.formula(data.name))
#> x ~ w + t

mf <- model.frame(form, data = wb)
mf[[2]] <- factor(mf[[2]])
mf[[3]] <- factor(mf[[3]])
head(mf)
#>          x w t
#> 1 44.55556 A L
#> 2 28.22222 B L
#> 3 24.00000 A M
#> 4 28.77778 B M
#> 5 24.55556 A H
#> 6 18.77778 B H

^{Created on 2023-07-04 with reprex v2.0.2}

However....

Note that all formula interfaces to htest also have optional subset and na.action arguments that don't get saved and would also have to be passed (to model.frame).

Jul 06 '23 10:07 mattansb

~~Most of~~ All these cases now work:

t-test

## Two 
library(effectsize)
model <- t.test(mpg ~ am, data = mtcars)
effectsize(model, data = mtcars)
#> Cohen's d |         95% CI
#> --------------------------
#> -1.41     | [-2.26, -0.53]
#> 
#> - Estimated using un-pooled SD.

## One
model <- t.test(mpg ~ 1, data = mtcars)
effectsize(model, data = mtcars)
#> Cohen's d |       95% CI
#> ------------------------
#> 3.33      | [2.43, 4.22]

## Paired
model <- t.test(Pair(mpg, hp) ~ 1, data = mtcars)
effectsize(model, data = mtcars)
#> Cohen's d |       95% CI
#> ------------------------
#> 1.04      | [0.74, 1.34]

One Way

model <- oneway.test(mpg ~ gear, data = mtcars)
effectsize(model, data = mtcars)
#> `var.equal = FALSE` - effect size is an approximation.
#> Eta2 |       95% CI
#> -------------------
#> 0.70 | [0.31, 1.00]
#> 
#> - One-sided CIs: upper bound fixed at [1.00].

Wilcox test

## rank sum
model <- wilcox.test(mpg ~ am, data = mtcars)
#> Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
#> compute exact p-value with ties
effectsize(model, data = mtcars)
#> r (rank biserial) |         95% CI
#> ----------------------------------
#> -0.66             | [-0.84, -0.36]

## signed rank (1)
model <- wilcox.test(mpg ~ 1, data = mtcars)
#> Warning in wilcox.test.default(x = respVar, ...): cannot compute exact p-value
#> with ties
effectsize(model, data = mtcars)
#> r (rank biserial) |       95% CI
#> --------------------------------
#> 1.00              | [1.00, 1.00]

## signed rank (2)
model <- wilcox.test(Pair(mpg, hp) ~ 1, data = mtcars)
#> Warning in wilcox.test.default(x = respVar[, 1L], y = respVar[, 2L], paired =
#> TRUE, : cannot compute exact p-value with ties
effectsize(model, data = mtcars)
#> r (rank biserial) |       95% CI
#> --------------------------------
#> 1.00              | [1.00, 1.00]

Kruskal-Wallis Rank Sum Test

model <- kruskal.test(Ozone ~ Month, data = airquality)
effectsize(model, data = airquality)
#> Epsilon2 (rank) |       95% CI
#> ------------------------------
#> 0.25            | [0.14, 1.00]
#> 
#> - One-sided CIs: upper bound fixed at [1.00].

Friedman Rank Sum Test

wb <- aggregate(warpbreaks$breaks,
                by = list(w = warpbreaks$wool,
                          t = warpbreaks$tension),
                FUN = mean)

model <- friedman.test(x ~ w | t, data = wb)
effectsize(model, data = wb)
#> Kendall's W |       95% CI
#> --------------------------
#> 0.11        | [0.11, 1.00]
#> 
#> - One-sided CIs: upper bound fixed at [1.00].

Edge cases

library(testthat)
#> Warning: package 'testthat' was built under R version 4.3.1

tt1 <- t.test(mpg ~ I(am + cyl == 4), data = mtcars)
dd1 <- cohens_d(mpg ~ I(am + cyl == 4), data = mtcars, pooled_sd = FALSE)

dat <- mtcars
tt2 <- t.test(dat$mpg[dat$am==1], dat$mpg[dat$am==0])
dd2 <- cohens_d(dat$mpg[dat$am==1], dat$mpg[dat$am==0], pooled_sd = FALSE)
rm("dat")

col.y <- "mpg"
tt3 <- t.test(mtcars[[col.y]])
dd3 <- cohens_d(mtcars[[col.y]])
rm("col.y")

expect_equal(effectsize(tt1, data = mtcars)[[1]], dd1[[1]])

expect_equal(effectsize(tt2, data = mtcars)[[1]], dd2[[1]])

expect_equal(effectsize(tt3, data = mtcars)[[1]], dd3[[1]])

~~This last one will never work I think because the col_y object has been removed, which was specifying which column of mtcars to use, whereas the data argument can only specify the larger data set, not the specific column to extract. Therefore in this case I believe it is acceptable to return the warning about returning an approximate effect size.~~

Fixed.

^{Created on 2023-07-06 with reprex v2.0.2}

I will address the subset and na.action next.

Jul 06 '23 17:07 rempsyc

Note that sometimes the data is not easily recoverable through stats::model.frame:

wb <- aggregate(warpbreaks$breaks, by = list(
  w = warpbreaks$wool, t = warpbreaks$tension
), FUN = mean)

form <- stats::as.formula(x ~ w | t)
form
#> x ~ w | t

# result
stats::model.frame(form, data = wb)
#> Warning in Ops.factor(w, t): '|' not meaningful for factors
#> [1] x     w | t
#> <0 rows> (or 0-length row.names)

wb2 <- aggregate(warpbreaks$breaks, by = list(
  w = as.character(warpbreaks$wool), t = as.character(warpbreaks$tension)
), FUN = mean)

# operations are possible only for numeric, logical or complex types
stats::model.frame(form, data = wb2)
#> Error in w | t: operations are possible only for numeric, logical or complex types

^{Created on 2023-07-06 with reprex v2.0.2}

In this case, we currently use the following (previous) workaround:

data <- dots$data[, vars_split]

However, because we don't use stats::model.frame, we can't pass arguments na.action and subset, so would need to (eventually) account for them manually, or better yet find a way to make the above example work with stats::model.frame.

Jul 06 '23 18:07 rempsyc

I think that's just an issue with | not being defined in formula syntax in base R, only in lme4 and similar. What is that model supposed to be doing?

Jul 06 '23 23:07 bwiernik

Thanks Brenton. This is for a Friedman Rank Sum Test, and the example below is taken from ?friedman.test

wb <- aggregate(warpbreaks$breaks,
                by = list(w = warpbreaks$wool,
                          t = warpbreaks$tension),
                FUN = mean)

model <- friedman.test(x ~ w | t, data = wb)
effectsize(model, data = wb)
#> Kendall's W |       95% CI
#> --------------------------
#> 0.11        | [0.11, 1.00]
#> 
#> - One-sided CIs: upper bound fixed at [1.00].

Any workaround you could think of?

Jul 07 '23 01:07 rempsyc

You're working too hard! We already have functions that clean up data and formula

https://github.com/easystats/effectsize/blob/main/R/utils_validate_input_data.R

Jul 07 '23 03:07 mattansb

You're working too hard! We already have functions that clean up data and formula

Ah damn, should have started with that 😂

Jul 07 '23 14:07 rempsyc

You're working too hard! We already have functions that clean up data and formula

Ok, so am I right that the correct internal function to use with that Friedman test example would be effectsize:::.resolve_formula?

And am using it correctly here?

wb <- aggregate(warpbreaks$breaks,
                by = list(w = warpbreaks$wool,
                          t = warpbreaks$tension),
                FUN = mean)

form <- stats::as.formula(x ~ w | t)
form
#> x ~ w | t

effectsize:::.resolve_formula(formula = form, data = wb)
#> Warning in Ops.factor(w, t): '|' not meaningful for factors
#>          x w | t
#> 1 44.55556    NA
#> 2 28.22222    NA
#> 3 24.00000    NA
#> 4 28.77778    NA
#> 5 24.55556    NA
#> 6 18.77778    NA

^{Created on 2023-07-07 with reprex v2.0.2}

Seems like w and t are not detected correctly. Is it because I should pass an extra argument to .resolve_formula?

Since .resolve_formula internally calls stats::model.frame, I suppose there is no way around that issue, as demonstrated earlier?

https://github.com/easystats/effectsize/blob/9298110c467e2e0de88bccef1cb44ca997673d5d/R/utils_validate_input_data.R#L305

Jul 07 '23 14:07 rempsyc

@rempsyc I sincerely ask you - is this worth it? For this to work - even if just for very standard formula interface - the user will need to supply data, subset, na.action (and even then it might not work due to scoping issues) and at that point... they should just run the relevant effect-size function directly.

Here is a fully working example for a two sample t.test. For this to work you will literally need to re-write about 60% of the effectsize.htest.R file (which we would have to bring the coverage there to 100% first to make sure we don't break anything). And even then it will only maybe work for some people sometime.


get_two_sample_cohens_d_from_formula_t.test <- function(model, ...) {
  # Recover parameters:
  is_paired <- grepl("One", model$method)
  mu <- as.vector(model$null.value)
  pooled_sd <- !grepl("One", model$method) && !grepl("Welch", model$method)
  alt <- model$alternative
  conf.level <- attr(model$conf.int,"conf.level")
  
  # Recover data
  data.name <- model$data.name
  data.name <- gsub("by", "~", data.name, fixed = TRUE)
  form <- as.formula(data.name)
  
  
  data_for_test <- effectsize:::.get_data_2_samples(x = form, 
                                                    paired = is_paired,
                                                    ...)
  effectsize::cohens_d(data_for_test$x, 
                       data_for_test$y, 
                       mu = mu,
                       alternative = alt,
                       ci = conf.level,
                       paired = is_paired,
                       pooled_sd = pooled_sd)
}

# Set up data
some_data <- mtcars
some_data$mpg[1] <- NA

  
tt <- t.test(mpg ~ am, 
             data = some_data,
             
             alternative = "less", 
             mu = 1,
             var.equal = TRUE,
             
             subset = cyl == 4, 
             na.action = na.omit)

d1 <- get_two_sample_cohens_d_from_formula_t.test(tt, 
                                                  data = some_data,    # USER MUST PASS THIS
                                                  subset = cyl == 4,   # USER MUST PASS THIS
                                                  na.action = na.omit) # USER MUST PASS THIS
d1
#> Cohen's d |        95% CI
#> -------------------------
#> -1.54     | [-Inf, -0.24]
#> 
#> - Deviation from a difference of 1.
#> - Estimated using pooled SD.
#> - One-sided CIs: lower bound fixed at [-Inf].


# Manually:
d2 <- effectsize::cohens_d(mpg ~ am, 
                           data = some_data,
                           alternative = "less", 
                           mu = 1,
                           pooled_sd = TRUE,
                           subset = cyl == 4, 
                           na.action = na.omit)

all.equal(d1, d2)
#> [1] TRUE

^{Created on 2023-07-08 with reprex v2.0.2}

Here is the most concise way a user would do this manually:

args <- alist(mpg ~ am, 
              data = some_data,
              alternative = "less", 
              mu = 1,
              pooled_sd = TRUE,
              var.equal = TRUE,
              subset = cyl == 4, 
              na.action = na.omit)

do.call(args, what = t.test)
#> 
#> 	Two Sample t-test
#> 
#> data:  mpg by am
#> t = -2.2727, df = 9, p-value = 0.02457
#> alternative hypothesis: true difference in means between group 0 and group 1 is less than 1
#> 95 percent confidence interval:
#>        -Inf -0.1944731
#> sample estimates:
#> mean in group 0 mean in group 1 
#>          22.900          28.075
#>          

do.call(args, what = effectsize::cohens_d)
#> Cohen's d |        95% CI
#> -------------------------
#> -1.54     | [-Inf, -0.24]
#> 
#> - Deviation from a difference of 1.
#> - Estimated using pooled SD.
#> - One-sided CIs: lower bound fixed at [-Inf].

We cannot achieve such simplicity given the current architecture of the htest class. Arguably, your time would be better spent trying to convince the R developers to make htest class better! Maybe have a call slot? Maybe save the input data (like chisq.test() does).

Jul 08 '23 18:07 mattansb

For this to work you will literally need to re-write about 60% of the effectsize.htest.R file

Are you sure? Below I reproduce your example. Do you think it needs something more?

library(effectsize)

some_data <- mtcars
some_data$mpg[1] <- NA

args <- alist(
  mpg ~ am,
  data = some_data,
  alternative = "less", 
  mu = 1,
  pooled_sd = TRUE,
  var.equal = TRUE,
  subset = cyl == 4, 
  na.action = na.omit
)

tt <- do.call(args, what = t.test)

d1 <- do.call(args, what = cohens_d)

d2 <- do.call(c(alist(tt), args[-1]), what = effectsize)

all.equal(d1$Cohens_d, d2$Cohens_d)
#> [1] TRUE

^{Created on 2023-08-19 with reprex v2.0.2}

we would have to bring the coverage there to 100% first to make sure we don't break anything

I'm happy to start on that first, with your approval. Worst case scenario it's not wasted even if we don't end up merging this PR.

For this to work - even if just for very standard formula interface - the user will need to supply data, subset, na.action

True, but in my experience these arguments are not really used, but even so, we support them if the user really wants them, so I think that's ok.

is this worth it?

I think so, with the above example and all new tests I included, I think it adds value for the user. But maybe there is more that I am missing and hopefully you can point it out to me :P

Aug 20 '23 02:08 rempsyc

Without looking at the code, but looking at the tests, this looks pretty good. Unfortunately I won't really have time to delve into this summer - but I think maybe we can have this done for a fall release?

Aug 21 '23 12:08 mattansb

Mattan, there are still some (new) lints, but I can address them with your approval (since it involves changing a lot of variable names that collide with base R functions)

Apr 28 '24 12:04 rempsyc

I still just don't like this as a solution - it requires the user to do so much work that they might as well just re-run the code to not use a formula (all of these tests run in a fraction of a second), no?

Does this fix anything upstream (e.g. in report)?

If you really think this is advantageous, go a head and make the fixes and merge (:

Apr 30 '24 20:04 mattansb

effectsize
effectsize copied to clipboard

Adds data argument to effectsize.htest/cohens_d

Codecov Report

t-test

One Way

Wilcox test

Kruskal-Wallis Rank Sum Test

Friedman Rank Sum Test

However....

t-test

One Way

Wilcox test

Kruskal-Wallis Rank Sum Test

Friedman Rank Sum Test

Edge cases

effectsize effectsize copied to clipboard

Adds data argument to effectsize.htest/cohens_d

Codecov Report

t-test

One Way

Wilcox test

Kruskal-Wallis Rank Sum Test

Friedman Rank Sum Test

However....

t-test

One Way

Wilcox test

Kruskal-Wallis Rank Sum Test

Friedman Rank Sum Test

Edge cases

effectsize
effectsize copied to clipboard