-
Notifications
You must be signed in to change notification settings - Fork 7
/
pvals.R
143 lines (120 loc) · 4.9 KB
/
pvals.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
#' @title Calculate Bootstrap p-values for fixed effects
#'
#' @description
#' Perform bootstrap tests based on the t-statistic for each
#' fixed effect in order to calculate approximate p-values.
#'
#' @details
#' The bootstrap test compares the fitted model specified by the user
#' to reduced models that eliminate a single fixed effect, the same
#' comparison summarized by the table of coefficients in the summary.
#' The bootstrap p-value is then calculated as
#' $(n_{extreme} + 1) / (B + 1)$.
#'
#' @return
#' A tibble giving the table of coefficients from the model summary with a column
#' appended containing bootstrap p-values.
#'
#' @references
#' Davison, A., & Hinkley, D. (1997). Tests. In Bootstrap Methods and their
#' Application (Cambridge Series in Statistical and Probabilistic Mathematics, pp. 136-190).
#' Cambridge: Cambridge University Press. doi:10.1017/CBO9780511802843.005
#'
#' @examples
#' \dontrun{
#' # This takes a while to run
#' bootstrap_pvals.merMod(jsp_mod, type = "wild", B = 1000, hccme = "hc2", aux.dist = "mammen")
#' }
#'
#' @inheritParams bootstrap
#' @export
bootstrap_pvals <- function(model, type, B, resample = NULL,
reb_type = NULL,
hccme = NULL,
aux.dist = NULL) {
if(!type %in% c("parametric", "residual", "case", "wild", "reb"))
stop("'type' must be one of 'parametric', 'residual', 'case', 'wild', or 'reb'")
if(!is.null(reb_type))
if(!reb_type %in% 0:2)
stop("'reb_type' must be either 0, 1, or 2")
UseMethod("bootstrap_pvals", model)
}
#' @rdname bootstrap_pvals
#' @export
#' @method bootstrap_pvals merMod
bootstrap_pvals.merMod <- function(model, type, B, resample = NULL,
reb_type = NULL,
hccme = NULL,
aux.dist = NULL) {
mod_summary <- summary(model)
tstats <- stats::coef(mod_summary)[, "t value"]
trms <- names(tstats)
# xmat <- lme4::getME(model, "X")
pvals <- rep(NA, length(tstats))
for(i in seq_along(pvals)) {
trm <- trms[i]
if(trm == "(Intercept)") null_mod <- update(model, . ~ . - 1)
else null_mod <- update(model, as.formula(paste(". ~ . -", trm)))
if(type == "case") {
gen_data <- bootstrap(null_mod, type = type, B = B, .refit = FALSE,
resample = resample, orig_data = model@frame)
} else {
gen_data <- bootstrap(null_mod, type = type, B = B, .refit = FALSE, resample = resample,
reb_type = reb_type, hccme = hccme, aux.dist = aux.dist)
}
extract_t <- function(x) stats::coef(summary(x))[i, "t value"]
if(type == "case") {
boot_ts <- purrr::map_dbl(gen_data, ~refit_to_newdf(model, newdata = .x, .f = extract_t))
} else{
boot_ts <- unlist(refit_merMod(gen_data, model, .f = extract_t)$tstar)
}
boot_ts <- boot_ts - mean(boot_ts, na.rm = TRUE)
pvals[i] <- (sum(abs(boot_ts) >= abs(tstats[i])) + 1) / (B + 1)
}
coef_tbl <- cbind(stats::coef(mod_summary), p.value = pvals) %>%
tibble::as_tibble(rownames = "term")
structure(
list(coefficients = coef_tbl, B = B, seed = .Random.seed, type = type),
class = "coef_tbl"
)
}
#' @rdname bootstrap_pvals
#' @export
#' @method bootstrap_pvals lme
bootstrap_pvals.lme <- function(model, type, B, resample = NULL,
reb_type = NULL,
hccme = NULL,
aux.dist = NULL) {
mod_summary <- summary(model)
tstats <- stats::coef(mod_summary)[, "t-value"]
trms <- names(tstats)
pvals <- rep(NA, length(tstats))
for(i in seq_along(pvals)) {
trm <- trms[i]
if(trm == "(Intercept)") null_mod <- update(model, . ~ . - 1)
else null_mod <- update(model, as.formula(paste(". ~ . -", trm)))
if(type == "case") {
gen_data <- bootstrap(null_mod, type = type, B = B, .refit = FALSE,
resample = resample)
} else {
gen_data <- bootstrap(null_mod, type = type, B = B, .refit = FALSE, resample = resample,
reb_type = reb_type, hccme = hccme, aux.dist = aux.dist)
}
extract_t <- function(x) stats::coef(summary(x))[i, "t-value"]
if(type == "case") {
boot_ts <- purrr::map(gen_data, ~updated.model(model, new.data = .x)) %>%
purrr::map_dbl(~extract_t(.x))
} else{
boot_ts <- purrr::map(gen_data, ~updated.model(model, new.y = .x)) %>%
purrr::map_dbl(~extract_t(.x))
}
boot_ts <- boot_ts - mean(boot_ts, na.rm = TRUE)
pvals[i] <- (sum(abs(boot_ts) >= abs(tstats[i])) + 1) / (B + 1)
}
coef_tbl <- cbind(stats::coef(mod_summary), p.value = pvals) %>%
tibble::as_tibble(rownames = "term")
structure(
list(coefficients = coef_tbl, B = B, seed = .Random.seed, type = type),
class = "coef_tbl"
)
}