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
|
#' @param ss_function For Bayesian models, the function used to extract
#' sum-of-squares. Uses [`anova()`] by default, but can also be `car::Anova()`
#' for simple linear models.
#' @param draws For Bayesian models, an integer indicating the number of draws
#' from the posterior predictive distribution to return. Larger numbers take
#' longer to run, but provide estimates that are more stable.
#'
#' @export
#' @rdname eta_squared
eta_squared_posterior <- function(model,
partial = TRUE,
generalized = FALSE,
ss_function = stats::anova,
draws = 500,
verbose = TRUE,
...) {
UseMethod("eta_squared_posterior")
}
#' @export
#' @importFrom stats lm setNames
#' @importFrom insight find_formula get_predictors find_response check_if_installed
eta_squared_posterior.stanreg <- function(model,
partial = TRUE,
generalized = FALSE,
ss_function = stats::anova,
draws = 500,
verbose = TRUE,
...) {
insight::check_if_installed("rstantools")
mi <- .get_model_info(model, ...)
if (!mi$is_linear || mi$is_multivariate) {
insight::format_error("Computation of Eta Squared is only applicable to univariate linear models.")
}
if (mi$is_mixed) {
if (partial) {
if (verbose) {
insight::format_warning(
"Bayesian Partial Eta Squared not supported for mixed models.",
"Returning Eta Squared instead."
)
}
partial <- FALSE
# would need to account for random effects if present.
# Too hard right now.
}
if (isTRUE(generalized) || is.character(generalized)) {
if (verbose) {
insight::format_warning(
"Bayesian Generalized Eta Squared not supported for mixed models.",
"Returning Eta Squared instead."
)
}
generalized <- FALSE
}
}
## 1. get model data
f <- insight::find_formula(model)$conditional
X <- insight::get_predictors(model)
resp_name <- insight::find_response(model)
# test centered predictors
# if (verbose) .all_centered(X)
## 2. get ppd
ppd <- rstantools::posterior_predict(model,
draws = draws, # for rstanreg
nsamples = draws # for brms
)
## 3. Compute effect size...
if (verbose) {
insight::format_alert("Simulating effect size... This can take a while...")
}
res <- apply(ppd, 1, function(r) {
# sampled outcome + predictors
temp_dat <- X
temp_dat[[resp_name]] <- r
# fit a simple linear model
temp_fit <- stats::lm(f, temp_dat)
# compute effect size
if (isTRUE(verbose)) {
ANOVA <- ss_function(temp_fit, ...)
} else {
ANOVA <- suppressWarnings(ss_function(temp_fit, ...))
}
es <- eta_squared(ANOVA, ci = NULL, partial = partial, generalized = generalized, verbose = verbose)
es <- stats::setNames(
es[[if (partial) "Eta2_partial" else "Eta2"]],
es$Parameter
)
data.frame(t(es), check.names = FALSE)
})
res <- do.call("rbind", res)
attr(res, "generalized") <- generalized
return(res)
}
#' @export
eta_squared_posterior.brmsfit <- eta_squared_posterior.stanreg
#' #' @keywords internal
#' #' @importFrom stats contrasts
#' .all_centered <- function(X) {
#' numeric <- sapply(X, inherits, what = c("numeric", "integer"))
#' numerics <- colnames(X)[numeric]
#' factors <- colnames(X)[!numeric]
#'
#' numerics_centered <- factors_centered <- logical(0)
#'
#' if (length(numerics)) {
#' numerics_centered <- sapply(
#' X[, numerics, drop = FALSE],
#' function(xi) isTRUE(all.equal(mean(xi), 0))
#' )
#' }
#'
#'
#' of <- options()$contrasts
#' if (length(factors)) {
#' factors_centered <- sapply(X[, factors, drop = FALSE], function(xi) {
#' # if a contrast has negative and positive values, it is assumed to be one of:
#' # "contr.sum", "contr.helmert", "contr.poly", "contr.bayes"
#' (is.factor(xi) && (any(contrasts(xi) < 0) & any(contrasts(xi) > 0))) ||
#' # Or if it is not a factor, is the default method one of these?
#' (!is.factor(xi) && all(of %in% c("contr.sum", "contr.poly", "contr.bayes", "contr.helmert")))
#' })
#' }
#'
#'
#' if ((length(numerics_centered) && !all(numerics_centered)) ||
#' length(factors_centered) && !all(factors_centered)) {
#' non_centered <- !c(numerics_centered, factors_centered)
#' non_centered <- names(non_centered)[non_centered]
#' insight::format_warning(
#' "Not all variables are centered:\n ",
#' paste(non_centered, collapse = ", "),
#' "\n Results might be bogus if involved in interactions..."
#' )
#' }
#'
#' return(invisible(NULL))
#' }
|
