Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

emmeans objects for certain models with a different type than the default fail in apa_print. #592

Open
TheDom42 opened this issue Sep 17, 2024 · 3 comments

Comments

@TheDom42
Copy link

TheDom42 commented Sep 17, 2024

Describe the bug
emmeans objects for certain models with a different type than the default fail in apa_print. Examples that I have identified are type = unlink and type = response (https://cran.r-project.org/web/packages/emmeans/vignettes/transformations.html#tranlink) for Gamma models. There are probably more.

To Reproduce

library("lme4")
#> Loading required package: Matrix
library("broom")
library("emmeans")
#> Welcome to emmeans.
#> Caution: You lose important information if you filter this package's results.
#> See '? untidy'
library("papaja")
#> Loading required package: tinylabels
glmm <- glm(sqrt(breaks) ~ wool*tension, family = Gamma, data = warpbreaks)


#### regular link works
emm.glmm <- emmeans(glmm,  ~ tension)
#> NOTE: Results may be misleading due to involvement in interactions
tidy(emm.glmm, conf.int = TRUE)
#> # A tibble: 3 × 8
#>   tension estimate std.error    df conf.low conf.high statistic  p.value
#>   <chr>      <dbl>     <dbl> <dbl>    <dbl>     <dbl>     <dbl>    <dbl>
#> 1 L          0.172   0.00747    48    0.157     0.187      23.0 1.41e-27
#> 2 M          0.198   0.00856    48    0.181     0.215      23.1 1.13e-27
#> 3 H          0.219   0.00948    48    0.200     0.238      23.1 1.16e-27

apa_print(emm.glmm)
#> $estimate
#> $estimate$L
#> [1] "$1/M = 0.17$, 95\\% CI $[0.16, 0.19]$"
#> 
#> $estimate$M
#> [1] "$1/M = 0.20$, 95\\% CI $[0.18, 0.22]$"
#> 
#> $estimate$H
#> [1] "$1/M = 0.22$, 95\\% CI $[0.20, 0.24]$"
#> 
#> 
#> $statistic
#> $statistic$L
#> [1] "$t(48) = 23.01$, $p < .001$"
#> 
#> $statistic$M
#> [1] "$t(48) = 23.13$, $p < .001$"
#> 
#> $statistic$H
#> [1] "$t(48) = 23.11$, $p < .001$"
#> 
#> 
#> $full_result
#> $full_result$L
#> [1] "$1/M = 0.17$, 95\\% CI $[0.16, 0.19]$, $t(48) = 23.01$, $p < .001$"
#> 
#> $full_result$M
#> [1] "$1/M = 0.20$, 95\\% CI $[0.18, 0.22]$, $t(48) = 23.13$, $p < .001$"
#> 
#> $full_result$H
#> [1] "$1/M = 0.22$, 95\\% CI $[0.20, 0.24]$, $t(48) = 23.11$, $p < .001$"
#> 
#> 
#> $table
#> A data.frame with 6 labelled columns:
#> 
#>   tension estimate     conf.int statistic df p.value
#> L       L     0.17 [0.16, 0.19]     23.01 48  < .001
#> M       M     0.20 [0.18, 0.22]     23.13 48  < .001
#> H       H     0.22 [0.20, 0.24]     23.11 48  < .001
#> 
#> tension  : tension 
#> estimate : $1/M$ 
#> conf.int : 95\\% CI 
#> statistic: $t$ 
#> df       : $\\mathit{df}$ 
#> p.value  : $p$ 
#> attr(,"class")
#> [1] "apa_results" "list"


#### with unlink or other links, this fails
emm.glmm_unlk <- emmeans(glmm,  ~ tension, type = "unlink")
#> NOTE: Results may be misleading due to involvement in interactions
tidy(emm.glmm_unlk, conf.int = TRUE)
#> # A tibble: 3 × 9
#>   tension response std.error    df conf.low conf.high  null statistic  p.value
#>   <chr>      <dbl>     <dbl> <dbl>    <dbl>     <dbl> <dbl>     <dbl>    <dbl>
#> 1 L           5.82     0.253    48     5.35      6.38   Inf      23.0 1.41e-27
#> 2 M           5.05     0.218    48     4.65      5.53   Inf      23.1 1.13e-27
#> 3 H           4.56     0.197    48     4.20      5.00   Inf      23.1 1.16e-27

apa_print(emm.glmm_unlk)
#> Error in apa_num.default(tidy_x$estimate, ...): The parameter 'x' is NULL. Please provide a value for 'x'

Created on 2024-09-17 with reprex v2.1.1

Expected behavior
apa_results should have been returned with the column of response as the estimate.

Additional context
This links back to

# Hot fix while we wait on a merge of this PR: https://github.com/tidymodels/broom/pull/1047
tidy_x <- rename_column(
tidy_x
, c("prob", "rate", "ratio", "odds.ratio", "asymp.LCL", "asymp.UCL", "z.ratio")
, c("estimate", "estimate", "estimate", "estimate", "conf.low", "conf.high", "statistic")
)
and in turn the non-standard naming in broom due to the different link type.

Because the tidy_x$estimate cannot be found, this fails.

@crsh
Copy link
Owner

crsh commented Sep 17, 2024

Hi Dominik,

Thanks for reaching out and pointing me to the origin of the problem. It seems to me this is something that should ideally be addressed in broom, but we could implement a hot fix in papaja until this happens? Would you be willing to give a PR a shot? And could you open an issue on the broom repository about this to see whether a fix on their end is in the cards?

If not, we can take care of this, but it will almost certainly take longer. ;)

Cheers,
Frederik

@TheDom42
Copy link
Author

The earliest I would have time to try a PR would be October.

However, before trying the fix, I was also thinking about something else:

est_name <- switch(
link
, "logit" = "\\mathrm{logit}(p)"
, "probit" = "\\Phi^{-1}(p)"
, "cauchit" = "\\mathrm{cauchit}(p)"
, "cloglog" = "\\mathrm{cloglog}(p)"
, "log" = "\\log(M)"
, "sqrt" = "\\sqrt{M}"
, "log.o.r." = "\\log(\\mathit{OR})"
, "exp" = "\\exp(M)"
, "inverse" = "1/M"
, "M"

With this code, the label for the estimate column is decided. It defaults to $M$. Currently, I find it quite hard to imagine a good logic to decide, which label to use for a "generic" response column. This could lead to misinterpretations if the response/DV is not actually expressed as a mean.

@crsh
Copy link
Owner

crsh commented Sep 18, 2024

I agree that there are cases in which M is not ideal. We could think about a more generic default, such as \hat{y}. It is always possible to overwrite anything we output here:

apa_print(emm.glmm, est_name = "\\hat{y}")

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants