Skip to content

Commit 4e957f3

Browse files
strengejackestrengejacke
committed
refactor, tests
1 parent 4eb65a0 commit 4e957f3

3 files changed

Lines changed: 137 additions & 78 deletions

File tree

‎R/format_p_adjust.R‎

Lines changed: 109 additions & 75 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,8 @@ format_p_adjust <- function(method) {
3131
}
3232

3333

34+
# p-value adjustment -----
35+
3436
.p_adjust <- function(params, p_adjust, model = NULL, verbose = TRUE) {
3537
# check if we have any adjustment at all, and a p-column
3638
if (!is.null(p_adjust) && "p" %in% colnames(params) && p_adjust != "none") {
@@ -42,20 +44,7 @@ format_p_adjust <- function(method) {
4244
# for interaction terms, e.g. for "by" argument in emmeans
4345
# pairwise comparison, we have to adjust the rank resp. the
4446
# number of estimates in a comparison family
45-
rank_adjust <- tryCatch(
46-
{
47-
correction <- 1
48-
by_vars <- model@misc$by.vars
49-
if (!is.null(by_vars) && by_vars %in% colnames(params)) {
50-
correction <- insight::n_unique(params[[by_vars]])
51-
}
52-
correction
53-
},
54-
error = function(e) {
55-
1
56-
}
57-
)
58-
47+
rank_adjust <- .p_adjust_rank(model, params)
5948

6049
# only proceed if valid argument-value
6150
if (tolower(p_adjust) %in% tolower(all_methods)) {
@@ -70,74 +59,18 @@ format_p_adjust <- function(method) {
7059
params$p <- stats::p.adjust(params$p, method = p_adjust)
7160
} else if (tolower(p_adjust) == "tukey") {
7261
# tukey adjustment
73-
if ("df" %in% colnames(params) && length(stat_column) > 0) {
74-
params$p <- suppressWarnings(stats::ptukey(
75-
sqrt(2) * abs(params[[stat_column]]),
76-
nrow(params) / rank_adjust,
77-
params$df,
78-
lower.tail = FALSE
79-
))
80-
# for specific contrasts, ptukey might fail, and the tukey-adjustement
81-
# could just be simple p-value calculation
82-
if (all(is.na(params$p))) {
83-
params$p <- 2 * stats::pt(abs(params[[stat_column]]), df = params$df, lower.tail = FALSE)
84-
verbose <- FALSE
85-
}
86-
}
62+
result <- .p_adjust_tukey(params, stat_column, rank_adjust, verbose)
63+
params <- result$params
64+
verbose <- result$verbose
8765
} else if (tolower(p_adjust) == "scheffe" && !is.null(model)) {
8866
# scheffe adjustment
89-
if ("df" %in% colnames(params) && length(stat_column) > 0) {
90-
# 1st try
91-
scheffe_ranks <- try(qr(model@linfct)$rank, silent = TRUE)
92-
93-
# 2nd try
94-
if (inherits(scheffe_ranks, "try-error") || is.null(scheffe_ranks)) {
95-
scheffe_ranks <- try(model$qr$rank, silent = TRUE)
96-
}
97-
98-
if (inherits(scheffe_ranks, "try-error") || is.null(scheffe_ranks)) {
99-
scheffe_ranks <- nrow(params)
100-
}
101-
scheffe_ranks <- scheffe_ranks / rank_adjust
102-
params$p <- stats::pf(params[[stat_column]]^2 / scheffe_ranks,
103-
df1 = scheffe_ranks,
104-
df2 = params$df,
105-
lower.tail = FALSE
106-
)
107-
}
67+
params <- .p_adjust_scheffe(model, params, stat_column, rank_adjust)
10868
} else if (tolower(p_adjust) == "sidak") {
10969
# sidak adjustment
11070
params$p <- 1 - (1 - params$p)^(nrow(params) / rank_adjust)
11171
} else if (tolower(p_adjust) == "sup-t") {
11272
# sup-t adjustment
113-
insight::check_if_installed("mvtnorm")
114-
# get correlation matrix, based on the covariance matrix
115-
vc <- .safe(stats::cov2cor(insight::get_varcov(model)))
116-
if (is.null(vc)) {
117-
insight::format_warning("Could not calculate covariance matrix for `sup-t` adjustment.")
118-
return(params)
119-
}
120-
# get confidence interval level, or set default
121-
ci_level <- .safe(params$CI[1])
122-
if (is.null(ci_level)) {
123-
ci_level <- 0.95
124-
}
125-
# calculate updated confidence interval level, based on simultaenous
126-
# confidence intervals (https://onlinelibrary.wiley.com/doi/10.1002/jae.2656)
127-
crit <- mvtnorm::qmvt(ci_level, df = params$df[1], tail = "both.tails", corr = vc)$quantile
128-
ci_level <- 1 - 2 * stats::pt(-abs(crit), df = params$df[1])
129-
# update confidence intervals
130-
params$CI_low <- params$Coefficient - crit * params$SE
131-
params$CI_high <- params$Coefficient + crit * params$SE
132-
# udpate p-values
133-
for (i in 1:nrow(params)) {
134-
params$p[i] <- 1 - mvtnorm::pmvt(
135-
lower = rep(-abs(stats::qt(params$p[i] / 2, df = params$df[i])), nrow(vc)),
136-
upper = rep(abs(stats::qt(params$p[i] / 2, df = params$df[i])), nrow(vc)),
137-
corr = vc,
138-
df = params$df[i]
139-
)
140-
}
73+
params <- .p_adjust_supt(model, params)
14174
}
14275

14376
if (isTRUE(all(old_p_vals == params$p)) && !identical(p_adjust, "none") && verbose) {
@@ -149,3 +82,104 @@ format_p_adjust <- function(method) {
14982
}
15083
params
15184
}
85+
86+
87+
# calculate rank adjustment -----
88+
89+
.p_adjust_rank <- function(model, params) {
90+
tryCatch(
91+
{
92+
correction <- 1
93+
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]])
96+
}
97+
correction
98+
},
99+
error = function(e) {
100+
1
101+
}
102+
)
103+
}
104+
105+
106+
# tukey adjustment -----
107+
108+
.p_adjust_tukey <- function(params, stat_column, rank_adjust = 1, verbose = TRUE) {
109+
if ("df" %in% colnames(params) && length(stat_column) > 0) {
110+
params$p <- suppressWarnings(stats::ptukey(
111+
sqrt(2) * abs(params[[stat_column]]),
112+
nrow(params) / rank_adjust,
113+
params$df,
114+
lower.tail = FALSE
115+
))
116+
# for specific contrasts, ptukey might fail, and the tukey-adjustement
117+
# could just be simple p-value calculation
118+
if (all(is.na(params$p))) {
119+
params$p <- 2 * stats::pt(abs(params[[stat_column]]), df = params$df, lower.tail = FALSE)
120+
verbose <- FALSE
121+
}
122+
}
123+
list(params = params, verbose = verbose)
124+
}
125+
126+
127+
# scheffe adjustment -----
128+
129+
.p_adjust_scheffe <- function(model, params, stat_column, rank_adjust = 1) {
130+
if ("df" %in% colnames(params) && length(stat_column) > 0) {
131+
# 1st try
132+
scheffe_ranks <- try(qr(model@linfct)$rank, silent = TRUE)
133+
134+
# 2nd try
135+
if (inherits(scheffe_ranks, "try-error") || is.null(scheffe_ranks)) {
136+
scheffe_ranks <- try(model$qr$rank, silent = TRUE)
137+
}
138+
139+
if (inherits(scheffe_ranks, "try-error") || is.null(scheffe_ranks)) {
140+
scheffe_ranks <- nrow(params)
141+
}
142+
scheffe_ranks <- scheffe_ranks / rank_adjust
143+
params$p <- stats::pf(params[[stat_column]]^2 / scheffe_ranks,
144+
df1 = scheffe_ranks,
145+
df2 = params$df,
146+
lower.tail = FALSE
147+
)
148+
}
149+
params
150+
}
151+
152+
153+
# sup-t adjustment -----
154+
155+
.p_adjust_supt <- function(model, params) {
156+
insight::check_if_installed("mvtnorm")
157+
# get correlation matrix, based on the covariance matrix
158+
vc <- .safe(stats::cov2cor(insight::get_varcov(model)))
159+
if (is.null(vc)) {
160+
insight::format_warning("Could not calculate covariance matrix for `sup-t` adjustment.")
161+
return(params)
162+
}
163+
# get confidence interval level, or set default
164+
ci_level <- .safe(params$CI[1])
165+
if (is.null(ci_level)) {
166+
ci_level <- 0.95
167+
}
168+
# calculate updated confidence interval level, based on simultaenous
169+
# confidence intervals (https://onlinelibrary.wiley.com/doi/10.1002/jae.2656)
170+
crit <- mvtnorm::qmvt(ci_level, df = params$df[1], tail = "both.tails", corr = vc)$quantile
171+
ci_level <- 1 - 2 * stats::pt(-abs(crit), df = params$df[1])
172+
# update confidence intervals
173+
params$CI_low <- params$Coefficient - crit * params$SE
174+
params$CI_high <- params$Coefficient + crit * params$SE
175+
# udpate p-values
176+
for (i in 1:nrow(params)) {
177+
params$p[i] <- 1 - mvtnorm::pmvt(
178+
lower = rep(-abs(stats::qt(params$p[i] / 2, df = params$df[i])), nrow(vc)),
179+
upper = rep(abs(stats::qt(params$p[i] / 2, df = params$df[i])), nrow(vc)),
180+
corr = vc,
181+
df = params$df[i]
182+
)
183+
}
184+
params
185+
}

‎tests/testthat/_snaps/p_adjust.md‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
# model_parameters, simultaenous confidence intervals
2+
3+
Code
4+
print(out, zap_small = TRUE)
5+
Output
6+
Parameter | Coefficient | SE | 95% CI | t(29) | p
7+
------------------------------------------------------------------
8+
(Intercept) | 37.23 | 1.60 | [33.32, 41.14] | 23.28 | < .001
9+
wt | -3.88 | 0.63 | [-5.42, -2.33] | -6.13 | < .001
10+
hp | -0.03 | 0.01 | [-0.05, -0.01] | -3.52 | 0.003
11+
12+
p-value adjustment method: Simultaneous confidence bands
13+
Message
14+
15+
Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
16+
using a Wald t-distribution approximation.
17+

‎tests/testthat/test-p_adjust.R‎

Lines changed: 11 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -45,15 +45,23 @@ test_that("model_parameters, p-adjust after keep/drop", {
4545
)
4646
})
4747

48+
4849
test_that("model_parameters, emmeans, p-adjust", {
4950
skip_if_not_installed("emmeans")
5051
m <- pairs(emmeans::emmeans(aov(Sepal.Width ~ Species, data = iris), ~Species))
5152
mp <- model_parameters(m)
5253
expect_equal(mp$p, as.data.frame(m)$p.value, tolerance = 1e-4)
53-
})
54-
test_that("model_parameters, emmeans, p-adjust", {
55-
skip_if_not_installed("emmeans")
54+
5655
m <- pairs(emmeans::emmeans(aov(Sepal.Width ~ Species, data = iris), ~Species), adjust = "scheffe")
5756
mp <- model_parameters(m, p_adjust = "scheffe")
5857
expect_equal(mp$p, as.data.frame(m)$p.value, tolerance = 1e-4)
5958
})
59+
60+
61+
test_that("model_parameters, simultaenous confidence intervals", {
62+
skip_if_not_installed("mvtnorm")
63+
m <- lm(mpg ~ wt + hp, data = mtcars)
64+
set.seed(123)
65+
out <- model_parameters(m, p_adjust = "sup-t")
66+
expect_snapshot(print(out, zap_small = TRUE))
67+
})

0 commit comments

Comments
 (0)