Potential integration with papaja

[crsh]

Hi there,

great package! I recently found some time to take it for a spin and find it super useful. In fact, I'm considering to use or support it in my package papaja. papaja exports apa_print(), a function that facilitate reporting results of your analyses in accordance with APA guidelines. It has been my stated goal from the start that I would like to do no computing if possible---just formatting and typesetting. Unfortunately, many analysis objects do not include all desired effect size estimates. Because of this, many methods report only unstandardized effect sizes end up doing some computations.

Some recent internal changes provide an avenue to change this with, potentially, the help of effectsize. I'm looking for a way to let users pass a function that is used internally to calculate the desired effect size and confidence interval and include it in the output. For example,

library("dplyr")

apa_print.afex_aov <- function(x, es = effectsize::omega_squared) {
  aov_table <- apa_print(x, es = NULL)$table
  aov_es_table <- es(x) %>% 

    # The following would be replaced by a tidy-method
    data.frame() %>% 
    rename_with(~"estimate", matches("Eta|Epsilon|Omega")) %>% 
    rename(
      conf.low = CI_low
      , conf.high = CI_high
    ) %>%

    # This is where the typesetting happens
    mutate(estimate = printnum(estimate, digits = 3, gt1 = FALSE)) %>% 
    rowwise() %>% 
    mutate(conf.int = print_confint(conf.low, conf.high, digits = 3, gt1 = FALSE)) %>% 
    select(estimate, conf.int) %>% 
    label_variables(
      estimate = "\\hat\\eta^2_p"
      , conf.int = "90% CI"
    )

  aov_table$estimate <- aov_es_table$estimate
  aov_table$conf.int <- aov_es_table$conf.int

  # Assembly of the strings for reporting follows
  # ...
}

Such es functions should work on the analysis output x and return a data.frame in a standardized format: There should be the columns estimate, conf.low, and conf.high plus additional columns such as term where applicable (we largely follow the conventions defined in broom). Such a standardized format would make it very easy for us to process the results.

From what I can tell, effectsize is very close to meeting our needs. Functions process common analysis outputs and return a predictable data.frame with a proper S3 class. As such, it seems like a great candidate to develop this workflow with.

How could we do it?
To really make this work, it would be ideal for us, if two things could be added to effectsize:

1. As far as I know, broom is starting to request package owners to ship tidy()-methods for their own objects. Hence, it would be great, if effectsize could provide such methods (which essentially would only rename column names and remove the added class tag).
2. For ANOVA, we currently default to generalized eta-squared (without CI). As far as I can tell, effectsize does not support this effect size, yet. If it could be added, we could delegate all ANOVA-related effect size computations to effectsize. It seems like you are planning to add this anyway.
3. For comparisons of linear regression models, we use a bootstrap-method to estimate confidence intervals of delta R-squared. If you would be willing to adopt this method into effectsize I would finally be able to avoid doing any computations within papaja.

Do you feel that these changes are within the scope of effectsize and would you consider PRs that implement those changes?

[DominiqueMakowski]

That seems like a great idea to me ☺️ for 1) afaik it's feasible (@strengejacke we should maybe a generic tidy() method in insight that would be reexported in all easystats packages? 'cause I reckon it will be quite similar) and for 3) maybe it would be a feature that could be better added in performance?

[strengejacke]

for 1) There is parameters::standardize_names(), maybe we could see that it covers all relevant columns in effectsize as well, though I don't know which names are "standard" in broom. Maybe someone can comment on / open a PR? The related code in question is here:

https://github.com/easystats/parameters/blob/ba887ff912c17e7326affc995af075e65fcfda64/R/standardize_names.R#L57

[strengejacke]

3. would fit into effectsize since we already have Add Cohen's f for R2 change #151. Not sure if CIs aren't already supported? Though limited to lm() only at the moment. But indeed @mattansb suggested adding this to performance::compare_performance() as well.

[mattansb]

2. I hope to have generalized eta squared working in the next few days.

would fit into effectsize since we already have #151. Not sure if CIs aren't already supported? Though limited to lm() only at the moment. But indeed @mattansb suggested adding this to performance::compare_performance() as well.

If implemented in performance with the various R2 functions, and some bootstrapping, it could definitely work (the reason effectsize only support lms for now is because its easy to get the F, and dfs for a parameteric CI, but if we're going the bootstrap way...).

for 1) There is parameters::standardize_names(), maybe we could see that it covers all relevant columns in effectsize as well,

All names from effectsize can be found here:
https://github.com/easystats/effectsize/blob/master/R/is_effectsize_name.R

Perhaps standardize_names can be moved to insight to support across the whole easyvrse?

[IndrajeetPatil]

All names from effectsize can be found here:

Should all these be renamed to estimate in broom convention or something like effsize.estimate?
I don't think broom has any convention about effect size names.

[strengejacke]

@crsh

The only thing I'm unsure about is "Effects"

This indicates a model component, i.e. if parameters are fixed or random effects, so this could remain as is.

[crsh]

All names from effectsize can be found here:

Should all these be renamed to estimate in broom convention or something like effsize.estimate?
I don't think broom has any convention about effect size names.

Considering that unstandardized effect size estimates are also effect size estimates, personally, I would opt for estimate.

The only thing I'm unsure about is "Effects"

This indicates a model component, i.e. if parameters are fixed or random effects, so this could remain as is.

Agreed, broom.mixed uses effect to distinguish between fixed and random. I updated my draft above.

[strengejacke]

Are you sure about cn[cn %in% c("df_residual", "df_error")] <- "df.den"?

[mattansb]

@crsh Generalized Eta squared is now implemented in the package for lm, aov, aovlist, afex_aov models via eta_squared(model, generalized = ...). For the first 3, can pass the observed variables, else all are assumed to be manipulated (for afex_aov they are extracted from the model).
See example:

https://easystats.github.io/effectsize/articles/anovaES.html#generalized-eta2

[strengejacke]

And all effect size names listed in https://github.com/easystats/effectsize/blob/master/R/is_effectsize_name.R should be renamed to estimate, right?

[strengejacke]

broom::tidy(chisq.test(x))

I'm not sure, does broom::tidy() not have any df-column for its tidiers?

[crsh]

It is my experience with other aspects of broom that the conventions are not always super consistently used across tidiers. As with the one.way example above, here is another where df is used:

tidy(aov(yield ~ block + N * P + K, npk))

# A tibble: 6 x 6
  term         df  sumsq meansq statistic  p.value
  <chr>     <dbl>  <dbl>  <dbl>     <dbl>    <dbl>
1 block         5 343.    68.7      4.39   0.0130 
2 N             1 189.   189.      12.1    0.00368
3 P             1   8.40   8.40     0.537  0.476  
4 K             1  95.2   95.2      6.09   0.0271 
5 N:P           1  21.3   21.3      1.36   0.263  
6 Residuals    14 219.    15.6     NA     NA      

I think df.den is actually den.df

Yes, you are right, that was a typo. I updated the draft above.

In general the glossary lists df, num.df, and den.df specifically for tidy(). df, num.df, den.df, null.df, and df.residual can be used for glance(). So, I think, it depends on the model context what the most appropriate would be.

[strengejacke]

I thought about adding a tidy() method, but we can't define the generics (unless we override them), and I'm not willing to import generics in insight. So should we just stick to standardize_names()?

[strengejacke]

Looking at the initial issue:

1. is solved (via insight::standardize_names() resp. insight::tidy())
2. is solved (Potential integration with papaja #156 (comment))
3. still missing. I'm still not sure that a delta-R2 is a high-priority issue for compare_performance(), in particular since we allow different types of models. How should we address this? If all models are of same type, add a delta-R2 column? This would really only be R2(delta2,1) = R2(model2) - R2(model1) and R2(delta3,2) = R2(model3) - R2(model2) etc.?

[DominiqueMakowski]

How should we address this? If all models are of same type, add a delta-R2 column?

Let's open an issue in performance. But it would be cool if it worked across models type, like if model 1 (the reference) includes any type of R2, then compute the difference between it and all the other models that include any type of R2. Sure it will have limitations that some of the R2 will not be of the same nature and all, but in the end all the R2 aim to quantify the same dimension (whether it can be called prop of variance explained or not), and people do already compare different R2s anyway so...

[strengejacke]

Let's open an issue in performance.

easystats/performance#28

[strengejacke]

but in the end all the R2 aim to quantify the same dimension (whether it can be called prop of variance explained or not), and people do already compare different R2s anyway so...

That's why I try to convince my colleagues not to report pseudo-R2 ;-)

[mattansb]

This is also why whatever function ends up doing this should validate that all models have the same outcome + number of observations (similar to what anova() does).

[crsh]

Should we ignore possible warnings about namespace conflicts and define a tidy() generic?

Hmm, wouldn't this potentially interfere with the method dispatch for other objects when broom is loaded prior to insight? I'm not sure what your reasons are for not importing generics, but if this is a hard decision I think the only clean way to add a method to tidy() is to offer it to broom, no?

Looking at the initial issue:

1. is solved (via insight::standardize_names() resp. insight::tidy())

2. is solved (#156 (comment))

3. still missing. I'm still not sure that a delta-R2 is a high-priority issue for compare_performance(), in particular since we allow different types of models. How should we address this? If all models are of same type, add a delta-R2 column? This would really only be R2(delta2,1) = R2(model2) - R2(model1) and R2(delta3,2) = R2(model3) - R2(model2) etc.?

Outsourcing our delta-R2 method is not a hard precondition for this integration. Rather, I thought it may fit well with one of the packages of your ecosystem. I think if there is a nice way to add the tidy()-method, I would feel good about starting the work on our end.

[strengejacke]

Hmm, wouldn't this potentially interfere with the method dispatch for other objects when broom is loaded prior to insight?

Not sure, since both have the same generic (i.e. tidy(x, ...) ).

but if this is a hard decision I think the only clean way to add a method to tidy() is to offer it to broom, no?

Yeah, we can do that, I'll open a PR.

[strengejacke]

Anyone willing to open a PR at the broom-repo, to add a tidy() method that call insight::standardize_names()? Not sure what kind of PR template is required...

@crsh Would it be enough to provide standardize_names() instead of tidy(), to avoid namespace conflicts?

[crsh]

Anyone willing to open a PR at the broom-repo

I can do it. I can check back with you before opening the PR, so you can ensure that we are on the same page.

Would it be enough to provide standardize_names() instead of tidy(), to avoid namespace conflicts?

Hmm, I don't think broom will accept a tidy method under a different name (i.e. standardize_names()).

Just to make sure I understand correctly, you have in mind something like the following?

tidy <- function(x, ...) {
  insight::standardize_names(x, ...)
}

[strengejacke]

Would it be enough to provide standardize_names() instead of tidy(), to avoid namespace conflicts?

No, this was addressed to you for the papaja package, not for the broom PR ;-)

Just to make sure I understand correctly, you have in mind something like the following?

Yeah, but with style = "broom":

tidy <- function(x, ...) {
  insight::standardize_names(x, style = "broom", ...)
}

The current working draft is here:
https://github.com/easystats/insight/blob/master/WIP/tidy_easystats.R

You may remove the note in the docs.

[crsh]

No, this was addressed to you for the papaja package, not for the broom PR ;-)

Ah, gotcha. 😅 I would like to keep the interface general and enable users to pass other functions to compute effect sizes as well. For this reason, I need a predictable format of the return object for the subsequent internal workflow. Calling broom::tidy() on the return object would be a widely applicable solution. A more lightweight solution, however, could be to first check if the output already is in the desired format and in that case skip tidying.

I could imagine something along these lines:

aov_es_table <- es(x)
aov_es_cols <- colnames(aov_es_table)
if(
  !is.data.frame(aov_es_table) |
  !"estimate" %in% aov_es_cols |
  !"conf.low" %in% aov_es_cols |
  !"conf.high" %in% aov_es_cols
) {
  aov_es_table <- broom::tidy(aov_es_table)
}

If, for example, eta_sqaured() could be made to return a data.frame with broom-style standardized names, this could be an alternative to adding a tidy-method to broom. Would this work better for you?

[strengejacke]

For easystats packages / functions, you could just check the class-attribute:

if(inherits(c(aov_es_table, "effectsize_table", "parameters_model", "parameters_distribution"))) {
  aov_es_table <- insight::standardize_names(aov_es_table, style = "broom")
} else {
  aov_es_table <- broom::tidy(aov_es_table)
}

[crsh]

After some more discussion with @mariusbarth, we think that this approach probably strikes a good balance between streamlining our internals, ease of use for the user, and flexibility. So, I think, with the broom-style output of insight::standardize_names, we should be good to go. Thanks everyone!

[crsh]

Hi again. We have started to explore our implementation and another issue has come up. Currently, we support testing and reporting results for intercepts in ANOVA models. It seems that effectsize does not support this and explicitly excludes the intercept term:

effectsize/R/eta_squared.R

Line 892 in b4f8408

model <- model[model$Parameter != "(Intercept)", , drop = FALSE]

Given that several packages allow users to test the intercept term (e.g. afex), we felt that users who do so should also be able to report the results including effect sizes. We were thus wondering whether you would accept a PR that allows users to optionally receive effect size estimates (and CI) for this term. Any thoughts on this?

[mattansb]

This should actually be a small (I think) fix.
The reason for dropping it was that it messed up the computation of total effect sizes (non partial eta square, etc.).
Unless you already started working on this PR, I can take a jab at it (:

[crsh]

Great, glad to hear it. There's nothing in the works, yet, so feel free. 👍

[mattansb]

So I've added an optional argument include_intercept (that defaults to FALSE):

library(effectsize)

model <- aov(mpg ~ factor(am) * factor(cyl), data = mtcars)


eta_squared(car::Anova(model, type = 3), include_intercept = TRUE)
#> Parameter              | Eta2 (partial) |       90% CI
#> ------------------------------------------------------
#> (Intercept)            |           0.87 | [0.78, 0.91]
#> factor(am)             |           0.20 | [0.02, 0.41]
#> factor(cyl)            |           0.41 | [0.15, 0.58]
#> factor(am):factor(cyl) |           0.10 | [0.00, 0.27]

however, it does depend on the anova table having an intercept to begin with,
so these don’t get one:

eta_squared(car::Anova(model, type = 2), include_intercept = TRUE)
#> Parameter              | Eta2 (partial) |       90% CI
#> ------------------------------------------------------
#> factor(am)             |           0.13 | [0.00, 0.34]
#> factor(cyl)            |           0.66 | [0.45, 0.77]
#> factor(am):factor(cyl) |           0.10 | [0.00, 0.27]
eta_squared(model, include_intercept = TRUE)
#> Parameter              | Eta2 (partial) |       90% CI
#> ------------------------------------------------------
#> factor(am)             |           0.63 | [0.42, 0.75]
#> factor(cyl)            |           0.66 | [0.45, 0.77]
#> factor(am):factor(cyl) |           0.10 | [0.00, 0.27]

Neither does this:

data(md_12.1, package = "afex")
m <- aov(rt ~ angle * noise + Error(id / (angle * noise)),
         data = md_12.1)
eta_squared(m, include_intercept = TRUE)
#> Group          |   Parameter | Eta2 (partial) |       90% CI
#> ------------------------------------------------------------
#> id:angle       |       angle |           0.82 | [0.66, 0.89]
#> id:noise       |       noise |           0.79 | [0.49, 0.89]
#> id:angle:noise | angle:noise |           0.83 | [0.69, 0.90]

But afex_aov objects genereally will[*]:

a <- afex::aov_ez("id", "rt", md_12.1, within = c("angle", "noise")) # same model!
eta_squared(a, include_intercept = TRUE)
#> Parameter   | Eta2 (partial) |       90% CI
#> -------------------------------------------
#> (Intercept) |           0.99 | [0.96, 0.99]
#> angle       |           0.82 | [0.66, 0.89]
#> noise       |           0.79 | [0.49, 0.89]
#> angle:noise |           0.83 | [0.69, 0.90]

[*] how effect sizes are extracted from afex_aov models is a bit convoluted. So you can get the intercept in the following cases:

• When the design is fully between - always possible.
• When the desing is mixed / within:
  • For partial eta squared.
  • For generalized eta squared.
  • For partial epsilon squared.
  • For partial Cohen's f / f squared.

Let me know if this is good (or good enough?).

[mattansb]

@crsh I think I've addressed all of your points. Feel free to chime back in with any thoughts.