| Title: | Robust Effect Size Index (RESI) Estimation |
|---|---|
| Description: | Summarize model output using a robust effect size index. The index is introduced in Vandekar, Tao, & Blume (2020) <doi:10.1007/s11336-020-09698-2>. |
| Authors: | Megan Jones [aut, cre], Kaidi Kang [aut], Simon Vandekar [aut] |
| Maintainer: | Megan Jones <[email protected]> |
| License: | GPL-3 |
| Version: | 1.2.4 |
| Built: | 2024-10-31 22:08:48 UTC |
| Source: | https://github.com/statimagcoll/resi |
After running the resi function on a fitted model, this function can be used to print the Anova-style table component. If the resi function was run with the ‘store.boot = TRUE' option to store the full matrix of bootstrapped estimates, the user can specify a different alpha level for this function’s confidence intervals.
## S3 method for class 'resi' anova(object, alpha = NULL, ...)## S3 method for class 'resi' anova(object, alpha = NULL, ...)
object |
an object resulting from resi function |
alpha |
an optional new specification for the confidence level. Can be vector-valued |
... |
ignored |
The resi function uses the car::Anova function to compute the Anova table.
Returns an 'anova' object containing the computed Anova-style table
# fit a model mod = lm(charges ~ bmi + sex, data = RESI::insurance) # run resi with the store.boot = TRUE option resi.obj = resi(mod, nboot = 100, store.boot = TRUE, alpha = 0.01) # run anova, specifying a different alpha level if desired anova(resi.obj, alpha = 0.05)# fit a model mod = lm(charges ~ bmi + sex, data = RESI::insurance) # run resi with the store.boot = TRUE option resi.obj = resi(mod, nboot = 100, store.boot = TRUE, alpha = 0.01) # run anova, specifying a different alpha level if desired anova(resi.obj, alpha = 0.05)
This function computes the robust effect size index from Vandekar, Tao, & Blume (2020). Vector arguments are accepted. If different length arguments are passed they are dealt with in the usual way of R. For mixed effects models, RESI is conditional on the average correlation structure within subjects.
chisq2S(chisq, df, n)chisq2S(chisq, df, n)
chisq |
The chi-square statistic for the parameter of interest. |
df |
Number of degrees of freedom of the chi-square statistic. |
n |
Number of independent samples. |
The formula for converting a Chi-square statistic to RESI is:
Returns a scalar or vector argument of the the robust effect size index estimate.
# obtain Chi-sq value by fitting an lm and running a Wald test mod = lm(charges ~ region * age + bmi + sex, data = RESI::insurance) # run a Wald test with robust variance wt = lmtest::waldtest(mod, vcov = sandwich::vcovHC, test = "Chisq") # get Chi-sq value and degrees of freedom chisq = wt$Chisq[2] df = abs(wt$Df[2]) # run chisq2S to convert to RESI chisq2S(chisq, df = df, n = nrow(mod$model))# obtain Chi-sq value by fitting an lm and running a Wald test mod = lm(charges ~ region * age + bmi + sex, data = RESI::insurance) # run a Wald test with robust variance wt = lmtest::waldtest(mod, vcov = sandwich::vcovHC, test = "Chisq") # get Chi-sq value and degrees of freedom chisq = wt$Chisq[2] df = abs(wt$Df[2]) # run chisq2S to convert to RESI chisq2S(chisq, df = df, n = nrow(mod$model))
Converts Cohen's d robust effect size index (S) using the formula from Vandekar, Tao, & Blume (2020).
d2S(d, pi = 0.5)d2S(d, pi = 0.5)
d |
Numeric, value of Cohen's d. |
pi |
Numeric, the sampling proportions. |
The pi parameter comes from the fact that Cohen's d doesn't account for unequal sample proportions in the population, but S does.
The default is set to a natural value 1/2, which corresponds to a case control design, for example, where sampling proportions always are controlled by the experimenter.
The formula to convert Cohen's d to S is:
Returns an estimate the robust effect size index
# Consider an experiment with equal sampling proportions and a medium effect size # corresponding to a Cohen's d of 0.5. # convert to RESI (S) d2S(d = 0.5) # This corresponds to a RESI of 0.25.# Consider an experiment with equal sampling proportions and a medium effect size # corresponding to a Cohen's d of 0.5. # convert to RESI (S) d2S(d = 0.5) # This corresponds to a RESI of 0.25.
A longitudinal dataset comparing two treatments for depression.
depressiondepression
A data frame with 1020 rows and 5 variables:
diagnosed depression severity
treatment; standard or new
patient id
time point of treatment
depression response at time of treatment. 1 = Normal, 0 = Abnormal
http://static.lib.virginia.edu/statlab/materials/data/depression.csv
Agresti, A. (2002). Categorical Data Analysis. Wiley, 2nd Edition.
This function computes the robust effect size index from Vandekar, Tao, & Blume (2020). Vector arguments are accepted. If different length arguments are passed they are dealt with in the usual way of R.
f2S(f, df, rdf, n)f2S(f, df, rdf, n)
f |
The F statistic for the parameter of interest. |
df |
Number of degrees of freedom of the F statistic. |
rdf |
Model residual degrees of freedom. |
n |
Number of independent samples. |
The formula for converting an F statistic to S is:
The estimator is derived by setting the statistic equal to the expected value of the test statistic and solving for S.
Returns a scalar or vector argument of the the robust effect size index estimate.
# to obtain example F values, first fit a lm mod = lm(charges ~ region * age + bmi + sex, data = RESI::insurance) # run Anova, using a robust variance-covariance function # get the F values and Df values fs = car::Anova(mod, vcov. = sandwich::vcovHC)[1:5, "F"] dfs = car::Anova(mod, vcov. = sandwich::vcovHC)[1:5, "Df"] # get RESI estimates f2S(fs, df = dfs, rdf = mod$df.residual, n = nrow(RESI::insurance))# to obtain example F values, first fit a lm mod = lm(charges ~ region * age + bmi + sex, data = RESI::insurance) # run Anova, using a robust variance-covariance function # get the F values and Df values fs = car::Anova(mod, vcov. = sandwich::vcovHC)[1:5, "F"] dfs = car::Anova(mod, vcov. = sandwich::vcovHC)[1:5, "Df"] # get RESI estimates f2S(fs, df = dfs, rdf = mod$df.residual, n = nrow(RESI::insurance))
Converts Cohen's f^2 to robust effect size index (S) using the formula from Vandekar, Tao, & Blume (2020).
fsq2S(fsq)fsq2S(fsq)
fsq |
Numeric, value of Cohen's f^2. |
The formula for the conversion is:
Returns an estimate the robust effect size index
# consider a moderate effect size of f^2 = 0.3 fsq2S(0.3) # This corresponds to a RESI of 0.5477226# consider a moderate effect size of f^2 = 0.3 fsq2S(0.3) # This corresponds to a RESI of 0.5477226
This function uses ggplot2 graphics to plot robust effect size (RESI) estimates and confidence intervals from 'resi', 'summary_resi', and 'anova_resi' objects.
## S3 method for class 'resi' ggplot(data, mapping, alpha = NULL, error.bars = TRUE, ..., environment)## S3 method for class 'resi' ggplot(data, mapping, alpha = NULL, error.bars = TRUE, ..., environment)
data |
Object of 'resi', 'summary_resi', or 'anova_resi' class |
mapping |
Ignored, included for consistency with 'ggplot' generic |
alpha |
Numeric, desired alpha level for confidence intervals |
error.bars |
Logical, whether to include end caps on the confidence intervals. Default = 'TRUE' |
... |
Ignored |
environment |
Ignored, included for consistency with 'ggplot' generic |
Returns a ggplot of RESI point estimates
# create a resi object resi_obj <- resi(lm(charges ~ region * age + bmi + sex, data = RESI::insurance), nboot = 10) # plot ANOVA table ggplot2::ggplot(anova(resi_obj))# create a resi object resi_obj <- resi(lm(charges ~ region * age + bmi + sex, data = RESI::insurance), nboot = 10) # plot ANOVA table ggplot2::ggplot(anova(resi_obj))
A dataset with 1338 observations on health insurance charges and demographic factors.
insuranceinsurance
A data frame with 1338 rows and 7 variables:
age of primary beneficiary in years
insurance contractor sex, male/female
body mass index
number of dependents
smoker/non-smoker
beneficiary's region of US
individual medical costs billed by health insurance
https://www.kaggle.com/datasets/teertha/ushealthinsurancedataset
After running the resi function on a fitted model, this function
can be used to print the overall Wald test component. If the resi function was run
with the 'store.boot = TRUE' option to store the full matrix of bootstrapped estimates,
the user can specify a different alpha level for this function's confidence intervals.
omnibus(object, alpha = NULL, ...)omnibus(object, alpha = NULL, ...)
object |
an object resulting from resi function |
alpha |
an optional new specification for the confidence level. Can be vector-valued |
... |
ignored |
Returns a 'omnibus_resi' object containing the computed omnibus Wald test
# fit a model mod = lm(charges ~ bmi + sex, data = RESI::insurance) # run resi with the store.boot = TRUE option resi_obj = resi(mod, nboot = 100, store.boot = TRUE, alpha = 0.01) # run summary, specifying a different alpha level if desired omnibus(resi_obj, alpha = 0.05)# fit a model mod = lm(charges ~ bmi + sex, data = RESI::insurance) # run resi with the store.boot = TRUE option resi_obj = resi(mod, nboot = 100, store.boot = TRUE, alpha = 0.01) # run summary, specifying a different alpha level if desired omnibus(resi_obj, alpha = 0.05)
This function uses base graphics to plot robust effect size (RESI) estimates and confidence intervals from 'resi', 'summary_resi', and 'anova_resi' objects.
## S3 method for class 'resi' plot( x, alpha = NULL, ycex.axis = NULL, yaxis.args = list(), automar = TRUE, ... )## S3 method for class 'resi' plot( x, alpha = NULL, ycex.axis = NULL, yaxis.args = list(), automar = TRUE, ... )
x |
Object of 'resi', 'summary_resi', or 'anova_resi' class |
alpha |
Numeric, desired alpha level for confidence intervals |
ycex.axis |
Numeric, scale specifically for the variable name labels |
yaxis.args |
List, other arguments to be passed to |
automar |
Logical, whether to automatically adjust the plotting margins to accommodate variable names. Default = 'TRUE' |
... |
This function creates a forest-like plot with RESI estimates for each variable or factor. The size of the left margin will be automatically adjusted (and returned to original after plotting) unless 'automar = FALSE'. Additional graphics parameters will be passed to the main plot function, the confidence intervals. Arguments specifically for the y-axis (variable names) can be specified using 'yaxis.args'. To manually adjust the size of the y-axis labels without affecting the x-axis, the user can specify a value for 'ycex.axis'.
Returns a plot of RESI point estimates
# create a resi object resi_obj <- resi(lm(charges ~ region * age + bmi + sex, data = RESI::insurance), nboot = 10) # plot coefficients table, changing size of labels for both axes in the usual way plot(resi_obj, cex.axis = 0.7) # plot ANOVA table, changing the size of just the y-axis plot(resi_obj, ycex.axis = 0.8)# create a resi object resi_obj <- resi(lm(charges ~ region * age + bmi + sex, data = RESI::insurance), nboot = 10) # plot coefficients table, changing size of labels for both axes in the usual way plot(resi_obj, cex.axis = 0.7) # plot ANOVA table, changing the size of just the y-axis plot(resi_obj, ycex.axis = 0.8)
This function will estimate the robust effect size (RESI) from Vandekar, Tao, & Blume (2020) and its confidence interval in various ways for a fitted model object. The overall RESI is estimated via a Wald test. RESI is (optionally) estimated for each factor in coefficients-style table. RESI is (optionally) estimated for each variable/interaction in an Anova-style table for models with existing Anova methods. CIs can be calculated using either non-parametric or Bayesian bootstrapping.
resi(model.full, ...) ## Default S3 method: resi( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, nboot = 1000, boot.method = "nonparam", vcovfunc = sandwich::vcovHC, alpha = 0.05, store.boot = FALSE, Anova.args = list(), vcov.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), long = FALSE, clvar = NULL, ... ) ## S3 method for class 'glm' resi( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, nboot = 1000, vcovfunc = sandwich::vcovHC, alpha = 0.05, store.boot = FALSE, Anova.args = list(), vcov.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'lm' resi( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, nboot = 1000, boot.method = "nonparam", vcovfunc = sandwich::vcovHC, alpha = 0.05, store.boot = FALSE, Anova.args = list(), vcov.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'nls' resi( model.full, model.reduced = NULL, data, coefficients = TRUE, overall = TRUE, nboot = 1000, boot.method = "nonparam", anova = FALSE, vcovfunc = r_nlshc, alpha = 0.05, store.boot = FALSE, vcov.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'survreg' resi( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, nboot = 1000, vcovfunc = vcov, alpha = 0.05, store.boot = FALSE, Anova.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'coxph' resi( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, nboot = 1000, vcovfunc = vcov, alpha = 0.05, store.boot = FALSE, Anova.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'hurdle' resi( model.full, model.reduced = NULL, data, coefficients = TRUE, overall = TRUE, nboot = 1000, vcovfunc = sandwich::sandwich, anova = FALSE, alpha = 0.05, store.boot = FALSE, vcov.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'zeroinfl' resi( model.full, model.reduced = NULL, data, coefficients = TRUE, overall = TRUE, nboot = 1000, vcovfunc = sandwich::sandwich, anova = FALSE, alpha = 0.05, store.boot = FALSE, vcov.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'geeglm' resi( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, nboot = 1000, alpha = 0.05, store.boot = FALSE, unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'gee' resi( model.full, data, nboot = 1000, alpha = 0.05, store.boot = FALSE, unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'lme' resi( model.full, alpha = 0.05, nboot = 1000, anova = TRUE, vcovfunc = clubSandwich::vcovCR, vcov.args = list(), ... ) ## S3 method for class 'lmerMod' resi( model.full, alpha = 0.05, nboot = 1000, anova = TRUE, vcovfunc = clubSandwich::vcovCR, vcov.args = list(), ... )resi(model.full, ...) ## Default S3 method: resi( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, nboot = 1000, boot.method = "nonparam", vcovfunc = sandwich::vcovHC, alpha = 0.05, store.boot = FALSE, Anova.args = list(), vcov.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), long = FALSE, clvar = NULL, ... ) ## S3 method for class 'glm' resi( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, nboot = 1000, vcovfunc = sandwich::vcovHC, alpha = 0.05, store.boot = FALSE, Anova.args = list(), vcov.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'lm' resi( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, nboot = 1000, boot.method = "nonparam", vcovfunc = sandwich::vcovHC, alpha = 0.05, store.boot = FALSE, Anova.args = list(), vcov.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'nls' resi( model.full, model.reduced = NULL, data, coefficients = TRUE, overall = TRUE, nboot = 1000, boot.method = "nonparam", anova = FALSE, vcovfunc = r_nlshc, alpha = 0.05, store.boot = FALSE, vcov.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'survreg' resi( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, nboot = 1000, vcovfunc = vcov, alpha = 0.05, store.boot = FALSE, Anova.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'coxph' resi( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, nboot = 1000, vcovfunc = vcov, alpha = 0.05, store.boot = FALSE, Anova.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'hurdle' resi( model.full, model.reduced = NULL, data, coefficients = TRUE, overall = TRUE, nboot = 1000, vcovfunc = sandwich::sandwich, anova = FALSE, alpha = 0.05, store.boot = FALSE, vcov.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'zeroinfl' resi( model.full, model.reduced = NULL, data, coefficients = TRUE, overall = TRUE, nboot = 1000, vcovfunc = sandwich::sandwich, anova = FALSE, alpha = 0.05, store.boot = FALSE, vcov.args = list(), unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'geeglm' resi( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, nboot = 1000, alpha = 0.05, store.boot = FALSE, unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'gee' resi( model.full, data, nboot = 1000, alpha = 0.05, store.boot = FALSE, unbiased = TRUE, parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", 1L), ... ) ## S3 method for class 'lme' resi( model.full, alpha = 0.05, nboot = 1000, anova = TRUE, vcovfunc = clubSandwich::vcovCR, vcov.args = list(), ... ) ## S3 method for class 'lmerMod' resi( model.full, alpha = 0.05, nboot = 1000, anova = TRUE, vcovfunc = clubSandwich::vcovCR, vcov.args = list(), ... )
model.full |
|
... |
Ignored. |
model.reduced |
Fitted model object of same type as model.full. By default 'NULL'; the same model as the full model but only having intercept. |
data |
Data.frame or object coercible to data.frame of model.full data (required for some model types). |
anova |
Logical, whether to produce an Anova table with the RESI columns added. By default = 'TRUE'. |
coefficients |
Logical, whether to produce a coefficients (summary) table with the RESI columns added. By default = 'TRUE'. |
overall |
Logical, whether to produce an overall Wald test comparing full to reduced model with RESI columns added. By default = 'TRUE'. |
nboot |
Numeric, the number of bootstrap replicates. By default, 1000. |
boot.method |
String, which type of bootstrap to use: 'nonparam' = non-parametric bootstrap (default); 'bayes' = Bayesian bootstrap. |
vcovfunc |
The variance estimator function for constructing the Wald test statistic. By default, sandwich::vcovHC (the robust (sandwich) variance estimator). |
alpha |
Numeric, significance level of the constructed CIs. By default, 0.05. |
store.boot |
Logical, whether to store all the bootstrapped estimates. By default, 'FALSE'. |
Anova.args |
List, additional arguments to be passed to Anova function. |
vcov.args |
List, additional arguments to be passed to vcovfunc. |
unbiased |
Logical, whether to use the unbiased or alternative T/Z statistic to RESI conversion. By default, 'TRUE'. See details. |
parallel |
See documentation for |
ncpus |
See documentation for |
long |
Logical, whether the data is longitudinal/clustered. By default, 'FALSE'. |
clvar |
Character, the name of the cluster/id variable if data is clustered. By default, 'NULL'. |
The RESI, denoted as S, is applicable across many model types. It is a unitless
index and can be easily be compared across models. The RESI can also be
converted to Cohen's d (S2d) under model homoskedasticity.
This function computes the RESI point estimates and bootstrapped confidence
intervals based on Chi-square, F, T, or Z statistics. The robust (sandwich)
variance is used by default, allowing for consistency under
model-misspecification. The RESI is related to the non-centrality parameter
of the test statistic. The RESI estimate is consistent for all four
(Chi-square, F, T, and Z) types of statistics used. The Chi-square and F-based
calculations rely on asymptotic theory, so they may be biased in small samples.
When possible, the T and Z statistics are used. There are two formulas for both
the T and Z statistic conversion. The first (default, unbiased = TRUE)
are based on solving the expected value of the T or Z statistic for the RESI.
The alternative is based on squaring the T or Z statistic and using the
F or Chi-square statistic conversion. Both of these methods are consistent, but
the alternative exhibits a notable amount of finite sample bias. The alternative
may be appealing because its absolute value will be equal to the RESI based on
the F or Chi-square statistic. The RESI based on the Chi-Square and F statistics
is always greater than or equal to 0. The type of statistic
used is listed with the output. See f2S, chisq2S,
t2S, and z2S for more details on the formulas.
For GEE (geeglm) models, a longitudinal RESI (L-RESI) and a cross-sectional,
per-measurement RESI (CS-RESI) is estimated. The longitudinal RESI takes the
specified clustering into account, while the cross-sectional RESI is estimated
using a model where each measurement is its own cluster.
For most lm and nls model types, there is a Bayesian bootstrap
option available as an alternative to the default, standard non-parametric
bootstrap. The interpretation of a Bayesian bootstrapped interval is similar to
that of a credible interval.
Certain model types require the data used for the model be entered as an argument.
These are: nls, survreg, and coxph. Additionally, if a model
includes certain functions (splines, factor, I), the data needs to be provided.
If running into convergence issues with nls models, it is advised to refit the
nls model with starting values equal to the estimates provided by the model
and then try rerunning resi.
Returns a list that includes function arguments, RESI point estimates, and confidence intervals in coefficients/anova-style tables
resi(default): RESI point and interval estimation for models
resi(glm): RESI point and interval estimation for models
resi(lm): RESI point and interval estimation for lm models
resi(nls): RESI point and interval estimation for nls models
resi(survreg): RESI point and interval estimation for survreg models
resi(coxph): RESI point and interval estimation for coxph models
resi(hurdle): RESI point and interval estimation for hurdle models
resi(zeroinfl): RESI point and interval estimation for zeroinfl models
resi(geeglm): RESI point and interval estimation for GEE models
resi(gee): RESI point and interval estimation for GEE models
resi(lme): RESI point and interval estimation for LME (nlme) models
resi(lmerMod): RESI point and interval estimation for lmerMod models
Vandekar S, Tao R, Blume J. A Robust Effect Size Index. Psychometrika. 2020 Mar;85(1):232-246. doi: 10.1007/s11336-020-09698-2.
Kang, K., Armstrong, K., Avery, S., McHugo, M., Heckers, S., & Vandekar, S. (2021). Accurate confidence interval estimation for non-centrality parameters and effect size indices. arXiv preprint arXiv:2111.05966.
resi_pe, vcovHC,
f2S, chisq2S, z2S, t2S
## for timing purposes, a small number of bootstrap replicates is used in the ## examples. Run them with a higher or default `nboot` argument for better performance ## RESI on a linear model # fit linear model mod = lm(charges ~ region * age + bmi + sex, data = RESI::insurance) # run resi on fitted model with desired number of bootstrap replicates # store bootstrap results for calculating different CIs later resi_obj = resi(mod, nboot = 50, store.boot = TRUE) # print output resi_obj # fit a reduced model for comparison mod_red = lm(charges ~ bmi, data = RESI::insurance) # running resi and including the reduced model will provide almost the exact same # output as not including a reduced model. The difference is that the "overall" # portion of the output will compare the full model to the reduced model. # The "summary" and "anova" RESI estimates will be the same. (The bootstrapped # confidence intervals may differ.) resi(model.full = mod, model.reduced = mod_red, nboot = 10) # used stored bootstrap results to get a different alpha-level confidence interval summary(resi_obj, alpha = c(0.01, 0.1)) car::Anova(resi_obj, alpha = c(0.01, 0.1)) # the result of resi, as well as the summary or Anova of a `resi` object can be plotted # if the resi object was created with the store.boot = `TRUE` option, any alpha # can be specified plot(resi_obj, alpha = 0.01) # if the variable names on the y-axis are too long, you can reduce their size with # the ycex.axis argument (or use regular common solutions like changing the margins) plot(resi_obj, alpha = 0.01, ycex.axis = 0.5) # for some model types and formula structures, data argument is required if(requireNamespace("splines")){ # fit logistic regression model with splines mod = glm(smoker ~ splines::ns(age, df = 3) + region, data = RESI::insurance, family = "binomial") # specify additional arguments to the variance-covariance function via vcov.args resi_obj = resi(mod, data = RESI::insurance, alpha = 0.01, vcov.args = list(type = "HC0"), nboot = 25) summary(resi_obj) car::Anova(resi_obj)} ## RESI on a survival model with alternate Z2S if(requireNamespace("survival")){ # fit coxph model on example data from survival package # Note: for survival models, you need to specify robust variance in the model # creation. resi will ignore the vcovfunc argument for this reason. mod.coxph = survival::coxph(survival::Surv(time, status) ~ age + sex + wt.loss, data=survival::lung, robust = TRUE) # run resi on the model # to use the alternative Z to RESI formula (which is equal in absolute value to the # chi-square to RESI (S) formula), specify unbiased = FALSE. resi(mod.coxph, data = survival::lung, unbiased = FALSE, nboot = 10)}## for timing purposes, a small number of bootstrap replicates is used in the ## examples. Run them with a higher or default `nboot` argument for better performance ## RESI on a linear model # fit linear model mod = lm(charges ~ region * age + bmi + sex, data = RESI::insurance) # run resi on fitted model with desired number of bootstrap replicates # store bootstrap results for calculating different CIs later resi_obj = resi(mod, nboot = 50, store.boot = TRUE) # print output resi_obj # fit a reduced model for comparison mod_red = lm(charges ~ bmi, data = RESI::insurance) # running resi and including the reduced model will provide almost the exact same # output as not including a reduced model. The difference is that the "overall" # portion of the output will compare the full model to the reduced model. # The "summary" and "anova" RESI estimates will be the same. (The bootstrapped # confidence intervals may differ.) resi(model.full = mod, model.reduced = mod_red, nboot = 10) # used stored bootstrap results to get a different alpha-level confidence interval summary(resi_obj, alpha = c(0.01, 0.1)) car::Anova(resi_obj, alpha = c(0.01, 0.1)) # the result of resi, as well as the summary or Anova of a `resi` object can be plotted # if the resi object was created with the store.boot = `TRUE` option, any alpha # can be specified plot(resi_obj, alpha = 0.01) # if the variable names on the y-axis are too long, you can reduce their size with # the ycex.axis argument (or use regular common solutions like changing the margins) plot(resi_obj, alpha = 0.01, ycex.axis = 0.5) # for some model types and formula structures, data argument is required if(requireNamespace("splines")){ # fit logistic regression model with splines mod = glm(smoker ~ splines::ns(age, df = 3) + region, data = RESI::insurance, family = "binomial") # specify additional arguments to the variance-covariance function via vcov.args resi_obj = resi(mod, data = RESI::insurance, alpha = 0.01, vcov.args = list(type = "HC0"), nboot = 25) summary(resi_obj) car::Anova(resi_obj)} ## RESI on a survival model with alternate Z2S if(requireNamespace("survival")){ # fit coxph model on example data from survival package # Note: for survival models, you need to specify robust variance in the model # creation. resi will ignore the vcovfunc argument for this reason. mod.coxph = survival::coxph(survival::Surv(time, status) ~ age + sex + wt.loss, data=survival::lung, robust = TRUE) # run resi on the model # to use the alternative Z to RESI formula (which is equal in absolute value to the # chi-square to RESI (S) formula), specify unbiased = FALSE. resi(mod.coxph, data = survival::lung, unbiased = FALSE, nboot = 10)}
This function will estimate the robust effect size (RESI) from Vandekar, Tao, & Blume (2020).
The overall RESI is estimated via a Wald test. RESI is (optionally) estimated for each factor in coefficients-style table.
RESI is (optionally) estimated for each variable/interaction in an Anova-style table
for models with existing Anova methods. This function is the building block for the resi function.
resi_pe(model.full, ...) ## Default S3 method: resi_pe( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, vcovfunc = sandwich::vcovHC, Anova.args = list(), vcov.args = list(), unbiased = TRUE, waldtype = 0, ... ) ## S3 method for class 'glm' resi_pe( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, vcovfunc = sandwich::vcovHC, Anova.args = list(), vcov.args = list(), unbiased = TRUE, waldtype = 0, ... ) ## S3 method for class 'lm' resi_pe( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, vcovfunc = sandwich::vcovHC, Anova.args = list(), vcov.args = list(), unbiased = TRUE, overall = TRUE, ... ) ## S3 method for class 'nls' resi_pe( model.full, model.reduced = NULL, data, coefficients = TRUE, anova = FALSE, vcovfunc = r_nlshc, vcov.args = list(), unbiased = TRUE, overall = TRUE, ... ) ## S3 method for class 'survreg' resi_pe( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, vcovfunc = vcov, Anova.args = list(), unbiased = TRUE, overall = TRUE, ... ) ## S3 method for class 'coxph' resi_pe( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, vcovfunc = vcov, Anova.args = list(), unbiased = TRUE, overall = TRUE, ... ) ## S3 method for class 'hurdle' resi_pe( model.full, model.reduced = NULL, data, coefficients = TRUE, anova = TRUE, vcovfunc = sandwich::sandwich, vcov.args = list(), unbiased = TRUE, overall = TRUE, ... ) ## S3 method for class 'zeroinfl' resi_pe( model.full, model.reduced = NULL, data, coefficients = TRUE, anova = TRUE, vcovfunc = sandwich::sandwich, vcov.args = list(), unbiased = TRUE, overall = TRUE, ... ) ## S3 method for class 'geeglm' resi_pe( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, unbiased = TRUE, ... ) ## S3 method for class 'gee' resi_pe(model.full, data, unbiased = TRUE, ...) ## S3 method for class 'lme' resi_pe( model.full, anova = TRUE, vcovfunc = clubSandwich::vcovCR, Anova.args = list(), vcov.args = list(), ... ) ## S3 method for class 'lmerMod' resi_pe( model.full, anova = TRUE, vcovfunc = clubSandwich::vcovCR, Anova.args = list(), vcov.args = list(), ... )resi_pe(model.full, ...) ## Default S3 method: resi_pe( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, vcovfunc = sandwich::vcovHC, Anova.args = list(), vcov.args = list(), unbiased = TRUE, waldtype = 0, ... ) ## S3 method for class 'glm' resi_pe( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, vcovfunc = sandwich::vcovHC, Anova.args = list(), vcov.args = list(), unbiased = TRUE, waldtype = 0, ... ) ## S3 method for class 'lm' resi_pe( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, vcovfunc = sandwich::vcovHC, Anova.args = list(), vcov.args = list(), unbiased = TRUE, overall = TRUE, ... ) ## S3 method for class 'nls' resi_pe( model.full, model.reduced = NULL, data, coefficients = TRUE, anova = FALSE, vcovfunc = r_nlshc, vcov.args = list(), unbiased = TRUE, overall = TRUE, ... ) ## S3 method for class 'survreg' resi_pe( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, vcovfunc = vcov, Anova.args = list(), unbiased = TRUE, overall = TRUE, ... ) ## S3 method for class 'coxph' resi_pe( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, vcovfunc = vcov, Anova.args = list(), unbiased = TRUE, overall = TRUE, ... ) ## S3 method for class 'hurdle' resi_pe( model.full, model.reduced = NULL, data, coefficients = TRUE, anova = TRUE, vcovfunc = sandwich::sandwich, vcov.args = list(), unbiased = TRUE, overall = TRUE, ... ) ## S3 method for class 'zeroinfl' resi_pe( model.full, model.reduced = NULL, data, coefficients = TRUE, anova = TRUE, vcovfunc = sandwich::sandwich, vcov.args = list(), unbiased = TRUE, overall = TRUE, ... ) ## S3 method for class 'geeglm' resi_pe( model.full, model.reduced = NULL, data, anova = TRUE, coefficients = TRUE, overall = TRUE, unbiased = TRUE, ... ) ## S3 method for class 'gee' resi_pe(model.full, data, unbiased = TRUE, ...) ## S3 method for class 'lme' resi_pe( model.full, anova = TRUE, vcovfunc = clubSandwich::vcovCR, Anova.args = list(), vcov.args = list(), ... ) ## S3 method for class 'lmerMod' resi_pe( model.full, anova = TRUE, vcovfunc = clubSandwich::vcovCR, Anova.args = list(), vcov.args = list(), ... )
model.full |
|
... |
Ignored. |
model.reduced |
Fitted model object of same type as model.full. By default 'NULL'; the same model as the full model but only having intercept. |
data |
Data.frame or object coercible to data.frame of model.full data (required for some model types). |
anova |
Logical, whether to produce an Anova table with the RESI columns added. By default = 'TRUE'. |
coefficients |
Logical, whether to produce a coefficients (summary) table with the RESI columns added. By default = 'TRUE'. |
overall |
Logical, whether to produce an overall Wald test comparing full to reduced model with RESI columns added. By default = 'TRUE'. |
vcovfunc |
The variance estimator function for constructing the Wald test statistic. By default, sandwich::vcovHC (the robust (sandwich) variance estimator). |
Anova.args |
List, additional arguments to be passed to Anova function. |
vcov.args |
List, additional arguments to be passed to vcovfunc. |
unbiased |
Logical, whether to use the unbiased or alternative T/Z statistic to RESI conversion. By default, 'TRUE'. See details. |
waldtype |
Numeric, indicates which function to use for overall Wald test. 0 (default) = lmtest::waldtest Chi-square, 1 = lmtest::waldtest F, 2 = aod::wald.test |
The Robust Effect Size Index (RESI) is an effect size measure based on M-estimators.
This function is called by resi a specified number of times to
form bootstrapped confidence intervals. Called by itself, this function will
only calculate point estimates.
The RESI, denoted as S, is applicable across many model types. It is a unitless
index and can be easily be compared across models. The RESI can also be
converted to Cohen's d (S2d) under model homoskedasticity.
The RESI is related to the non-centrality parameter
of the test statistic. The RESI estimate is consistent for all four
(Chi-square, F, T, and Z) types of statistics used. The Chi-square and F-based
calculations rely on asymptotic theory, so they may be biased in small samples.
When possible, the T and Z statistics are used. There are two formulas for both
the T and Z statistic conversion. The first (default, unbiased = TRUE)
are based on solving the expected value of the T or Z statistic for the RESI.
The alternative is based on squaring the T or Z statistic and using the
F or Chi-square statistic conversion. Both of these methods are consistent, but
the alternative exhibits a notable amount of finite sample bias. The alternative
may be appealing because its absolute value will be equal to the RESI based on
the F or Chi-square statistic. The RESI based on the Chi-Square and F statistics
is always greater than or equal to 0. The type of statistic
used is listed with the output. See f2S, chisq2S,
t2S, and z2S for more details on the formulas.
For GEE (geeglm) models, a longitudinal RESI (L-RESI) and a cross-sectional,
per-measurement RESI (CS-RESI) is estimated. The longitudinal RESI takes the
specified clustering into account, while the cross-sectional RESI is estimated
using a model where each measurement is its own cluster.
Returns a list containing RESI point estimates
resi_pe(default): RESI point estimation
resi_pe(glm): RESI point estimation for generalized linear models
resi_pe(lm): RESI point estimation for linear models
resi_pe(nls): RESI point estimation for nonlinear least squares models
resi_pe(survreg): RESI point estimation for survreg
resi_pe(coxph): RESI point estimation for coxph models
resi_pe(hurdle): RESI point estimation for hurdle models
resi_pe(zeroinfl): RESI point estimation for zeroinfl models
resi_pe(geeglm): RESI point estimation for geeglm object
resi_pe(gee): RESI point estimation for gee object
resi_pe(lme): RESI point estimation for lme object
resi_pe(lmerMod): RESI point estimation for lmerMod object
Vandekar S, Tao R, Blume J. A Robust Effect Size Index. Psychometrika. 2020 Mar;85(1):232-246. doi: 10.1007/s11336-020-09698-2.
# This function produces point estimates for the RESI. The resi function will # provide the same point estimates but adds confidence intervals. See resi for # more detailed examples. ## resi_pe for a linear model # fit linear model mod <- lm(charges ~ region * age + bmi + sex, data = RESI::insurance) # run resi_pe on the model resi_pe(mod) # if you want to have RESI estimates in the coefficient table that are equal in absolute # value to those in the Anova table (except for those with >1 df and/or included in other # interaction terms), you can specify unbiased = FALSE to use the alternate conversion. resi_pe(mod, unbiased = FALSE)# This function produces point estimates for the RESI. The resi function will # provide the same point estimates but adds confidence intervals. See resi for # more detailed examples. ## resi_pe for a linear model # fit linear model mod <- lm(charges ~ region * age + bmi + sex, data = RESI::insurance) # run resi_pe on the model resi_pe(mod) # if you want to have RESI estimates in the coefficient table that are equal in absolute # value to those in the Anova table (except for those with >1 df and/or included in other # interaction terms), you can specify unbiased = FALSE to use the alternate conversion. resi_pe(mod, unbiased = FALSE)
Converts R^2, the partial coefficient of determination, to robust effect size index (S) using the formula from Vandekar, Tao, & Blume (2020).
Rsq2S(Rsq)Rsq2S(Rsq)
Rsq |
Numeric, R^2 |
The formula for the conversion is:
Returns an estimate of R^2 based on the RESI
# consider a moderate effect size of R^2 = 0.1 Rsq2S(0.1) # this corresponds to a RESI of 0.333# consider a moderate effect size of R^2 = 0.1 Rsq2S(0.1) # this corresponds to a RESI of 0.333
Converts the robust effect size index (S) to Cohen's d using the formula from Vandekar, Tao, & Blume (2020).
S2d(S, pi = 0.5)S2d(S, pi = 0.5)
S |
Numeric, the robust effect size index. |
pi |
Numeric, the sampling proportions. |
The pi parameter comes from the fact that Cohen's d doesn't account for unequal sample proportions in the population, but S does.
The default is set to a natural value 1/2, which corresponds to a case control design, for example, where sampling proportions always are controlled by the experimenter.
The formula for the conversion is:
Returns an estimate of Cohen's d based on the RESI.
# fit a simple linear regression with a binary predictor mod = lm(charges ~ sex, data = RESI::insurance) # calculate t-value t = summary(mod)$coefficients[2, "t value"] # calculate RESI (S) S = t2S(t, n = 1338, rdf = 1336) # determine sample proportions pi = length(which(RESI::insurance[,"sex"]=="male"))/1338 # convert S to Cohen's d S2d(S = S, pi = pi)# fit a simple linear regression with a binary predictor mod = lm(charges ~ sex, data = RESI::insurance) # calculate t-value t = summary(mod)$coefficients[2, "t value"] # calculate RESI (S) S = t2S(t, n = 1338, rdf = 1336) # determine sample proportions pi = length(which(RESI::insurance[,"sex"]=="male"))/1338 # convert S to Cohen's d S2d(S = S, pi = pi)
Converts robust effect size index (S) to Cohen's f^2 (effect size for multiple regression) using the formula from Vandekar, Tao, & Blume (2020).
S2fsq(S)S2fsq(S)
S |
Numeric,the robust effect size index. |
The formula for the conversion is:
Returns an estimate of Cohen's f^2 based on the RESI
# fit a linear regression model with continuous outcome and predictor mod = lm(charges ~ age, data = RESI::insurance) # obtain t value for calculating RESI t = summary(mod)$coefficients[2, "t value"] # calculate RESI S = t2S(t, n = 1338, rdf = 1336) # convert to f^2 S2fsq(S)# fit a linear regression model with continuous outcome and predictor mod = lm(charges ~ age, data = RESI::insurance) # obtain t value for calculating RESI t = summary(mod)$coefficients[2, "t value"] # calculate RESI S = t2S(t, n = 1338, rdf = 1336) # convert to f^2 S2fsq(S)
Converts robust effect size index (S) to R^2, the partial coefficient of determination, using the formula from Vandekar, Tao, & Blume (2020).
S2Rsq(S)S2Rsq(S)
S |
Numeric, the robust effect size index. |
The formula for the conversion is:
Returns an estimate of R^2 based on the RESI
# fit a simple linear regression with a binary predictor mod = lm(charges ~ sex, data = RESI::insurance) # calculate t-value t = summary(mod)$coefficients[2, "t value"] # calculate RESI (S) S = t2S(t, n = 1338, rdf = 1336) # convert S to R^2 S2Rsq(S)# fit a simple linear regression with a binary predictor mod = lm(charges ~ sex, data = RESI::insurance) # calculate t-value t = summary(mod)$coefficients[2, "t value"] # calculate RESI (S) S = t2S(t, n = 1338, rdf = 1336) # convert S to R^2 S2Rsq(S)
After running the resi function on a fitted model, this function can
be used to print the coefficients table component. If the resi function was run with
the 'store.boot = TRUE' option to store the full matrix of bootstrapped estimates,
the user can specify a different alpha level for this function's confidence intervals.
## S3 method for class 'resi' summary(object, alpha = NULL, ...)## S3 method for class 'resi' summary(object, alpha = NULL, ...)
object |
an object resulting from resi function |
alpha |
an optional new specification for the confidence level. Can be vector-valued |
... |
ignored |
Returns a 'summary_resi' object containing the computed coefficients table
# fit a model mod = lm(charges ~ bmi + sex, data = RESI::insurance) # run resi with the store.boot = TRUE option resi_obj = resi(mod, nboot = 100, store.boot = TRUE, alpha = 0.01) # run summary, specifying a different alpha level if desired summary(resi_obj, alpha = 0.05)# fit a model mod = lm(charges ~ bmi + sex, data = RESI::insurance) # run resi with the store.boot = TRUE option resi_obj = resi(mod, nboot = 100, store.boot = TRUE, alpha = 0.01) # run summary, specifying a different alpha level if desired summary(resi_obj, alpha = 0.05)
This function computes the robust effect size index from Vandekar, Tao, & Blume (2020). Vector arguments are accepted. If different length arguments are passed they are dealt with in the usual way of R.
t2S(t, rdf, n, unbiased = TRUE)t2S(t, rdf, n, unbiased = TRUE)
t |
The t statistic for the parameter of interest. |
rdf |
Model residual degrees of freedom/degrees of freedom of the t statistic. |
n |
Number of independent samples. |
unbiased |
Logical, whether to use unbiased or alternative estimator. See details. |
This function computes S, the RESI, from a t statistic. The formula for the unbiased estimator (default) is derived by solving the expected value of the t statistic for S. It is unbiased and consistent.
The formula for the unbiased conversion is:
The formula for the alternative estimator is derived by squaring the t statistic
and using the f2S formula. This estimator may be appealing for its
intuitive relationship to the F statistic; the absolute value of RESI estimates
using this formula will be equal to a RESI estimate using an F statistic for
the same model. However, this estimator does have finite sample bias, which is an
important consideration for the coverage of the bootstrapping that resi uses.
The formula for the alternative conversion is:
Returns a scalar or vector argument of the the robust effect size index estimate.
# to obtain t values, first fit a lm mod = lm(charges ~ region * age + bmi + sex, data = RESI::insurance) # run lmtest::coeftest to get t values, using a robust variance-covariance formula ts = lmtest::coeftest(mod, vcov. = sandwich::vcovHC)[,'t value'] # get RESI estimates using unbiased estimator t2S(ts, n = nrow(RESI::insurance), rdf = mod$df.residual) # get RESI estimates using alternative estimator t2S(ts, n = nrow(RESI::insurance), rdf = mod$df.residual, unbiased = FALSE)# to obtain t values, first fit a lm mod = lm(charges ~ region * age + bmi + sex, data = RESI::insurance) # run lmtest::coeftest to get t values, using a robust variance-covariance formula ts = lmtest::coeftest(mod, vcov. = sandwich::vcovHC)[,'t value'] # get RESI estimates using unbiased estimator t2S(ts, n = nrow(RESI::insurance), rdf = mod$df.residual) # get RESI estimates using alternative estimator t2S(ts, n = nrow(RESI::insurance), rdf = mod$df.residual, unbiased = FALSE)
This function computes the robust effect size index from Vandekar, Tao, & Blume (2020). Vector arguments are accepted. If different length arguments are passed they are dealt with in the usual way of R.
z2S(z, n, unbiased = TRUE)z2S(z, n, unbiased = TRUE)
z |
The Z statistic for the parameter of interest. |
n |
Number of independent samples. |
unbiased |
Logical, whether to use unbiased or alternative estimator. See details. |
This function computes S, the RESI, from a Z statistic. The formula for the unbiased estimator (default) is derived by solving the expected value of the Z statistic for S. It is unbiased and consistent.
The formula for the unbiased conversion is:
The formula for the alternative estimator is derived by squaring the Z statistic
and using the chisq2S formula. This estimator may be appealing for its
intuitive relationship to the Chi-square statistic; the absolute value of RESI estimates
using this formula will be equal to a RESI estimate using a Chi-square statistic for
the same model. However, this estimator does have finite sample bias, which is an
important consideration for the coverage of the bootstrapping that resi uses.
The formula for the alternative conversion is:
Returns a scalar or vector argument of the the robust effect size index estimate.
# to obtain example z values, first fit a glm mod = glm(charges ~ region * age + bmi + sex, data = RESI::insurance) # run coeftest to get z values using a robust variance-covariance function zs = lmtest::coeftest(mod, vcov. = sandwich::vcovHC)[,'z value'] # get RESI estimates using unbiased estimator z2S(zs, n = nrow(RESI::insurance)) # get RESI estimates usng alternative estimator z2S(zs, n = nrow(RESI::insurance), unbiased = FALSE)# to obtain example z values, first fit a glm mod = glm(charges ~ region * age + bmi + sex, data = RESI::insurance) # run coeftest to get z values using a robust variance-covariance function zs = lmtest::coeftest(mod, vcov. = sandwich::vcovHC)[,'z value'] # get RESI estimates using unbiased estimator z2S(zs, n = nrow(RESI::insurance)) # get RESI estimates usng alternative estimator z2S(zs, n = nrow(RESI::insurance), unbiased = FALSE)