Skip to content

Commit 806d308

Browse files
strengejackestrengejacke
authored
mismatch between {emmeans}\s and model_parameters()'s adjusted p-values (#1194)
* mismatch between `{emmeans}`\s and `model_parameters()`'s adjusted p-values Fixes #1103 * fix, tests * fix * fix tests * fix * fix test * fix test * suppress warnings from summary * typo
1 parent 28b7d74 commit 806d308

File tree

9 files changed

+274
-158
lines changed

9 files changed

+274
-158
lines changed

‎DESCRIPTION‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
Type: Package
22
Package: parameters
33
Title: Processing of Model Parameters
4-
Version: 0.28.3.3
4+
Version: 0.28.3.4
55
Authors@R:
66
c(person(given = "Daniel",
77
family = "Lüdecke",

‎NEWS.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,8 @@
1717
* Fixed issue where wrong (non-robust) standard errors were calculated for
1818
`coxph` and `svycoxph` objects.
1919

20+
* Fixed issues with Tukey-p-value adjustment for *emmeans* objects.
21+
2022
# parameters 0.28.3
2123

2224
* fixed bug in `standardize_info(<fixest>)` that was preventing

‎R/format_p_adjust.R‎

Lines changed: 71 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,8 @@
1414
format_p_adjust <- function(method) {
1515
method <- tolower(method)
1616

17-
switch(method,
17+
switch(
18+
method,
1819
holm = "Holm (1979)",
1920
hochberg = "Hochberg (1988)",
2021
hommel = "Hommel (1988)",
@@ -44,7 +45,7 @@ format_p_adjust <- function(method) {
4445
# for interaction terms, e.g. for "by" argument in emmeans
4546
# pairwise comparison, we have to adjust the rank resp. the
4647
# number of estimates in a comparison family
47-
rank_adjust <- .p_adjust_rank(model, params)
48+
rank_adjust <- .p_adjust_rank(model, params, tolower(p_adjust))
4849

4950
# only proceed if valid argument-value
5051
if (tolower(p_adjust) %in% tolower(all_methods)) {
@@ -68,13 +69,19 @@ format_p_adjust <- function(method) {
6869
} else if (tolower(p_adjust) == "sidak") {
6970
# sidak adjustment
7071
params$p <- 1 - (1 - params$p)^(nrow(params) / rank_adjust)
71-
} else if (tolower(p_adjust) == "sup-t") {
72+
} else if (tolower(p_adjust) == "sup-t") {
7273
# sup-t adjustment
7374
params <- .p_adjust_supt(model, params)
7475
}
7576

76-
if (isTRUE(all(old_p_vals == params$p)) && !identical(p_adjust, "none") && verbose) {
77-
insight::format_warning(paste0("Could not apply ", p_adjust, "-adjustment to p-values. Either something went wrong, or the non-adjusted p-values were already very large.")) # nolint
77+
if (
78+
isTRUE(all(old_p_vals == params$p)) && !identical(p_adjust, "none") && verbose
79+
) {
80+
insight::format_warning(paste0(
81+
"Could not apply ",
82+
p_adjust,
83+
"-adjustment to p-values. Either something went wrong, or the non-adjusted p-values were already very large."
84+
)) # nolint
7885
}
7986
} else if (verbose) {
8087
insight::format_alert(paste0("`p_adjust` must be one of ", toString(all_methods)))
@@ -86,14 +93,31 @@ format_p_adjust <- function(method) {
8693

8794
# calculate rank adjustment -----
8895

89-
.p_adjust_rank <- function(model, params) {
96+
.p_adjust_rank <- function(model, params, adjust = "tukey") {
9097
tryCatch(
9198
{
9299
correction <- 1
93100
by_vars <- model@misc$by.vars
94-
if (!is.null(by_vars) && by_vars %in% colnames(params)) {
95-
correction <- insight::n_unique(params[[by_vars]])
101+
if (
102+
!is.null(by_vars) && length(by_vars) > 0 && all(by_vars %in% colnames(params))
103+
) {
104+
if (length(by_vars) == 1) {
105+
correction <- insight::n_unique(params[[by_vars]])
106+
} else {
107+
# compute unique combinations of multiple by-variables
108+
by_data <- params[by_vars]
109+
by_groups <- interaction(by_data, drop = TRUE)
110+
correction <- insight::n_unique(by_groups)
111+
}
112+
} else if (identical(adjust, "tukey")) {
113+
# correction <- .safe(prod(vapply(model@model.info$xlev, length, numeric(1))))
114+
correction <- .safe(insight::n_unique(unlist(strsplit(
115+
model@levels$contrast,
116+
" - ",
117+
fixed = TRUE
118+
))))
96119
}
120+
97121
correction
98122
},
99123
error = function(e) {
@@ -106,22 +130,26 @@ format_p_adjust <- function(method) {
106130
# tukey adjustment -----
107131

108132
.p_adjust_tukey <- function(params, stat_column, rank_adjust = 1, verbose = TRUE) {
109-
df_column <- colnames(params)[stats::na.omit(match(c("df", "df_error"), colnames(params)))][1]
133+
df_column <- colnames(params)[stats::na.omit(match(
134+
c("df", "df_error"),
135+
colnames(params)
136+
))][1]
110137
if (!is.na(df_column) && length(stat_column)) {
111138
params$p <- suppressWarnings(stats::ptukey(
112139
sqrt(2) * abs(params[[stat_column]]),
113-
nmeans = nrow(params) / rank_adjust,
140+
nmeans = rank_adjust,
114141
df = params[[df_column]],
115142
lower.tail = FALSE
116143
))
117144
# for specific contrasts, ptukey might fail, and the tukey-adjustement
118145
# could just be simple p-value calculation
119146
if (all(is.na(params$p))) {
120-
params$p <- 2 * stats::pt(
121-
abs(params[[stat_column]]),
122-
df = params[[df_column]],
123-
lower.tail = FALSE
124-
)
147+
params$p <- 2 *
148+
stats::pt(
149+
abs(params[[stat_column]]),
150+
df = params[[df_column]],
151+
lower.tail = FALSE
152+
)
125153
verbose <- FALSE
126154
}
127155
}
@@ -132,7 +160,10 @@ format_p_adjust <- function(method) {
132160
# scheffe adjustment -----
133161

134162
.p_adjust_scheffe <- function(model, params, stat_column, rank_adjust = 1) {
135-
df_column <- colnames(params)[stats::na.omit(match(c("df", "df_error"), colnames(params)))][1]
163+
df_column <- colnames(params)[stats::na.omit(match(
164+
c("df", "df_error"),
165+
colnames(params)
166+
))][1]
136167
if (!is.na(df_column) && length(stat_column)) {
137168
# 1st try
138169
scheffe_ranks <- try(qr(model@linfct)$rank, silent = TRUE)
@@ -146,7 +177,8 @@ format_p_adjust <- function(method) {
146177
scheffe_ranks <- nrow(params)
147178
}
148179
scheffe_ranks <- scheffe_ranks / rank_adjust
149-
params$p <- stats::pf(params[[stat_column]]^2 / scheffe_ranks,
180+
params$p <- stats::pf(
181+
params[[stat_column]]^2 / scheffe_ranks,
150182
df1 = scheffe_ranks,
151183
df2 = params[[df_column]],
152184
lower.tail = FALSE
@@ -182,7 +214,9 @@ format_p_adjust <- function(method) {
182214
# get correlation matrix, based on the covariance matrix
183215
vc <- .safe(stats::cov2cor(insight::get_varcov(model, component = component)))
184216
if (is.null(vc)) {
185-
insight::format_warning("Could not calculate covariance matrix for `sup-t` adjustment.")
217+
insight::format_warning(
218+
"Could not calculate covariance matrix for `sup-t` adjustment."
219+
)
186220
return(params)
187221
}
188222
# get confidence interval level, or set default
@@ -197,18 +231,30 @@ format_p_adjust <- function(method) {
197231
}
198232
# calculate updated confidence interval level, based on simultaenous
199233
# confidence intervals (https://onlinelibrary.wiley.com/doi/10.1002/jae.2656)
200-
crit <- mvtnorm::qmvt(ci_level, df = params[[df_column]][1], tail = "both.tails", corr = vc)$quantile
234+
crit <- mvtnorm::qmvt(
235+
ci_level,
236+
df = params[[df_column]][1],
237+
tail = "both.tails",
238+
corr = vc
239+
)$quantile
201240
# update confidence intervals
202241
params$CI_low <- params$Coefficient - crit * params$SE
203242
params$CI_high <- params$Coefficient + crit * params$SE
204243
# udpate p-values
205244
for (i in 1:nrow(params)) {
206-
params$p[i] <- 1 - mvtnorm::pmvt(
207-
lower = rep(-abs(stats::qt(params$p[i] / 2, df = params[[df_column]][i])), nrow(vc)),
208-
upper = rep(abs(stats::qt(params$p[i] / 2, df = params[[df_column]][i])), nrow(vc)),
209-
corr = vc,
210-
df = params[[df_column]][i]
211-
)
245+
params$p[i] <- 1 -
246+
mvtnorm::pmvt(
247+
lower = rep(
248+
-abs(stats::qt(params$p[i] / 2, df = params[[df_column]][i])),
249+
nrow(vc)
250+
),
251+
upper = rep(
252+
abs(stats::qt(params$p[i] / 2, df = params[[df_column]][i])),
253+
nrow(vc)
254+
),
255+
corr = vc,
256+
df = params[[df_column]][i]
257+
)
212258
}
213259
params
214260
}

‎R/methods_glmmTMB.R‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -240,7 +240,7 @@ model_parameters.glmmTMB <- function(
240240
}
241241

242242
# fix argument, if model has only conditional component
243-
cs <- stats::coef(summary(model))
243+
cs <- suppressWarnings(stats::coef(summary(model)))
244244
has_zeroinf <- modelinfo$is_zero_inflated
245245
has_disp <- is.list(cs) && !is.null(cs$disp)
246246

@@ -616,7 +616,7 @@ standard_error.glmmTMB <- function(
616616
return(se_kenward(model, component = "conditional"))
617617
}
618618

619-
cs <- insight::compact_list(stats::coef(summary(model)))
619+
cs <- suppressWarnings(insight::compact_list(stats::coef(summary(model))))
620620
x <- lapply(names(cs), function(i) {
621621
.data_frame(
622622
Parameter = insight::find_parameters(

‎tests/testthat/_snaps/glmmTMB.md‎

Lines changed: 60 additions & 60 deletions
Original file line numberDiff line numberDiff line change
@@ -7,78 +7,78 @@
77
[2] ""
88
[3] "Parameter | Log-Mean | SE | 95% CI | z | p"
99
[4] "----------------------------------------------------------------"
10-
[5] "child | -1.09 | 0.10 | [-1.28, -0.90] | -11.09 | < .001"
11-
[6] "camper [1] | 0.27 | 0.10 | [ 0.07, 0.47] | 2.70 | 0.007 "
10+
[5] "child | -1.04 | 0.09 | [-1.22, -0.87] | -11.60 | < .001"
11+
[6] "camper [1] | 0.27 | 0.09 | [ 0.08, 0.45] | 2.83 | 0.005 "
1212

1313
# print-model_parameters glmmTMB digits
1414

1515
Code
1616
out[-c(5, 14)]
1717
Output
18-
[1] "# Fixed Effects (Count Model)"
19-
[2] ""
20-
[3] "Parameter | Log-Mean | SE | 95% CI | z | p"
21-
[4] "--------------------------------------------------------------------------"
22-
[5] "child | -1.0875 | 0.0981 | [-1.27967, -0.89528] | -11.0901 | < .001"
23-
[6] "camper [1] | 0.2723 | 0.1009 | [ 0.07461, 0.46999] | 2.6997 | 0.007 "
24-
[7] ""
25-
[8] "# Fixed Effects (Zero-Inflation Component)"
26-
[9] ""
27-
[10] "Parameter | Log-Odds | SE | 95% CI | z | p"
28-
[11] "-----------------------------------------------------------------------"
29-
[12] "(Intercept) | 1.8896 | 0.6642 | [ 0.58780, 3.19147] | 2.8449 | 0.004"
30-
[13] "camper [1] | -0.1701 | 0.3869 | [-0.92836, 0.58822] | -0.4396 | 0.660"
31-
[14] ""
32-
[15] "# Random Effects Variances"
33-
[16] ""
34-
[17] "Parameter | Coefficient | 95% CI"
35-
[18] "---------------------------------------------------------------"
36-
[19] "SD (Intercept: persons) | 3.4056 | [ 1.64567, 7.04777]"
37-
[20] "SD (xb: persons) | 1.2132 | [ 0.59190, 2.48650]"
38-
[21] "Cor (Intercept~xb: persons) | -1.0000 | [-1.00000, 1.00000]"
39-
[22] ""
40-
[23] "# Random Effects (Zero-Inflation Component)"
41-
[24] ""
42-
[25] "Parameter | Coefficient | 95% CI"
43-
[26] "---------------------------------------------------------------"
44-
[27] "SD (Intercept: persons) | 2.7358 | [ 1.16329, 6.43414]"
45-
[28] "SD (zg: persons) | 1.5683 | [ 0.64246, 3.82852]"
46-
[29] "Cor (Intercept~zg: persons) | 1.0000 | [-1.00000, 1.00000]"
18+
[1] "# Fixed Effects (Count Model)"
19+
[2] ""
20+
[3] "Parameter | Log-Mean | SE | 95% CI | z | p"
21+
[4] "--------------------------------------------------------------------------"
22+
[5] "child | -1.0426 | 0.0899 | [-1.21878, -0.86645] | -11.5996 | < .001"
23+
[6] "camper [1] | 0.2684 | 0.0948 | [ 0.08262, 0.45421] | 2.8315 | 0.005 "
24+
[7] ""
25+
[8] "# Fixed Effects (Zero-Inflation Component)"
26+
[9] ""
27+
[10] "Parameter | Log-Odds | SE | 95% CI | z | p"
28+
[11] "-----------------------------------------------------------------------------"
29+
[12] "(Intercept) | -4.0420 | 0.2038 | [-4.44154, -3.64250] | -19.8292 | < .001"
30+
[13] "camper [1] | -4.8290 | 0.0002 | [-4.82949, -4.82861] | -21474.3719 | < .001"
31+
[14] ""
32+
[15] "# Random Effects Variances"
33+
[16] ""
34+
[17] "Parameter | Coefficient | 95% CI"
35+
[18] "--------------------------------------------------------------"
36+
[19] "SD (Intercept: persons) | 2.1800 | [1.15990, 4.09711]"
37+
[20] "SD (xb: persons) | 1.0592 | [0.59457, 1.88686]"
38+
[21] "Cor (Intercept~xb: persons) | -0.9883 | "
39+
[22] ""
40+
[23] "# Random Effects (Zero-Inflation Component)"
41+
[24] ""
42+
[25] "Parameter | Coefficient | 95% CI"
43+
[26] "----------------------------------------------------------------"
44+
[27] "SD (Intercept: persons) | 5.8635 | [ 3.83067, 8.97520]"
45+
[28] "SD (zg: persons) | 12.5180 | [10.94956, 14.31115]"
46+
[29] "Cor (Intercept~zg: persons) | 0.9790 | "
4747

4848
---
4949

5050
Code
5151
out[-c(5, 14)]
5252
Output
53-
[1] "# Fixed Effects (Count Model)"
54-
[2] ""
55-
[3] "Parameter | Log-Mean | SE | 95% CI | z | p"
56-
[4] "--------------------------------------------------------------------------"
57-
[5] "child | -1.0875 | 0.0981 | [-1.27967, -0.89528] | -11.0901 | < .001"
58-
[6] "camper [1] | 0.2723 | 0.1009 | [ 0.07461, 0.46999] | 2.6997 | 0.007 "
59-
[7] ""
60-
[8] "# Fixed Effects (Zero-Inflation Component)"
61-
[9] ""
62-
[10] "Parameter | Log-Odds | SE | 95% CI | z | p"
63-
[11] "-----------------------------------------------------------------------"
64-
[12] "(Intercept) | 1.8896 | 0.6642 | [ 0.58780, 3.19147] | 2.8449 | 0.004"
65-
[13] "camper [1] | -0.1701 | 0.3869 | [-0.92836, 0.58822] | -0.4396 | 0.660"
66-
[14] ""
67-
[15] "# Random Effects Variances"
68-
[16] ""
69-
[17] "Parameter | Coefficient | 95% CI"
70-
[18] "---------------------------------------------------------------"
71-
[19] "SD (Intercept: persons) | 3.4056 | [ 1.64567, 7.04777]"
72-
[20] "SD (xb: persons) | 1.2132 | [ 0.59190, 2.48650]"
73-
[21] "Cor (Intercept~xb: persons) | -1.0000 | [-1.00000, 1.00000]"
74-
[22] ""
75-
[23] "# Random Effects (Zero-Inflation Component)"
76-
[24] ""
77-
[25] "Parameter | Coefficient | 95% CI"
78-
[26] "---------------------------------------------------------------"
79-
[27] "SD (Intercept: persons) | 2.7358 | [ 1.16329, 6.43414]"
80-
[28] "SD (zg: persons) | 1.5683 | [ 0.64246, 3.82852]"
81-
[29] "Cor (Intercept~zg: persons) | 1.0000 | [-1.00000, 1.00000]"
53+
[1] "# Fixed Effects (Count Model)"
54+
[2] ""
55+
[3] "Parameter | Log-Mean | SE | 95% CI | z | p"
56+
[4] "--------------------------------------------------------------------------"
57+
[5] "child | -1.0426 | 0.0899 | [-1.21878, -0.86645] | -11.5996 | < .001"
58+
[6] "camper [1] | 0.2684 | 0.0948 | [ 0.08262, 0.45421] | 2.8315 | 0.005 "
59+
[7] ""
60+
[8] "# Fixed Effects (Zero-Inflation Component)"
61+
[9] ""
62+
[10] "Parameter | Log-Odds | SE | 95% CI | z | p"
63+
[11] "-----------------------------------------------------------------------------"
64+
[12] "(Intercept) | -4.0420 | 0.2038 | [-4.44154, -3.64250] | -19.8292 | < .001"
65+
[13] "camper [1] | -4.8290 | 0.0002 | [-4.82949, -4.82861] | -21474.3719 | < .001"
66+
[14] ""
67+
[15] "# Random Effects Variances"
68+
[16] ""
69+
[17] "Parameter | Coefficient | 95% CI"
70+
[18] "--------------------------------------------------------------"
71+
[19] "SD (Intercept: persons) | 2.1800 | [1.15990, 4.09711]"
72+
[20] "SD (xb: persons) | 1.0592 | [0.59457, 1.88686]"
73+
[21] "Cor (Intercept~xb: persons) | -0.9883 | "
74+
[22] ""
75+
[23] "# Random Effects (Zero-Inflation Component)"
76+
[24] ""
77+
[25] "Parameter | Coefficient | 95% CI"
78+
[26] "----------------------------------------------------------------"
79+
[27] "SD (Intercept: persons) | 5.8635 | [ 3.83067, 8.97520]"
80+
[28] "SD (zg: persons) | 12.5180 | [10.94956, 14.31115]"
81+
[29] "Cor (Intercept~zg: persons) | 0.9790 | "
8282

8383
# print-model_parameters glmmTMB CI alignment
8484

0 commit comments

Comments
 (0)