power_marginaleffect {postcard} | R Documentation |
Power approximation for estimating marginal effects in GLMs
Description
The functions implements the algorithm for power estimation described in Powering RCTs for marginal effects with GLMs using prognostic score adjustment by HĂžjbjerre-Frandsen et. al (2025). See a bit more context in details and all details in reference.
Usage
power_marginaleffect(
response,
predictions,
target_effect,
exposure_prob,
var1 = NULL,
kappa1_squared = NULL,
estimand_fun = "ate",
estimand_fun_deriv0 = NULL,
estimand_fun_deriv1 = NULL,
inv_estimand_fun = NULL,
margin = estimand_fun(1, 1),
alpha = 0.05,
tolerance = 0.001,
verbose = options::opt("verbose"),
...
)
Arguments
response |
a vector of mode |
predictions |
a vector of mode |
target_effect |
a |
exposure_prob |
a |
var1 |
a |
kappa1_squared |
a |
estimand_fun |
a |
estimand_fun_deriv0 |
a |
estimand_fun_deriv1 |
a |
inv_estimand_fun |
(optional) a |
margin |
a |
alpha |
a |
tolerance |
passed to all.equal when comparing calculated |
verbose |
|
... |
arguments passed to stats::uniroot, which is used to find the
inverse of |
Details
The reference in the description shows in its "Prospective power" section a derivation of a variance bound
v_\infty^{\uparrow 2} = r_0'^{\, 2}\sigma_0^2+
r_1'^{\, 2}\sigma_1^2+
\pi_0\pi_1\left(\frac{|r_0'|\kappa_0}{\pi_0} +
\frac{|r_1'|\kappa_1}{\pi_1} \right)^2
where r_a'
is the derivative of the estimand_fun
with respect to
\Psi_a
, \sigma_a^2
is the variance of the potential outcome corresponding to
group a
, \pi_a
is the probability of being assigned to group a
,
and \kappa_a
is the expected mean-squared error when predicting the
potential outcome corresponding to group a
.
The variance bound is then used for calculating a lower bound of the power using
the distributions corresponding to the null and alternative hypotheses
\mathcal{H}_0: \hat{\Psi} \sim F_0 = \mathcal{N}(\Delta ,v_\infty^{\uparrow 2} / n)
and
\mathcal{H}_1: \hat{\Psi} \sim F_1 = \mathcal{N}(\Psi,v_\infty^{\uparrow 2} / n)
.
See more details in the reference.
Relating arguments to formulas
-
response
: Used to obtain both\sigma_0^2
(by taking the sample variance of the response) and\kappa_0
. -
predictions
: Used when calculating the MSE\kappa_0
. -
var1
:\sigma_1^2
. As a default, chosen to be the same as\sigma_0^2
. Can specify differently through this argument fx. byInflating or deflating the value of
\sigma_0^2
by a scalar according to prior beliefs. Fx. specifyvar1 = function(x) 1.2 * x
to inflate\sigma_0^2
by 1.2.If historical data is available for group 1, an estimate of the variance from that data can be provided here.
-
kappa1_squared
:\kappa_1
. Same as forvar1
, default is to assume the same value askappa0_squared
, which is found by using theresponse
andpredictions
vectors. Adjust the value according to prior beliefs if relevant. -
target_effect
:\Psi
. -
exposure_prob
:\pi_1
Value
a numeric
with the estimated power.
See Also
See power_linear for power approximation functionalities for linear models.
Examples
# Generate a data set to use as an example
n <- 100
exposure_prob <- .5
dat_gaus <- glm_data(Y ~ 1+2*X1-X2+3*A+1.6*A*X2,
X1 = rnorm(n),
X2 = rgamma(n, shape = 2),
A = rbinom(n, size = 1, prob = exposure_prob),
family = gaussian())
# ---------------------------------------------------------------------------
# Obtain out-of-sample (OOS) prediction using glm model
# ---------------------------------------------------------------------------
gaus1 <- dat_gaus[1:(n/2), ]
gaus2 <- dat_gaus[(n/2+1):n, ]
glm1 <- glm(Y ~ X1 + X2 + A, data = gaus1)
glm2 <- glm(Y ~ X1 + X2 + A, data = gaus2)
preds_glm1 <- predict(glm2, newdata = gaus1, type = "response")
preds_glm2 <- predict(glm1, newdata = gaus2, type = "response")
preds_glm <- c(preds_glm1, preds_glm2)
# Obtain power
power_marginaleffect(
response = dat_gaus$Y,
predictions = preds_glm,
target_effect = 2,
exposure_prob = exposure_prob
)
# ---------------------------------------------------------------------------
# Get OOS predictions using discrete super learner and adjust variance
# ---------------------------------------------------------------------------
learners <- list(
mars = list(
model = parsnip::set_engine(
parsnip::mars(
mode = "regression", prod_degree = 3
),
"earth"
)
),
lm = list(
model = parsnip::set_engine(
parsnip::linear_reg(),
"lm"
)
)
)
lrnr1 <- fit_best_learner(preproc = list(mod = Y ~ X1 + X2 + A),
data = gaus1,
learners = learners)
lrnr2 <- fit_best_learner(preproc = list(mod = Y ~ X1 + X2 + A),
data = gaus2,
learners = learners)
preds_lrnr1 <- dplyr::pull(predict(lrnr2, new_data = gaus1))
preds_lrnr2 <- dplyr::pull(predict(lrnr1, new_data = gaus2))
preds_lrnr <- c(preds_lrnr1, preds_lrnr2)
# Estimate the power AND adjust the assumed variance in the "unknown"
# group 1 to be 20% larger than in group 0
power_marginaleffect(
response = dat_gaus$Y,
predictions = preds_lrnr,
target_effect = 2,
exposure_prob = exposure_prob,
var1 = function(var0) 1.2 * var0
)