/
svalue.R
311 lines (265 loc) · 11.6 KB
/
svalue.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
#' Severity of publication bias needed to "explain away" results
#'
#' Estimates the S-value, defined as the severity of publication bias (i.e., the
#' ratio by which affirmative studies are more likely to be published than
#' nonaffirmative studies) that would be required to shift the pooled point
#' estimate or its confidence interval limit to the value `q`.
#'
#' @inheritParams metabias::params
#' @inheritParams pubbias_meta
#' @param selection_ratio_max The largest value of `selection_ratio` that should
#' be included in the grid search. This argument is only needed when
#' `model_type = "robust"`.
#'
#' @details To illustrate interpretation of the S-value, if the S-value for the
#' point estimate is 30 with `q=0`, this indicates that affirmative studies
#' (i.e., those with a "statistically significant" and positive estimate)
#' would need to be 30-fold more likely to be published than nonaffirmative
#' studies (i.e., those with a "nonsignificant" or negative estimate) to
#' attenuate the pooled point estimate to `q`.
#'
#' If `favor_positive = FALSE`, such that publication bias is assumed to favor
#' negative rather than positive estimates, the signs of `yi` will be reversed
#' prior to performing analyses. The returned number of affirmative and
#' nonaffirmative studies will reflect the recoded signs.
#'
#' If `return_worst_meta = TRUE`, also returns the worst-case meta-analysis of
#' only the nonaffirmative studies. If `model_type = "fixed"`, the worst-case
#' meta-analysis is fit by `metafor::rma.uni()`. If `model_type = "robust"`,
#' it is fit by `robumeta::robu()`. Note that in the latter case, custom
#' inverse-variance weights are used, which are the inverse of the sum of the
#' study's variance and a heterogeneity estimate from a naive random-effects
#' meta-analysis (Mathur & VanderWeele, 2020). This is done for consistency
#' with the results of `pubbias_meta()`, which is used to determine `sval_est`
#' and `sval_ci`. Therefore, the worst-case meta-analysis results may differ
#' slightly from what you would obtain if you simply fit `robumeta::robu()` on
#' the nonaffirmative studies with the default weights.
#'
#' @return An object of class [metabias::metabias()], a list containing:
#' \describe{
#' \item{data}{A tibble with one row per study and the columns
#' `r meta_names_str("data")`.}
#' \item{values}{A list with the elements `r meta_names_str("values")`.}
#' \item{stats}{A tibble with the columns `r meta_names_str("stats")`.
#' `sval_est` represents the amount of publication bias required
#' to attenuate the pooled point estimate to `q`; `sval_ci`
#' represents the amount of publication bias required to
#' attenuate the confidence interval limit of the pooled point
#' estimate to `q`.}
#' \item{fit}{A list of fitted models, if any.}
#' }
#'
#' @references
#' \insertRef{mathur2020}{metabias}
#'
#' @export
#' @example inst/examples/pubbias_svalue.R
pubbias_svalue <- function(yi, # data
vi,
sei,
cluster = 1:length(yi),
q = 0, # params
model_type = "robust", # opts
favor_positive = TRUE,
alpha_select = 0.05,
ci_level = 0.95,
small = TRUE,
selection_ratio_max = 200,
return_worst_meta = FALSE) {
# error if selection_ratio_max doesn't make sense
if (selection_ratio_max < 1) stop("selection_ratio_max must be at least 1.")
# resolve vi and sei
if (missing(vi)) {
if (missing(sei)) {
stop("Must specify 'vi' or 'sei' argument.")
}
vi <- sei ^ 2
}
# number of point estimates
k <- length(yi)
# calculate alpha for inference on point estimate
alpha <- 1 - ci_level
# warn if clusters are present but model_type == "fixed"
nclusters <- length(unique(cluster))
if (nclusters < k && model_type == "fixed") {
warning("You indicated there are clusters, but these will be ignored due to fixed-effects specification. To accommodate clusters, instead choose model_type = robust.")
}
# fit uncorrected model
m0 <- pubbias_meta(yi = yi,
vi = vi,
sei = sei,
cluster = cluster,
selection_ratio = 1,
selection_tails = 1,
model_type = model_type,
favor_positive = favor_positive,
ci_level = ci_level,
small = small)
# error if q is on wrong side of null
est0 <- m0$stats$estimate
q_error <- function(dir) {
glue("The uncorrected pooled point estimate is {est0}. q must be {dir} than this value (i.e., closer to zero).")
}
if (est0 > 0 && q > est0) stop(q_error("less"))
if (est0 < 0 && q < est0) stop(q_error("greater"))
##### Flip Estimate Signs If Needed #####
if (!favor_positive) {
yi <- -yi
q <- -q
}
# 2-sided p-values for each study even if 1-tailed selection
pvals <- 2 * (1 - pnorm(abs(yi) / sqrt(vi)))
# affirmative indicator under 1-tailed selection
affirm <- (pvals < alpha_select) & (yi > 0)
k_affirmative <- sum(affirm)
k_nonaffirmative <- k - k_affirmative
# error if there are zero affirmative or zero nonaffirmative studies
k_zero_msg <- \(dir) glue("There are zero {dir} studies. Model estimation cannot proceed.")
if (k_affirmative == 0) stop(k_zero_msg("affirmative"))
if (k_nonaffirmative == 0) stop(k_zero_msg("nonaffirmative"))
dat <- data.frame(yi, vi, affirm, cluster)
fits <- list()
# fit worst-case meta-analysis of only nonaffirmative studies
# always need model for svalue search initilization
meta_worst <- fit_meta_worst(dat, model_type = model_type,
ci_level = ci_level, small = small)
if (return_worst_meta) fits$meta_worst <- meta_worst$meta
##### Fixed-Effects Model #####
if (model_type == "fixed") {
strat <- dat |>
group_by(.data$affirm) |>
summarise(nu = sum(1 / .data$vi), ybar = sum(.data$yi / .data$vi))
# components of bias-corrected estimate by affirmative status
ybar_n <- strat$ybar[!strat$affirm]
ybar_a <- strat$ybar[strat$affirm]
nu_n <- strat$nu[!strat$affirm]
nu_a <- strat$nu[strat$affirm]
# S-value for point estimate
sval_est <- (nu_a * q - ybar_a) / (ybar_n - nu_n * q)
# S-value for CI (to shift it to q)
# match term names used in Wolfram Alpha
a <- ybar_n
b <- ybar_a
c <- nu_n
d <- nu_a
k <- if (!small) qnorm(1 - (alpha / 2)) else qt(1 - (alpha / 2),
df = k - 1)
term_a <- k ^ 2 * (a ^ 2 * d -
(2 * c * d * q) * (a + b) +
b ^ 2 * c +
q ^ 2 * (c ^ 2 * d + d ^ 2 * c) -
c * d * k ^ 2)
term_b <- -a * b + a * d * q + b * c * q - c * d * q ^ 2
term_c <- a ^ 2 - 2 * a * c * q + c ^ 2 * q ^ 2 - c * k ^ 2
sval_ci <- (-sqrt(term_a) + term_b) / term_c
if (sval_ci < 0) sval_ci <- (sqrt(term_a) + term_b) / term_c
} # end fixed = TRUE
##### Robust Independent and Robust Clustered #####
if (model_type == "robust") {
.meta_fun <- function(.selection_ratio) {
pubbias_meta(yi = yi,
vi = vi,
sei = sei,
cluster = cluster,
selection_ratio = .selection_ratio,
selection_tails = 1,
model_type = model_type,
favor_positive = TRUE, # always TRUE because we've already flipped signs if needed
ci_level = ci_level,
small = small)
}
.svalue_fun <- function(.param) {
find_svalue(param = .param, meta_fun = .meta_fun,
meta_worst = meta_worst, q = q,
selection_ratio_max = selection_ratio_max)
}
sval_est <- .svalue_fun("estimate")
sval_ci <- .svalue_fun("ci_lower")
}
# s-values less than 1 indicate complete robustness
# is.numeric is in case we have a "< XXX" string instead of a number
.pos_sval <- \(sval) if (is.numeric(sval) & sval <= 1) "Not possible" else sval
sval_est <- .pos_sval(sval_est)
sval_ci <- .pos_sval(sval_ci)
# check if naive confidence interval contains q
# m0 was fit BEFORE flipping signs
# but q has now been flipped in the latter case in "or" statement below
if ((est0 > 0 && m0$stats$ci_lower < q) ||
(est0 < 0 && m0$stats$ci_upper > -q)) {
# important: Shiny website assumes that this exact string ("--") for CI can
# be interpreted as the naive CI's already containing q
sval_ci <- "--"
message("sval_ci is not applicable because the naive confidence interval already contains q")
}
values <- list(
q = q,
model_type = model_type,
favor_positive = favor_positive,
alpha_select = alpha_select,
ci_level = ci_level,
small = small,
selection_ratio_max = selection_ratio_max,
k = k,
k_affirmative = k_affirmative,
k_nonaffirmative = k_nonaffirmative
)
stats <- tibble(sval_est = sval_est, sval_ci = sval_ci)
metabias::metabias(data = dat, values = values, stats = stats, fits = fits)
}
#' @keywords internal
find_svalue <- function(param, meta_fun, meta_worst, q, selection_ratio_max,
tol = 0.0001) { # param is estimate or ci_lower
if (meta_worst$stats[[param]] > q) return("Not possible")
# define the function we need to minimize
# i.e., distance between corrected estimate and the target value of q
func <- function(.selection_ratio) {
corrected <- suppressWarnings(meta_fun(.selection_ratio))
return(abs(corrected$stats[[param]] - q))
}
opt <- optimize(f = func,
interval = c(1, selection_ratio_max),
maximum = FALSE)
sval <- opt$minimum
# discrepancy between the corrected estimate and the s-value
diff <- opt$objective
# if the optimal value is very close to the upper range of grid search
# AND we're still not very close to the target q,
# that means the optimal value was above selection_ratio_max
if (abs(sval - selection_ratio_max) < tol && diff > tol)
sval <- paste(">", selection_ratio_max)
return(sval)
}
#' @rdname pubbias_svalue
#' @param clustervar (deprecated) see cluster
#' @param model (deprecated) see model_type
#' @param alpha.select (deprecated) see alpha_select
#' @param eta.grid.hi (deprecated) see selection_ratio_max
#' @param favor.positive (deprecated) see favor_positive
#' @param CI.level (deprecated) see ci_level
#' @param return.worst.meta (deprecated) see return_worst_meta
#' @keywords internal
#' @export
svalue <- function(yi,
vi,
q,
clustervar = 1:length(yi),
model,
alpha.select = 0.05,
eta.grid.hi = 200,
favor.positive,
CI.level = 0.95,
small = TRUE,
return.worst.meta = FALSE) {
lifecycle::deprecate_warn("2.3.0", "svalue()", "pubbias_svalue()")
pubbias_svalue(yi = yi,
vi = vi,
cluster = clustervar,
q = q,
model_type = model,
favor_positive = favor.positive,
alpha_select = alpha.select,
ci_level = CI.level,
small = small,
selection_ratio_max = eta.grid.hi,
return_worst_meta = return.worst.meta)
}