For fitted models, this command calculates (1) how much bias there must be in an estimate to invalidate/sustain an inference; (2) the impact of an omitted variable necessary to invalidate/sustain an inference for a regression coefficient. Currently works for: models created with lm() (linear models).
konfound(
model_object,
tested_variable,
alpha = 0.05,
tails = 2,
index = "RIR",
to_return = "print",
test_all = FALSE,
two_by_two = FALSE,
n_treat = NULL,
switch_trm = TRUE,
replace = "control"
)output from a model (currently works for: lm)
Variable associated with the unstandardized beta coefficient to be tested
probability of rejecting the null hypothesis (defaults to 0.05)
integer whether hypothesis testing is one-tailed (1) or two-tailed (2; defaults to 2)
whether output is RIR or IT (impact threshold); defaults to "RIR"
whether to return a data.frame (by specifying this argument to equal "raw_output" for use in other analyses) or a plot ("plot"); default is to print ("print") the output to the console; can specify a vector of output to return
whether to carry out the sensitivity test for all of the coefficients (defaults to FALSE)
whether or not the tested variable is a dichotomous variable in a GLM; if so, the 2X2 table approach is used; only works for single variables at present (so test_all = TRUE will return an error)
the number of cases associated with the treatment condition; applicable only when model_type = "logistic"
whether to switch the treatment and control cases; defaults to FALSE; applicable only when model_type = "logistic"
whether using entire sample or the control group to calculate the base rate; default is the control group
prints the bias and the number of cases that would have to be replaced with cases for which there is no effect to invalidate the inference
# using lm() for linear models
m1 <- lm(mpg ~ wt + hp, data = mtcars)
konfound(m1, wt)
#> Robustness of Inference to Replacement (RIR):
#> To invalidate an inference, 66.521 % of the estimate would have to be due to bias.
#> This is based on a threshold of -1.298 for statistical significance (alpha = 0.05).
#>
#> To invalidate an inference, 21 observations would have to be replaced with cases
#> for which the effect is 0 (RIR = 21).
#>
#> See Frank et al. (2013) for a description of the method.
#>
#> Citation: Frank, K.A., Maroulis, S., Duong, M., and Kelcey, B. (2013).
#> What would it take to change an inference?
#> Using Rubin's causal model to interpret the robustness of causal inferences.
#> Education, Evaluation and Policy Analysis, 35 437-460.
#> For more detailed output, consider setting `to_return` to table
#> To consider other predictors of interest, consider setting `test_all` to TRUE.
konfound(m1, wt, test_all = TRUE)
#> Note that this output is calculated based on the correlation-based approach used in mkonfound()
#> Warning: row names were found from a short variable and have been discarded
#> # A tibble: 2 × 8
#> var_name t df action inference pct_bias_to_change_i…¹ itcv r_con
#> <chr> <dbl> <dbl> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 wt -6.13 29 to_invalida… reject_n… 52.7 0.614 0.784
#> 2 hp -3.52 29 to_invalida… reject_n… 35.1 0.298 0.546
#> # ℹ abbreviated name: ¹pct_bias_to_change_inference
konfound(m1, wt, to_return = "table")
#> Dependent variable is mpg
#> Warning: Unknown or uninitialised column: `itcv`.
#> Warning: Unknown or uninitialised column: `impact`.
#> # A tibble: 3 × 7
#> term estimate std.error statistic p.value itcv impact
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 (Intercept) 37.2 1.60 23.3 0 NA NA
#> 2 wt -3.88 0.633 -6.13 0 0.625 NA
#> 3 hp -0.032 0.009 -3.52 0.001 NA 0.511
# using glm() for non-linear models
if (requireNamespace("forcats")) {
d <- forcats::gss_cat
d$married <- ifelse(d$marital == "Married", 1, 0)
m2 <- glm(married ~ age, data = d, family = binomial(link = "logit"))
konfound(m2, age)
}
#> Loading required namespace: forcats
#> Note that for a non-linear model, impact threshold should not be interpreted.
#> Note that this is only an approximation. For exact results in terms of the number of cases that must be switched from treatment success to treatment failure to invalidate the inference see: https://msu.edu/~kenfrank/non%20linear%20replacement%20treatment.xlsm
#> If a dichotomous independent variable is used, consider using the 2X2 table approach enabled with the argument two_by_two = TRUE
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'd' not found
# using lme4 for mixed effects (or multi-level) models
if (requireNamespace("lme4")) {
library(lme4)
m3 <- fm1 <- lme4::lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
konfound(m3, Days)
}
#> Loading required namespace: lme4
#> Loading required package: Matrix
#> Robustness of Inference to Replacement (RIR):
#> To invalidate an inference, 84.825 % of the estimate would have to be due to bias.
#> This is based on a threshold of 1.588 for statistical significance (alpha = 0.05).
#>
#> To invalidate an inference, 137 observations would have to be replaced with cases
#> for which the effect is 0 (RIR = 137).
#>
#> See Frank et al. (2013) for a description of the method.
#>
#> Citation: Frank, K.A., Maroulis, S., Duong, M., and Kelcey, B. (2013).
#> What would it take to change an inference?
#> Using Rubin's causal model to interpret the robustness of causal inferences.
#> Education, Evaluation and Policy Analysis, 35 437-460.
#> Note that the Kenward-Roger approximation is used to estimate degrees of freedom for the predictor(s) of interest. We are presently working to add other methods for calculating the degrees of freedom for the predictor(s) of interest. If you wish to use other methods now, consider those detailed here: https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#why-doesnt-lme4-display-denominator-degrees-of-freedomp-values-what-other-options-do-i-have. You can then enter degrees of freedom obtained from another method along with the coefficient, number of observations, and number of covariates to the pkonfound() function to quantify the robustness of the inference.
#> NULL
m4 <- glm(outcome ~ condition, data = binary_dummy_data, family = binomial(link = "logit"))
konfound(m4, condition, two_by_two = TRUE, n_treat = 55)
#> Conclusion:
#> User-entered Table:
#> Fail Success
#> Control 36 16
#> Treatment 18 37
#>
#> Note: Values have been rounded to the nearest integer.This may cause a little change to the estimated effect for the Implied Table.
#>
#> To invalidate the inference, you would need to replace 15 treatment success cases
#> with cases for which the probability of failure in the control group applies (RIR = 15).
#> This is equivalent to transferring 10 cases from treatment success to treatment failure,
#> as shown, from the Implied Table to the Transfer Table.
#>
#> Transfer Table:
#> Fail Success
#> Control 36 16
#> Treatment 28 27
#>
#> For the Implied Table, we have an estimate of 1.531, with a SE of 0.416 and a t-ratio of 3.684.
#> For the Transfer Table, we have an estimate of 0.775, with a SE of 0.404 and a t-ratio of 1.918.
#>
#> RIR:
#> RIR = 15
#> NULL