| Title: | Cluster-Robust (Sandwich) Variance Estimators with Small-Sample Corrections |
|---|---|
| Description: | Provides several cluster-robust variance estimators (i.e., sandwich estimators) for ordinary and weighted least squares linear regression models, including the bias-reduced linearization estimator introduced by Bell and McCaffrey (2002) <https://www150.statcan.gc.ca/n1/pub/12-001-x/2002002/article/9058-eng.pdf> and developed further by Pustejovsky and Tipton (2017) <DOI:10.1080/07350015.2016.1247004>. The package includes functions for estimating the variance- covariance matrix and for testing single- and multiple- contrast hypotheses based on Wald test statistics. Tests of single regression coefficients use Satterthwaite or saddle-point corrections. Tests of multiple- contrast hypotheses use an approximation to Hotelling's T-squared distribution. Methods are provided for a variety of fitted models, including lm() and mlm objects; glm(); geeglm() (from package 'geepack'); lm_robust(), lm_lin(), and iv_robust() (from package 'estimatr'); ivreg() (from package 'AER'); ivreg() (from package 'ivreg' when estimated by ordinary least squares); plm() (from package 'plm'); gls() and lme() (from 'nlme'); lmer() (from `lme4`); robu() (from 'robumeta'); rma.uni() and rma.mv() (from 'metafor'); and mmrm() (from 'mmrm'). |
| Authors: | James E. Pustejovsky [aut, cre] (ORCID: <https://orcid.org/0000-0003-0591-9465>), Hazim Izani [ctb], Samuel Pekofsky [ctb], Jingru Zhang [ctb] |
| Maintainer: | James E. Pustejovsky <[email protected]> |
| License: | GPL-3 |
| Version: | 0.7.0.9000 |
| Built: | 2026-06-03 06:55:33 UTC |
| Source: | https://github.com/jepusto/clubsandwich |
Data from a randomized trial of the Achievement Awards Demonstration program, reported in Angrist & Lavy (2009).
AchievementAwardsRCTAchievementAwardsRCT
A data frame with 16526 rows and 21 variables:
Fictitious school identification number
Factor identifying the school type (Arab religious, Jewish religious, Jewish secular)
Number of treatment pair. Note that 7 is a triple.
Indicator for whether school was in treatment group
Cohort year
Fictitious student identification number
Factor identifying student sex
Number of siblings
Indicator for immigrant status
Father's level of education
Mother's level of education
Indicator for Bagrut attainment
Number of Bagrut units attempted
Number of Bagrut units awarded
Indicator for satisfaction of math requirement
Indicator for satisfaction of English requirement
Indicator for satisfaction of Hebrew requirement
Lagged Bagrut score
Quartile within distribution of lagscore, calculated by cohort and sex
Lower or upper half within distribution of lagscore, calculated by cohort and sex
Angrist, J. D., & Lavy, V. (2009). The effects of high stakes high school achievement awards : Evidence from a randomized trial. American Economic Review, 99(4), 1384-1414. doi:10.1257/aer.99.4.1384
coef_test reports one- or two-sided t-tests for each coefficient
estimate in a fitted linear regression model, using a sandwich estimator for
the standard errors and (optionally) a small sample correction for the
p-value. Available small-sample corrections include Satterthwaite
approximation or a saddlepoint approximation. Coefficients can be tested
against non-zero null values by specifying null_constants.
coef_test( obj, vcov, test = "Satterthwaite", alternative = c("two-sided", "greater", "less"), coefs = "All", null_constants = 0, p_values = TRUE, ... )coef_test( obj, vcov, test = "Satterthwaite", alternative = c("two-sided", "greater", "less"), coefs = "All", null_constants = 0, p_values = TRUE, ... )
obj |
Fitted model for which to calculate t-tests. |
vcov |
Variance covariance matrix estimated using |
test |
Character vector specifying which small-sample corrections to
calculate. |
alternative |
Character string specifying the alternative hypothesis, with options "two-sided" (the default), "greater" or "less". |
coefs |
Character, integer, or logical vector specifying which
coefficients should be tested. The default value |
null_constants |
vector of null values for each coefficient to test.
Must have length equal to the number of coefficients specified in
|
p_values |
Logical indicating whether to report p-values. The default
value is |
... |
Further arguments passed to |
A data frame containing estimated regression coefficients, standard errors, specified values of null hypotheses, and test results. For the Satterthwaite approximation, degrees of freedom and a p-value are reported. For the saddlepoint approximation, the saddlepoint and a p-value are reported.
data("ChickWeight", package = "datasets") lm_fit <- lm(weight ~ Diet * Time, data = ChickWeight) diet_index <- grepl("Diet.:Time", names(coef(lm_fit))) coef_test(lm_fit, vcov = "CR2", cluster = ChickWeight$Chick, coefs = diet_index) V_CR2 <- vcovCR(lm_fit, cluster = ChickWeight$Chick, type = "CR2") coef_test(lm_fit, vcov = V_CR2, coefs = diet_index) # non-inferiority test whether time-by-diet interaction effects are 2 or greater coef_test(lm_fit, vcov = V_CR2, coefs = diet_index, null_constants = 2, alternative = "greater")data("ChickWeight", package = "datasets") lm_fit <- lm(weight ~ Diet * Time, data = ChickWeight) diet_index <- grepl("Diet.:Time", names(coef(lm_fit))) coef_test(lm_fit, vcov = "CR2", cluster = ChickWeight$Chick, coefs = diet_index) V_CR2 <- vcovCR(lm_fit, cluster = ChickWeight$Chick, type = "CR2") coef_test(lm_fit, vcov = V_CR2, coefs = diet_index) # non-inferiority test whether time-by-diet interaction effects are 2 or greater coef_test(lm_fit, vcov = V_CR2, coefs = diet_index, null_constants = 2, alternative = "greater")
conf_int reports confidence intervals for each coefficient estimate in
a fitted linear regression model, using a sandwich estimator for the standard
errors and a small sample correction for the critical values. The
small-sample correction is based on a Satterthwaite approximation.
conf_int( obj, vcov, level = 0.95, test = "Satterthwaite", coefs = "All", ..., p_values = FALSE )conf_int( obj, vcov, level = 0.95, test = "Satterthwaite", coefs = "All", ..., p_values = FALSE )
obj |
Fitted model for which to calculate confidence intervals. |
vcov |
Variance covariance matrix estimated using |
level |
Desired coverage level for confidence intervals. |
test |
Character vector specifying which small-sample corrections to
calculate. |
coefs |
Character, integer, or logical vector specifying which
coefficients should be tested. The default value |
... |
Further arguments passed to |
p_values |
Logical indicating whether to report p-values. The default
value is |
A data frame containing estimated regression coefficients, standard errors, confidence intervals, and (optionally) p-values.
data("ChickWeight", package = "datasets") lm_fit <- lm(weight ~ Diet * Time, data = ChickWeight) diet_index <- grepl("Diet.:Time", names(coef(lm_fit))) conf_int(lm_fit, vcov = "CR2", cluster = ChickWeight$Chick, coefs = diet_index) V_CR2 <- vcovCR(lm_fit, cluster = ChickWeight$Chick, type = "CR2") conf_int(lm_fit, vcov = V_CR2, level = .99, coefs = diet_index)data("ChickWeight", package = "datasets") lm_fit <- lm(weight ~ Diet * Time, data = ChickWeight) diet_index <- grepl("Diet.:Time", names(coef(lm_fit))) conf_int(lm_fit, vcov = "CR2", cluster = ChickWeight$Chick, coefs = diet_index) V_CR2 <- vcovCR(lm_fit, cluster = ChickWeight$Chick, type = "CR2") conf_int(lm_fit, vcov = V_CR2, level = .99, coefs = diet_index)
Helper functions to create common types of constraint matrices,
for use with Wald_test to conduct Wald-type tests of linear
contrasts from a fitted regression model.
constrain_zero(constraints, coefs, reg_ex = FALSE) constrain_equal(constraints, coefs, reg_ex = FALSE) constrain_pairwise(constraints, coefs, reg_ex = FALSE, with_zero = FALSE)constrain_zero(constraints, coefs, reg_ex = FALSE) constrain_equal(constraints, coefs, reg_ex = FALSE) constrain_pairwise(constraints, coefs, reg_ex = FALSE, with_zero = FALSE)
constraints |
Set of constraints to test. Can be logical (using
|
coefs |
Vector of coefficient estimates, used to determine the column
dimension of the constraint matrix. Can be omitted if the function is
called inside |
reg_ex |
Logical indicating whether |
with_zero |
Logical indicating whether coefficients should also be
compared to zero. Defaults to |
Constraints can be specified as character vectors, regular
expressions (with reg_ex = TRUE), integer vectors, or logical
vectors.
constrain_zero() Creates a matrix that constrains a specified set of
coefficients to all be equal to zero.
constrain_equal() Creates a matrix that constrains a specified set
of coefficients to all be equal.
constrain_pairwise() Creates a list of constraint matrices
consisting of all pairwise comparisons between a specified set of
coefficients. If with_zero = TRUE, then the list will also include a
set of constraint matrices comparing each coefficient to zero.
A matrix or list of matrices encoding the specified set of constraints.
if (requireNamespace("carData", quietly = TRUE)) withAutoprint({ data(Duncan, package = "carData") Duncan$cluster <- sample(LETTERS[1:8], size = nrow(Duncan), replace = TRUE) Duncan_fit <- lm(prestige ~ 0 + type + income + type:income + type:education, data=Duncan) # Note that type:income terms are interactions because main effect of income is included # but type:education terms are separate slopes for each unique level of type Duncan_coefs <- coef(Duncan_fit) # The following are all equivalent constrain_zero(constraints = c("typeprof:income","typewc:income"), coefs = Duncan_coefs) constrain_zero(constraints = ":income", coefs = Duncan_coefs, reg_ex = TRUE) constrain_zero(constraints = 5:6, coefs = Duncan_coefs) constrain_zero(constraints = c(FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE), coefs = Duncan_coefs) # The following are all equivalent constrain_equal(c("typebc:education","typeprof:education","typewc:education"), Duncan_coefs) constrain_equal(":education", Duncan_coefs, reg_ex = TRUE) constrain_equal(7:9, Duncan_coefs) constrain_equal(c(FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,TRUE,TRUE,TRUE), Duncan_coefs) # Test pairwise equality of the education slopes constrain_pairwise(":education", Duncan_coefs, reg_ex = TRUE) # Test pairwise equality of the income slopes, plus compare against zero constrain_pairwise(":income", Duncan_coefs, reg_ex = TRUE, with_zero = TRUE) })if (requireNamespace("carData", quietly = TRUE)) withAutoprint({ data(Duncan, package = "carData") Duncan$cluster <- sample(LETTERS[1:8], size = nrow(Duncan), replace = TRUE) Duncan_fit <- lm(prestige ~ 0 + type + income + type:income + type:education, data=Duncan) # Note that type:income terms are interactions because main effect of income is included # but type:education terms are separate slopes for each unique level of type Duncan_coefs <- coef(Duncan_fit) # The following are all equivalent constrain_zero(constraints = c("typeprof:income","typewc:income"), coefs = Duncan_coefs) constrain_zero(constraints = ":income", coefs = Duncan_coefs, reg_ex = TRUE) constrain_zero(constraints = 5:6, coefs = Duncan_coefs) constrain_zero(constraints = c(FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, FALSE, FALSE), coefs = Duncan_coefs) # The following are all equivalent constrain_equal(c("typebc:education","typeprof:education","typewc:education"), Duncan_coefs) constrain_equal(":education", Duncan_coefs, reg_ex = TRUE) constrain_equal(7:9, Duncan_coefs) constrain_equal(c(FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,TRUE,TRUE,TRUE), Duncan_coefs) # Test pairwise equality of the education slopes constrain_pairwise(":education", Duncan_coefs, reg_ex = TRUE) # Test pairwise equality of the income slopes, plus compare against zero constrain_pairwise(":income", Duncan_coefs, reg_ex = TRUE, with_zero = TRUE) })
A dataset containing estimated effect sizes, variances, and covariates from a meta-analysis of dropout prevention/intervention program effects, conducted by Wilson et al. (2011). Missing observations were imputed.
dropoutPreventiondropoutPrevention
A data frame with 385 rows and 18 variables:
log-odds ratio measuring the intervention effect
estimated sampling variance of the log-odds ratio
unique identifier for each study
unique identifier for each sample within a study
study design (randomized, matched, or non-randomized and unmatched)
outcome measure for the intervention effect is estimated (school dropout, school enrollment, graduation, graduation or GED receipt)
degree of evaluator independence (independent, indirect but influential, involved in planning but not delivery, involved in delivery)
level of implementation quality (clear problems, possible problems, no apparent problems)
Program delivery site (community, mixed, school classroom, school but outside of classroom)
Overall attrition (proportion)
pretest group-equivalence log-odds ratio
adjusted or unadjusted data used to calculate intervention effect
proportion of the sample that is male
proportion of the sample that is white
average age of the sample
program duration (in weeks)
program contact hours per week
indicator for the 32 studies with 3 or more effect sizes
Wilson, S. J., Lipsey, M. W., Tanner-Smith, E., Huang, C. H., & Steinka-Fry, K. T. (2011). Dropout prevention and intervention programs: Effects on school completion and dropout Among school-aged children and youth: A systematic review. _Campbell Systematic Reviews, 7_(1), 1-61. doi:10.4073/csr.2011.8
Wilson, S. J., Lipsey, M. W., Tanner-Smith, E., Huang, C. H., & Steinka-Fry, K. T. (2011). Dropout prevention and intervention programs: Effects on school completion and dropout Among school-aged children and youth: A systematic review. _Campbell Systematic Reviews, 7_(1), 1-61. doi:10.4073/csr.2011.8
Tipton, E., & Pustejovsky, J. E. (2015). Small-sample adjustments for tests of moderators and model fit using robust variance estimation in meta-regression. _Journal of Educational and Behavioral Statistics, 40_(6), 604-634. doi:10.3102/1076998615606099
findCluster.rma.mv returns a vector of ID variables for the highest level of clustering in a fitted rma.mv model.
findCluster.rma.mv(obj)findCluster.rma.mv(obj)
obj |
A fitted |
A a vector of ID variables for the highest level of clustering in obj.
if (requireNamespace("metafor", quietly = TRUE)) { library(metafor) data(dat.assink2016, package = "metadat") mfor_fit <- rma.mv(yi ~ year + deltype, V = vi, random = ~ 1 | study / esid, data = dat.assink2016) findCluster.rma.mv(mfor_fit) }if (requireNamespace("metafor", quietly = TRUE)) { library(metafor) data(dat.assink2016, package = "metadat") mfor_fit <- rma.mv(yi ~ year + deltype, V = vi, random = ~ 1 | study / esid, data = dat.assink2016) findCluster.rma.mv(mfor_fit) }
'r lifecycle::badge("superseded")'
This function is superseded by the vcalc provided by
the metafor package. Compared to impute_covariance_matrix,
vcalc provides many further features, includes a
data argument, and uses syntax that is consistent with other
functions in metafor.
impute_covariance_matrix calculates a block-diagonal covariance
matrix, given the marginal variances, the block structure, and an assumed
correlation structure. Can be used to create compound-symmetric structures,
AR(1) auto-correlated structures, or combinations thereof.
impute_covariance_matrix( vi, cluster, r, ti, ar1, smooth_vi = FALSE, subgroup = NULL, return_list = identical(as.factor(cluster), sort(as.factor(cluster))), check_PD = TRUE )impute_covariance_matrix( vi, cluster, r, ti, ar1, smooth_vi = FALSE, subgroup = NULL, return_list = identical(as.factor(cluster), sort(as.factor(cluster))), check_PD = TRUE )
vi |
Vector of variances |
cluster |
Vector indicating which effects belong to the same cluster. Effects with the same value of 'cluster' will be treated as correlated. |
r |
Vector or numeric value of assumed constant correlation(s) between effect size estimates from each study. |
ti |
Vector of time-points describing temporal spacing of effects, for use with auto-regressive correlation structures. |
ar1 |
Vector or numeric value of assumed AR(1) auto-correlation(s)
between effect size estimates from each study. If specified, then |
smooth_vi |
Logical indicating whether to smooth the marginal variances
by taking the average |
subgroup |
Vector of category labels describing sub-groups of effects. If non-null, effects that share the same category label and the same cluster will be treated as correlated, but effects with different category labels will be treated as uncorrelated, even if they come from the same cluster. |
return_list |
Optional logical indicating whether to return a list of matrices (with one entry per block) or the full variance-covariance matrix. |
check_PD |
Optional logical indicating whether to check whether each
covariance matrix is positive definite. If |
A block-diagonal variance-covariance matrix (possibly represented as
a list of matrices) with a specified structure. The structure depends on
whether the r argument, ar1 argument, or both arguments are
specified. Let denote the specified variance for effect
in cluster and be the covariance
between effects and in cluster
.
If only r is specified, each block of the variance-covariance
matrix will have a constant (compound symmetric) correlation, so that
where is the specified correlation
for cluster . If only a single value is given in r, then
it will be used for every cluster.
If only ar1 is specified, each block of the variance-covariance matrix will have an
AR(1) auto-correlation structure, so that
where is the specified auto-correlation
for cluster and and
are specified time-points corresponding to effects and
in cluster . If only a single value is given in
ar1, then it will be used for every cluster.
If both r and ar1 are specified, each block of the variance-covariance matrix will have combination of compound symmetric and an AR(1)
auto-correlation structures, so that
where is the specified constant correlation for cluster
, is the specified auto-correlation for
cluster and and are
specified time-points corresponding to effects and
in cluster . If only single values are given in
r or ar1, they will be used for every cluster.
If smooth_vi = TRUE, then all of the variances within cluster
will be set equal to the average variance of cluster
, i.e.,
for
and .
If cluster is appropriately sorted, then a list of matrices,
with one entry per cluster, will be returned by default. If cluster
is out of order, then the full variance-covariance matrix will be returned
by default. The output structure can be controlled with the optional
return_list argument.
if (requireNamespace("metafor", quietly = TRUE)) { library(metafor) # Constant correlation data(SATcoaching) V_list <- impute_covariance_matrix(vi = SATcoaching$V, cluster = SATcoaching$study, r = 0.66) MVFE <- rma.mv(d ~ 0 + test, V = V_list, data = SATcoaching) conf_int(MVFE, vcov = "CR2", cluster = SATcoaching$study) }if (requireNamespace("metafor", quietly = TRUE)) { library(metafor) # Constant correlation data(SATcoaching) V_list <- impute_covariance_matrix(vi = SATcoaching$V, cluster = SATcoaching$study, r = 0.66) MVFE <- rma.mv(d ~ 0 + test, V = V_list, data = SATcoaching) conf_int(MVFE, vcov = "CR2", cluster = SATcoaching$study) }
linear_contrast reports confidence intervals and (optionally) p-values
for linear contrasts of regression coefficients from a fitted model, using a
sandwich estimator for the standard errors and (optionally) a small sample
correction for the critical values. The default small-sample correction is
based on a Satterthwaite approximation.
linear_contrast( obj, vcov, contrasts, level = 0.95, test = "Satterthwaite", ..., p_values = FALSE, adjustment_method = "none" )linear_contrast( obj, vcov, contrasts, level = 0.95, test = "Satterthwaite", ..., p_values = FALSE, adjustment_method = "none" )
obj |
Fitted model for which to calculate confidence intervals. |
vcov |
Variance covariance matrix estimated using |
contrasts |
A contrast matrix, or a list of multiple contrast matrices to test. See details and examples. |
level |
Desired coverage level for confidence intervals. |
test |
Character vector specifying which small-sample corrections to
calculate. |
... |
Further arguments passed to |
p_values |
Logical indicating whether to report p-values. The default
value is |
adjustment_method |
Correction method, a |
Constraints can be specified directly as q X p matrices or
indirectly through constrain_pairwise,
constrain_equal, or constrain_zero.
A data frame containing estimated contrasts, standard errors, confidence intervals, and (optionally) p-values.
data("ChickWeight", package = "datasets") lm_fit <- lm(weight ~ 0 + Diet + Time:Diet, data = ChickWeight) # Pairwise comparisons of diet-by-time slopes linear_contrast(lm_fit, vcov = "CR2", cluster = ChickWeight$Chick, contrasts = constrain_pairwise("Diet.:Time", reg_ex = TRUE)) if (requireNamespace("carData", quietly = TRUE)) withAutoprint({ data(Duncan, package = "carData") Duncan$cluster <- sample(LETTERS[1:8], size = nrow(Duncan), replace = TRUE) Duncan_fit <- lm(prestige ~ 0 + type + income + type:income + type:education, data=Duncan) # Note that type:income terms are interactions because main effect of income is included # but type:education terms are separate slopes for each unique level of type # Pairwise comparisons of type-by-education slopes linear_contrast(Duncan_fit, vcov = "CR2", cluster = Duncan$cluster, contrasts = constrain_pairwise(":education", reg_ex = TRUE), test = "Satterthwaite") # Repeat the above, but with p-values displayed and multiple comparisons adjustment applied. linear_contrast(Duncan_fit, vcov = "CR2", cluster = Duncan$cluster, contrasts = constrain_pairwise(":education", reg_ex = TRUE), test = "Satterthwaite", p_values = TRUE, adjustment_method = "hochberg") # Pairwise comparisons of type-by-income interactions linear_contrast(Duncan_fit, vcov = "CR2", cluster = Duncan$cluster, contrasts = constrain_pairwise(":income", reg_ex = TRUE, with_zero = TRUE), test = "Satterthwaite") })data("ChickWeight", package = "datasets") lm_fit <- lm(weight ~ 0 + Diet + Time:Diet, data = ChickWeight) # Pairwise comparisons of diet-by-time slopes linear_contrast(lm_fit, vcov = "CR2", cluster = ChickWeight$Chick, contrasts = constrain_pairwise("Diet.:Time", reg_ex = TRUE)) if (requireNamespace("carData", quietly = TRUE)) withAutoprint({ data(Duncan, package = "carData") Duncan$cluster <- sample(LETTERS[1:8], size = nrow(Duncan), replace = TRUE) Duncan_fit <- lm(prestige ~ 0 + type + income + type:income + type:education, data=Duncan) # Note that type:income terms are interactions because main effect of income is included # but type:education terms are separate slopes for each unique level of type # Pairwise comparisons of type-by-education slopes linear_contrast(Duncan_fit, vcov = "CR2", cluster = Duncan$cluster, contrasts = constrain_pairwise(":education", reg_ex = TRUE), test = "Satterthwaite") # Repeat the above, but with p-values displayed and multiple comparisons adjustment applied. linear_contrast(Duncan_fit, vcov = "CR2", cluster = Duncan$cluster, contrasts = constrain_pairwise(":education", reg_ex = TRUE), test = "Satterthwaite", p_values = TRUE, adjustment_method = "hochberg") # Pairwise comparisons of type-by-income interactions linear_contrast(Duncan_fit, vcov = "CR2", cluster = Duncan$cluster, contrasts = constrain_pairwise(":income", reg_ex = TRUE, with_zero = TRUE), test = "Satterthwaite") })
A dataset containing state-level annual mortality rates for select causes of death, as well as data related to the minimum legal drinking age and alcohol consumption.
MortalityRatesMortalityRates
A data frame with 5508 rows and 12 variables:
Year of observation
identifier for state
Number of deaths
Population size
Proportion of 18-20 year-old population that is legally allowed to drink
Beer taxation rate
Beer consumption per capita
Wine consumption per capita
Spirits consumption per capita
Total alcohol consumption per capita
Mortality rate per 10,000
Cause of death
Mastering 'Metrics data archive
Angrist, J. D., and Pischke, J. S. (2014). _Mastering'metrics: the path from cause to effect_. Princeton University Press, 2014.
Carpenter, C., & Dobkin, C. (2011). The minimum legal drinking age and public health. _Journal of Economic Perspectives, 25_(2), 133-156. doi:10.1257/jep.25.2.133
'r lifecycle::badge("superseded")'
This function is superseded by the vcalc provided by
the metafor package. Compared to pattern_covariance_matrix,
vcalc provides many further features, includes a
data argument, and uses syntax that is consistent with other
functions in metafor.
pattern_covariance_matrix calculates a
block-diagonal covariance matrix, given the marginal variances, the block
structure, and an assumed correlation structure defined by a patterned
correlation matrix.
pattern_covariance_matrix( vi, cluster, pattern_level, r_pattern, r, smooth_vi = FALSE, subgroup = NULL, return_list = identical(as.factor(cluster), sort(as.factor(cluster))), check_PD = TRUE )pattern_covariance_matrix( vi, cluster, pattern_level, r_pattern, r, smooth_vi = FALSE, subgroup = NULL, return_list = identical(as.factor(cluster), sort(as.factor(cluster))), check_PD = TRUE )
vi |
Vector of variances |
cluster |
Vector indicating which effects belong to the same cluster. Effects with the same value of 'cluster' will be treated as correlated. |
pattern_level |
Vector of categories for each effect size, used to determine which entry of the pattern matrix will be used to impute a correlation. |
r_pattern |
Patterned correlation matrix with row and column names
corresponding to the levels of |
r |
Vector or numeric value of assumed constant correlation(s) between effect size estimates from each study. |
smooth_vi |
Logical indicating whether to smooth the marginal variances
by taking the average |
subgroup |
Vector of category labels describing sub-groups of effects. If non-null, effects that share the same category label and the same cluster will be treated as correlated, but effects with different category labels will be treated as uncorrelated, even if they come from the same cluster. |
return_list |
Optional logical indicating whether to return a list of matrices (with one entry per block) or the full variance-covariance matrix. |
check_PD |
Optional logical indicating whether to check whether each
covariance matrix is positive definite. If |
A block-diagonal variance-covariance matrix (possibly represented as
a list of matrices) with a specified correlation structure, defined by a
patterned correlation matrix. Let denote the specified
variance for effect in cluster and
be the covariance between effects and
in cluster . Let be the level
of the pattern variable for effect in cluster ,
taking a value in . A patterned correlation matrix
is defined as a set of correlations between pairs of effects taking each
possible combination of patterns. Formally, let be the
correlation between effects in categories and ,
respectively, where . Then the
covariance between effects and in cluster
is taken to be
Correlations between effect sizes within the same category are defined by the diagonal values of the pattern matrix, which may take values less than one.
Combinations of pattern levels that do not occur in the patterned correlation matrix will be set equal to r.
If smooth_vi = TRUE, then all of the variances within cluster
will be set equal to the average variance of cluster
, i.e.,
for
and .
If cluster is appropriately sorted, then a list of matrices,
with one entry per cluster, will be returned by default. If cluster
is out of order, then the full variance-covariance matrix will be returned
by default. The output structure can be controlled with the optional
return_list argument.
pkgs_available <- requireNamespace("metafor", quietly = TRUE) & requireNamespace("robumeta", quietly = TRUE) if (pkgs_available) { library(metafor) data(oswald2013, package = "robumeta") dat <- escalc(data = oswald2013, measure = "ZCOR", ri = R, ni = N) subset_ids <- unique(dat$Study)[1:20] dat <- subset(dat, Study %in% subset_ids) # make a patterned correlation matrix p_levels <- levels(dat$Crit.Cat) r_pattern <- 0.7^as.matrix(dist(1:length(p_levels))) diag(r_pattern) <- seq(0.75, 0.95, length.out = 6) rownames(r_pattern) <- colnames(r_pattern) <- p_levels # impute the covariance matrix using patterned correlations V_list <- pattern_covariance_matrix(vi = dat$vi, cluster = dat$Study, pattern_level = dat$Crit.Cat, r_pattern = r_pattern, smooth_vi = TRUE) # fit a model using imputed covariance matrix MVFE <- rma.mv(yi ~ 0 + Crit.Cat, V = V_list, random = ~ Crit.Cat | Study, data = dat) conf_int(MVFE, vcov = "CR2") }pkgs_available <- requireNamespace("metafor", quietly = TRUE) & requireNamespace("robumeta", quietly = TRUE) if (pkgs_available) { library(metafor) data(oswald2013, package = "robumeta") dat <- escalc(data = oswald2013, measure = "ZCOR", ri = R, ni = N) subset_ids <- unique(dat$Study)[1:20] dat <- subset(dat, Study %in% subset_ids) # make a patterned correlation matrix p_levels <- levels(dat$Crit.Cat) r_pattern <- 0.7^as.matrix(dist(1:length(p_levels))) diag(r_pattern) <- seq(0.75, 0.95, length.out = 6) rownames(r_pattern) <- colnames(r_pattern) <- p_levels # impute the covariance matrix using patterned correlations V_list <- pattern_covariance_matrix(vi = dat$vi, cluster = dat$Study, pattern_level = dat$Crit.Cat, r_pattern = r_pattern, smooth_vi = TRUE) # fit a model using imputed covariance matrix MVFE <- rma.mv(yi ~ 0 + Crit.Cat, V = V_list, random = ~ Crit.Cat | Study, data = dat) conf_int(MVFE, vcov = "CR2") }
Effect sizes from studies on the effects of SAT coaching, reported in Kalaian and Raudenbush (1996)
SATcoachingSATcoaching
A data frame with 67 rows and 11 variables:
Study identifier
Year of publication
Character string indicating whether effect size corresponds to outcome on verbal (SATV) or math (SATM) test
Effect size estimate (Standardized mean difference)
Variance of effect size estimate
Sample size in treatment condition
Sample size in control condition
Character string indicating whether study design used a matched, non-equivalent, or randomized control group
Hours of coaching
Indicator variable for Educational Testing Service
Indicator variable for homework
Kalaian, H. A. & Raudenbush, S. W. (1996). A multivariate mixed linear model for meta-analysis. Psychological Methods, 1(3), 227-235. doi:10.1037/1082-989X.1.3.227
This is a generic function, with specific methods defined for
lm, plm, glm,
gls, lme,
robu, rma.uni,
rma.mv, and lm_robust
objects.
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates.
vcovCR(obj, cluster, type, target, inverse_var, form, ...) ## Default S3 method: vcovCR( obj, cluster, type, target = NULL, inverse_var = FALSE, form = "sandwich", ... )vcovCR(obj, cluster, type, target, inverse_var, form, ...) ## Default S3 method: vcovCR( obj, cluster, type, target = NULL, inverse_var = FALSE, form = "sandwich", ... )
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Expression or vector indicating which observations belong to the same cluster. For some classes, the cluster will be detected automatically if not specified. |
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Optional logical indicating whether the weights used in
fitting the model are inverse-variance. If not specified, |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates.
Several different small sample corrections are available, which run
parallel with the "HC" corrections for heteroskedasticity-consistent
variance estimators, as implemented in vcovHC. The
"CR2" adjustment is recommended (Pustejovsky & Tipton, 2017; Imbens &
Kolesar, 2016). See Pustejovsky and Tipton (2017) and Cameron and Miller
(2015) for further technical details. Available options include:
is the original form of the sandwich estimator (Liang & Zeger, 1986), which does not make any small-sample correction.
multiplies CR0 by m / (m - 1), where m is the
number of clusters.
multiplies CR0 by m / (m - p), where m is the
number of clusters and p is the number of covariates.
multiplies CR0 by (m (N-1)) / [(m -
1)(N - p)], where m is the number of clusters, N is the
total number of observations, and p is the number of covariates.
Some Stata commands use this correction by default.
is the "bias-reduced linearization" adjustment proposed by Bell and McCaffrey (2002) and further developed in Pustejovsky and Tipton (2017). The adjustment is chosen so that the variance-covariance estimator is exactly unbiased under a user-specified working model.
approximates the leave-one-cluster-out jackknife variance estimator (Bell & McCaffrey, 2002).
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates. The matrix has several attributes:
indicates which small-sample adjustment was used
contains the factor vector that defines independent clusters
contains the bread matrix
constant used in scaling the sandwich estimator
contains a list of estimating matrices used to calculate the sandwich estimator
contains a list of adjustment matrices used to calculate the sandwich estimator
contains the working variance-covariance model used to calculate the adjustment matrices. This is needed for calculating small-sample corrections for Wald tests.
Bell, R. M., & McCaffrey, D. F. (2002). Bias reduction in standard errors for linear regression with multi-stage samples. Survey Methodology, 28(2), 169-181.
Cameron, A. C., & Miller, D. L. (2015). A Practitioner's Guide to Cluster-Robust Inference. Journal of Human Resources, 50(2), 317-372. doi:10.3368/jhr.50.2.317
Imbens, G. W., & Kolesar, M. (2016). Robust standard errors in small samples: Some practical advice. Review of Economics and Statistics, 98(4), 701-712. doi:10.1162/rest_a_00552
Liang, K.-Y., & Zeger, S. L. (1986). Longitudinal data analysis using generalized linear models. Biometrika, 73(1), 13-22. doi:10.1093/biomet/73.1.13
Pustejovsky, J. E. & Tipton, E. (2018). Small sample methods for cluster-robust variance estimation and hypothesis testing in fixed effects models. Journal of Business and Economic Statistics, 36(4), 672-683. doi:10.1080/07350015.2016.1247004
vcovCR.lm, vcovCR.plm,
vcovCR.glm, vcovCR.gls,
vcovCR.lme, vcovCR.lmerMod, vcovCR.robu,
vcovCR.rma.uni, vcovCR.rma.mv,
vcovCR.lm_robust
# simulate design with cluster-dependence m <- 8 cluster <- factor(rep(LETTERS[1:m], 3 + rpois(m, 5))) n <- length(cluster) X <- matrix(rnorm(3 * n), n, 3) nu <- rnorm(m)[cluster] e <- rnorm(n) y <- X %*% c(.4, .3, -.3) + nu + e dat <- data.frame(y, X, cluster, row = 1:n) # fit linear model lm_fit <- lm(y ~ X1 + X2 + X3, data = dat) vcov(lm_fit) # cluster-robust variance estimator with CR2 small-sample correction vcovCR(lm_fit, cluster = dat$cluster, type = "CR2") # compare small-sample adjustments CR_types <- paste0("CR",c("0","1","1S","2","3")) sapply(CR_types, function(type) sqrt(diag(vcovCR(lm_fit, cluster = dat$cluster, type = type))))# simulate design with cluster-dependence m <- 8 cluster <- factor(rep(LETTERS[1:m], 3 + rpois(m, 5))) n <- length(cluster) X <- matrix(rnorm(3 * n), n, 3) nu <- rnorm(m)[cluster] e <- rnorm(n) y <- X %*% c(.4, .3, -.3) + nu + e dat <- data.frame(y, X, cluster, row = 1:n) # fit linear model lm_fit <- lm(y ~ X1 + X2 + X3, data = dat) vcov(lm_fit) # cluster-robust variance estimator with CR2 small-sample correction vcovCR(lm_fit, cluster = dat$cluster, type = "CR2") # compare small-sample adjustments CR_types <- paste0("CR",c("0","1","1S","2","3")) sapply(CR_types, function(type) sqrt(diag(vcovCR(lm_fit, cluster = dat$cluster, type = type))))
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from an geeglm object.
## S3 method for class 'geeglm' vcovCR( obj, cluster, type, target = NULL, inverse_var = NULL, form = "sandwich", ... )## S3 method for class 'geeglm' vcovCR( obj, cluster, type, target = NULL, inverse_var = NULL, form = "sandwich", ... )
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Expression or vector indicating which observations belong to
the same cluster. Required for |
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Optional logical indicating whether the weights used in
fitting the model are inverse-variance. If not specified, |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
if (requireNamespace("geepack", quietly = TRUE)) { library(geepack) data(dietox, package = "geepack") dietox$Cu <- as.factor(dietox$Cu) mf <- formula(Weight ~ Cu * (Time + I(Time^2) + I(Time^3))) gee1 <- geeglm(mf, data=dietox, id=Pig, family=poisson("identity"), corstr="ar1") V_CR <- vcovCR(gee1, cluster = dietox$Pig, type = "CR2") coef_test(gee1, vcov = V_CR, test = "Satterthwaite") }if (requireNamespace("geepack", quietly = TRUE)) { library(geepack) data(dietox, package = "geepack") dietox$Cu <- as.factor(dietox$Cu) mf <- formula(Weight ~ Cu * (Time + I(Time^2) + I(Time^3))) gee1 <- geeglm(mf, data=dietox, id=Pig, family=poisson("identity"), corstr="ar1") V_CR <- vcovCR(gee1, cluster = dietox$Pig, type = "CR2") coef_test(gee1, vcov = V_CR, test = "Satterthwaite") }
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from an glm object.
## S3 method for class 'glm' vcovCR( obj, cluster, type, target = NULL, inverse_var = NULL, form = "sandwich", ... )## S3 method for class 'glm' vcovCR( obj, cluster, type, target = NULL, inverse_var = NULL, form = "sandwich", ... )
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Expression or vector indicating which observations belong to
the same cluster. Required for |
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Optional logical indicating whether the weights used in
fitting the model are inverse-variance. If not specified, |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
if (requireNamespace("geepack", quietly = TRUE)) { data(dietox, package = "geepack") dietox$Cu <- as.factor(dietox$Cu) weight_fit <- glm(Weight ~ Cu * poly(Time, 3), data=dietox, family = "quasipoisson") V_CR <- vcovCR(weight_fit, cluster = dietox$Pig, type = "CR2") coef_test(weight_fit, vcov = V_CR, test = "Satterthwaite") }if (requireNamespace("geepack", quietly = TRUE)) { data(dietox, package = "geepack") dietox$Cu <- as.factor(dietox$Cu) weight_fit <- glm(Weight ~ Cu * poly(Time, 3), data=dietox, family = "quasipoisson") V_CR <- vcovCR(weight_fit, cluster = dietox$Pig, type = "CR2") coef_test(weight_fit, vcov = V_CR, test = "Satterthwaite") }
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from a gls object.
## S3 method for class 'gls' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)## S3 method for class 'gls' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Optional expression or vector indicating which observations
belong to the same cluster. If not specified, will be set to
|
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Optional logical indicating whether the weights used in
fitting the model are inverse-variance. If not specified, |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
if (requireNamespace("nlme", quietly = TRUE)) { library(nlme) data(Ovary, package = "nlme") Ovary$time_int <- 1:nrow(Ovary) lm_AR1 <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), data = Ovary, correlation = corAR1(form = ~ time_int | Mare)) vcovCR(lm_AR1, type = "CR2") }if (requireNamespace("nlme", quietly = TRUE)) { library(nlme) data(Ovary, package = "nlme") Ovary$time_int <- 1:nrow(Ovary) lm_AR1 <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), data = Ovary, correlation = corAR1(form = ~ time_int | Mare)) vcovCR(lm_AR1, type = "CR2") }
estimatr::iv_robust object.vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from an
iv_robust object.
## S3 method for class 'iv_robust' vcovCR( obj, cluster, type, target = NULL, inverse_var = FALSE, form = "sandwich", ... )## S3 method for class 'iv_robust' vcovCR( obj, cluster, type, target = NULL, inverse_var = FALSE, form = "sandwich", ... )
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Vector indicating which observations belong to
the same cluster. Required for |
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Not used for |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
if (requireNamespace("estimatr", quietly = TRUE)) withAutoprint({ library(estimatr) set.seed(20261201) N <- 200 cluster <- factor(rep(1:20, each = 10)) z1 <- rnorm(N) z2 <- rnorm(N) x_exog <- rnorm(N) x_endog <- 0.5 * z1 + 0.3 * z2 + rnorm(N) y <- 1 + 2 * x_endog + 0.5 * x_exog + rnorm(N) dat <- data.frame(y, x_endog, x_exog, z1, z2, cluster) # Basic 2SLS with clusters and CR2 iv_fit <- iv_robust( y ~ x_endog + x_exog | x_exog + z1 + z2, data = dat, clusters = cluster, se_type = "CR2" ) vcovCR(iv_fit) conf_int(iv_fit, vcov = "CR2") # 2SLS with fixed effects iv_fe <- iv_robust( y ~ x_endog + x_exog | x_exog + z1 + z2, fixed_effects = ~ cluster, data = dat, clusters = cluster, se_type = "CR2" ) vcovCR(iv_fe) conf_int(iv_fe, vcov = "CR2") })if (requireNamespace("estimatr", quietly = TRUE)) withAutoprint({ library(estimatr) set.seed(20261201) N <- 200 cluster <- factor(rep(1:20, each = 10)) z1 <- rnorm(N) z2 <- rnorm(N) x_exog <- rnorm(N) x_endog <- 0.5 * z1 + 0.3 * z2 + rnorm(N) y <- 1 + 2 * x_endog + 0.5 * x_exog + rnorm(N) dat <- data.frame(y, x_endog, x_exog, z1, z2, cluster) # Basic 2SLS with clusters and CR2 iv_fit <- iv_robust( y ~ x_endog + x_exog | x_exog + z1 + z2, data = dat, clusters = cluster, se_type = "CR2" ) vcovCR(iv_fit) conf_int(iv_fit, vcov = "CR2") # 2SLS with fixed effects iv_fe <- iv_robust( y ~ x_endog + x_exog | x_exog + z1 + z2, fixed_effects = ~ cluster, data = dat, clusters = cluster, se_type = "CR2" ) vcovCR(iv_fe) conf_int(iv_fe, vcov = "CR2") })
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from an ivreg object fitted
from the AER package or the ivreg package.
## S3 method for class 'ivreg' vcovCR( obj, cluster, type, target = NULL, inverse_var = FALSE, form = "sandwich", ... )## S3 method for class 'ivreg' vcovCR( obj, cluster, type, target = NULL, inverse_var = FALSE, form = "sandwich", ... )
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Expression or vector indicating which observations belong to
the same cluster. Required for |
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Not used for |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
For any "ivreg" objects fitted via the ivreg
function from the ivreg package, only traditional 2SLS
regression method (method = "OLS") is supported.
clubSandwich currently cannot support robust-regression methods such as
M-estimation (method = "M") or MM-estimation (method = "MM").
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
if (requireNamespace("AER", quietly = TRUE)) withAutoprint({ library(AER) data("CigarettesSW") Cigs <- within(CigarettesSW, { rprice <- price/cpi rincome <- income/population/cpi tdiff <- (taxs - tax)/cpi }) iv_fit_AER <- AER::ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), data = Cigs) vcovCR(iv_fit_AER, cluster = Cigs$state, type = "CR2") coef_test(iv_fit_AER, vcov = "CR2", cluster = Cigs$state) }) pkgs_available <- requireNamespace("AER", quietly = TRUE) & requireNamespace("ivreg", quietly = TRUE) if (pkgs_available) withAutoprint ({ data("CigarettesSW") Cigs <- within(CigarettesSW, { rprice <- price/cpi rincome <- income/population/cpi tdiff <- (taxs - tax)/cpi }) iv_fit_ivreg <- ivreg::ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), data = Cigs) vcovCR(iv_fit_ivreg, cluster = Cigs$state, type = "CR2") coef_test(iv_fit_ivreg, vcov = "CR2", cluster = Cigs$state) })if (requireNamespace("AER", quietly = TRUE)) withAutoprint({ library(AER) data("CigarettesSW") Cigs <- within(CigarettesSW, { rprice <- price/cpi rincome <- income/population/cpi tdiff <- (taxs - tax)/cpi }) iv_fit_AER <- AER::ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), data = Cigs) vcovCR(iv_fit_AER, cluster = Cigs$state, type = "CR2") coef_test(iv_fit_AER, vcov = "CR2", cluster = Cigs$state) }) pkgs_available <- requireNamespace("AER", quietly = TRUE) & requireNamespace("ivreg", quietly = TRUE) if (pkgs_available) withAutoprint ({ data("CigarettesSW") Cigs <- within(CigarettesSW, { rprice <- price/cpi rincome <- income/population/cpi tdiff <- (taxs - tax)/cpi }) iv_fit_ivreg <- ivreg::ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), data = Cigs) vcovCR(iv_fit_ivreg, cluster = Cigs$state, type = "CR2") coef_test(iv_fit_ivreg, vcov = "CR2", cluster = Cigs$state) })
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from an lm object.
## S3 method for class 'lm' vcovCR( obj, cluster, type, target = NULL, inverse_var = NULL, form = "sandwich", ... )## S3 method for class 'lm' vcovCR( obj, cluster, type, target = NULL, inverse_var = NULL, form = "sandwich", ... )
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Expression or vector indicating which observations belong to
the same cluster. Required for |
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Optional logical indicating whether the weights used in
fitting the model are inverse-variance. If not specified, |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
data("ChickWeight", package = "datasets") lm_fit <- lm(weight ~ Time + Diet:Time, data = ChickWeight) vcovCR(lm_fit, cluster = ChickWeight$Chick, type = "CR2") if (requireNamespace("plm", quietly = TRUE)) withAutoprint({ data("Produc", package = "plm") lm_individual <- lm(log(gsp) ~ 0 + state + log(pcap) + log(pc) + log(emp) + unemp, data = Produc) individual_index <- !grepl("state", names(coef(lm_individual))) vcovCR(lm_individual, cluster = Produc$state, type = "CR2")[individual_index,individual_index] # compare to plm() plm_FE <- plm::plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = Produc, index = c("state","year"), effect = "individual", model = "within") vcovCR(plm_FE, type="CR2") })data("ChickWeight", package = "datasets") lm_fit <- lm(weight ~ Time + Diet:Time, data = ChickWeight) vcovCR(lm_fit, cluster = ChickWeight$Chick, type = "CR2") if (requireNamespace("plm", quietly = TRUE)) withAutoprint({ data("Produc", package = "plm") lm_individual <- lm(log(gsp) ~ 0 + state + log(pcap) + log(pc) + log(emp) + unemp, data = Produc) individual_index <- !grepl("state", names(coef(lm_individual))) vcovCR(lm_individual, cluster = Produc$state, type = "CR2")[individual_index,individual_index] # compare to plm() plm_FE <- plm::plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = Produc, index = c("state","year"), effect = "individual", model = "within") vcovCR(plm_FE, type="CR2") })
estimatr::lm_robust
object.vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from an
lm_robust object.
## S3 method for class 'lm_robust' vcovCR( obj, cluster, type, target = NULL, inverse_var = NULL, form = "sandwich", ... )## S3 method for class 'lm_robust' vcovCR( obj, cluster, type, target = NULL, inverse_var = NULL, form = "sandwich", ... )
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Expression or vector indicating which observations belong to
the same cluster. If not specified, will be detected from the
|
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Optional logical indicating whether the weights used in
fitting the model are inverse-variance. If not specified, |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
data("ChickWeight", package = "datasets") ChickWeight$Chick <- factor(ChickWeight$Chick, ordered = FALSE) if (requireNamespace("estimatr", quietly = TRUE)) withAutoprint({ library(estimatr) lm_fit <- lm_robust( weight ~ Time + Diet:Time, data = ChickWeight ) vcovCR(lm_fit, cluster = ChickWeight$Chick, type = "CR2") lm_fit_clust <- lm_robust( weight ~ Time + Diet:Time, data = ChickWeight, clusters = Chick ) conf_int(lm_fit_clust, vcov = "CR2") # similar model via lm_lin() lin_fit_clust <- lm_lin( weight ~ Diet, covariates = ~ Time, data = ChickWeight, clusters = Chick ) conf_int(lin_fit_clust, vcov = "CR2") lm_fit_fe <- lm_robust( weight ~ Time:Diet, data = ChickWeight, clusters = Chick, fixed_effects = ~ Chick ) vcovCR(lm_fit_fe) # two-way fixed effects model data("MortalityRates") MortalityRates <- subset(MortalityRates, cause == "Motor Vehicle") MortalityRates$state <- factor(MortalityRates$state) MortalityRates$year <- factor(MortalityRates$year) MLDA_fit <- lm_robust( mrate ~ legal + beertaxa + beerpercap + winepercap + spiritpercap, fixed_effects = ~ year + state, data = MortalityRates, cluster = state ) conf_int(MLDA_fit, vcov = "CR2") if (requireNamespace("plm", quietly = TRUE)) withAutoprint({ data("Produc", package = "plm") lm_individual <- lm_robust( log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = Produc, fixed_effects = ~ state, cluster = state ) vcovCR(lm_individual, type = "CR2") }) })data("ChickWeight", package = "datasets") ChickWeight$Chick <- factor(ChickWeight$Chick, ordered = FALSE) if (requireNamespace("estimatr", quietly = TRUE)) withAutoprint({ library(estimatr) lm_fit <- lm_robust( weight ~ Time + Diet:Time, data = ChickWeight ) vcovCR(lm_fit, cluster = ChickWeight$Chick, type = "CR2") lm_fit_clust <- lm_robust( weight ~ Time + Diet:Time, data = ChickWeight, clusters = Chick ) conf_int(lm_fit_clust, vcov = "CR2") # similar model via lm_lin() lin_fit_clust <- lm_lin( weight ~ Diet, covariates = ~ Time, data = ChickWeight, clusters = Chick ) conf_int(lin_fit_clust, vcov = "CR2") lm_fit_fe <- lm_robust( weight ~ Time:Diet, data = ChickWeight, clusters = Chick, fixed_effects = ~ Chick ) vcovCR(lm_fit_fe) # two-way fixed effects model data("MortalityRates") MortalityRates <- subset(MortalityRates, cause == "Motor Vehicle") MortalityRates$state <- factor(MortalityRates$state) MortalityRates$year <- factor(MortalityRates$year) MLDA_fit <- lm_robust( mrate ~ legal + beertaxa + beerpercap + winepercap + spiritpercap, fixed_effects = ~ year + state, data = MortalityRates, cluster = state ) conf_int(MLDA_fit, vcov = "CR2") if (requireNamespace("plm", quietly = TRUE)) withAutoprint({ data("Produc", package = "plm") lm_individual <- lm_robust( log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = Produc, fixed_effects = ~ state, cluster = state ) vcovCR(lm_individual, type = "CR2") }) })
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from a lme object.
## S3 method for class 'lme' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)## S3 method for class 'lme' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Optional expression or vector indicating which observations
belong to the same cluster. If not specified, will be set to
|
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Optional logical indicating whether the weights used in
fitting the model are inverse-variance. If not specified, |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
if (requireNamespace("nlme", quietly = TRUE)) { library(nlme) rat_weight <- lme(weight ~ Time * Diet, data=BodyWeight, ~ Time | Rat) vcovCR(rat_weight, type = "CR2") } pkgs_available <- requireNamespace("nlme", quietly = TRUE) & requireNamespace("mlmRev", quietly = TRUE) if (pkgs_available) { data(egsingle, package = "mlmRev") subset_ids <- levels(egsingle$schoolid)[1:10] egsingle_subset <- subset(egsingle, schoolid %in% subset_ids) math_model <- lme(math ~ year * size + female + black + hispanic, random = list(~ year | schoolid, ~ 1 | childid), data = egsingle_subset) vcovCR(math_model, type = "CR2") }if (requireNamespace("nlme", quietly = TRUE)) { library(nlme) rat_weight <- lme(weight ~ Time * Diet, data=BodyWeight, ~ Time | Rat) vcovCR(rat_weight, type = "CR2") } pkgs_available <- requireNamespace("nlme", quietly = TRUE) & requireNamespace("mlmRev", quietly = TRUE) if (pkgs_available) { data(egsingle, package = "mlmRev") subset_ids <- levels(egsingle$schoolid)[1:10] egsingle_subset <- subset(egsingle, schoolid %in% subset_ids) math_model <- lme(math ~ year * size + female + black + hispanic, random = list(~ year | schoolid, ~ 1 | childid), data = egsingle_subset) vcovCR(math_model, type = "CR2") }
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from merMod object.
## S3 method for class 'lmerMod' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)## S3 method for class 'lmerMod' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Optional expression or vector indicating which observations
belong to the same cluster. If not specified, will be set to
|
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Optional logical indicating whether the weights used in
fitting the model are inverse-variance. If not specified, |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
if (requireNamespace("lme4", quietly = TRUE)) { library(lme4) sleep_fit <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy) vcovCR(sleep_fit, type = "CR2") } pkgs_available <- requireNamespace("lme4", quietly = TRUE) & requireNamespace("mlmRev", quietly = TRUE) if (pkgs_available) { data(egsingle, package = "mlmRev") subset_ids <- levels(egsingle$schoolid)[1:10] math_model <- lmer(math ~ year * size + female + black + hispanic + (1 | schoolid) + (1 | childid), data = egsingle, subset = schoolid %in% subset_ids) vcovCR(math_model, type = "CR2") }if (requireNamespace("lme4", quietly = TRUE)) { library(lme4) sleep_fit <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy) vcovCR(sleep_fit, type = "CR2") } pkgs_available <- requireNamespace("lme4", quietly = TRUE) & requireNamespace("mlmRev", quietly = TRUE) if (pkgs_available) { data(egsingle, package = "mlmRev") subset_ids <- levels(egsingle$schoolid)[1:10] math_model <- lmer(math ~ year * size + female + black + hispanic + (1 | schoolid) + (1 | childid), data = egsingle, subset = schoolid %in% subset_ids) vcovCR(math_model, type = "CR2") }
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from an mlm object.
## S3 method for class 'mlm' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)## S3 method for class 'mlm' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Optional expression or vector indicating which observations belong to the same cluster. If not specified, each row of the data will be treated as a separate cluster. |
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Optional logical indicating whether the weights used in
fitting the model are inverse-variance. If not specified, |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
iris_fit <- lm(cbind(Sepal.Length, Sepal.Width) ~ Species + Petal.Length + Petal.Width, data = iris) Vcluster <- vcovCR(iris_fit, type = "CR2") Vclusteriris_fit <- lm(cbind(Sepal.Length, Sepal.Width) ~ Species + Petal.Length + Petal.Width, data = iris) Vcluster <- vcovCR(iris_fit, type = "CR2") Vcluster
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from an
mmrm object.
## S3 method for class 'mmrm' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)## S3 method for class 'mmrm' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Optional expression or vector indicating which observations
belong to the same cluster. If not specified, will be set to the subject
variable from the |
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Optional logical indicating whether the weights used in
fitting the model are inverse-variance. If not specified, |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
if (requireNamespace("mmrm", quietly = TRUE)) withAutoprint({ library(mmrm) data(fev_data, package = "mmrm") # Fit an mmrm model with unstructured covariance mmrm_fit <- mmrm( FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID), data = fev_data ) # CR2 cluster-robust variance estimator (cluster auto-detected) vcovCR(mmrm_fit, type = "CR2") # Coefficient tests with Satterthwaite degrees of freedom coef_test(mmrm_fit, vcov = "CR2", test = "Satterthwaite") # Fit a weighted mmrm model fev_data$wt <- 1 + rpois(nrow(fev_data), lambda = 3) mmrm_wt <- mmrm( FEV1 ~ RACE + SEX + ARMCD + us(AVISIT | USUBJID), data = fev_data, weights = fev_data$wt ) # CR2 works with weighted models vcovCR(mmrm_wt, type = "CR2") })if (requireNamespace("mmrm", quietly = TRUE)) withAutoprint({ library(mmrm) data(fev_data, package = "mmrm") # Fit an mmrm model with unstructured covariance mmrm_fit <- mmrm( FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID), data = fev_data ) # CR2 cluster-robust variance estimator (cluster auto-detected) vcovCR(mmrm_fit, type = "CR2") # Coefficient tests with Satterthwaite degrees of freedom coef_test(mmrm_fit, vcov = "CR2", test = "Satterthwaite") # Fit a weighted mmrm model fev_data$wt <- 1 + rpois(nrow(fev_data), lambda = 3) mmrm_wt <- mmrm( FEV1 ~ RACE + SEX + ARMCD + us(AVISIT | USUBJID), data = fev_data, weights = fev_data$wt ) # CR2 works with weighted models vcovCR(mmrm_wt, type = "CR2") })
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from a plm
object.
## S3 method for class 'plm' vcovCR( obj, cluster, type, target, inverse_var, form = "sandwich", ignore_FE = FALSE, ... )## S3 method for class 'plm' vcovCR( obj, cluster, type, target, inverse_var, form = "sandwich", ignore_FE = FALSE, ... )
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Optional character string, expression, or vector indicating
which observations belong to the same cluster. For fixed-effect models that
include individual effects or time effects (but not both), the cluster will
be taken equal to the included fixed effects if not otherwise specified.
Clustering on individuals can also be obtained by specifying the name of
the individual index (e.g., |
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Optional logical indicating whether the weights used in
fitting the model are inverse-variance. If not specified, |
form |
Controls the form of the returned matrix. The default
|
ignore_FE |
Optional logical controlling whether fixed effects are ignored when calculating small-sample adjustments in models where fixed effects are estimated through absorption. |
... |
Additional arguments available for some classes of objects. |
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
if (requireNamespace("plm", quietly = TRUE)) withAutoprint({ library(plm) # fixed effects data("Produc", package = "plm") plm_FE <- plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = Produc, index = c("state","year","region"), effect = "individual", model = "within") vcovCR(plm_FE, type="CR2") vcovCR(plm_FE, type = "CR2", cluster = Produc$region) # clustering on region # random effects plm_RE <- update(plm_FE, model = "random") vcovCR(plm_RE, type = "CR2") vcovCR(plm_RE, type = "CR2", cluster = Produc$region) # clustering on region # nested random effects plm_nested <- update(plm_FE, effect = "nested", model = "random") vcovCR(plm_nested, type = "CR2") # clustering on region }) pkgs_available <- requireNamespace("plm", quietly = TRUE) & requireNamespace("AER", quietly = TRUE) if (pkgs_available) withAutoprint({ # first differencing data(Fatalities, package = "AER") Fatalities <- within(Fatalities, { frate <- 10000 * fatal / pop drinkagec <- cut(drinkage, breaks = 18:22, include.lowest = TRUE, right = FALSE) drinkagec <- relevel(drinkagec, ref = 4) }) plm_FD <- plm(frate ~ beertax + drinkagec + miles + unemp + log(income), data = Fatalities, index = c("state", "year"), model = "fd") vcovHC(plm_FD, method="arellano", type = "sss", cluster = "group") vcovCR(plm_FD, type = "CR1S") vcovCR(plm_FD, type = "CR2") })if (requireNamespace("plm", quietly = TRUE)) withAutoprint({ library(plm) # fixed effects data("Produc", package = "plm") plm_FE <- plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = Produc, index = c("state","year","region"), effect = "individual", model = "within") vcovCR(plm_FE, type="CR2") vcovCR(plm_FE, type = "CR2", cluster = Produc$region) # clustering on region # random effects plm_RE <- update(plm_FE, model = "random") vcovCR(plm_RE, type = "CR2") vcovCR(plm_RE, type = "CR2", cluster = Produc$region) # clustering on region # nested random effects plm_nested <- update(plm_FE, effect = "nested", model = "random") vcovCR(plm_nested, type = "CR2") # clustering on region }) pkgs_available <- requireNamespace("plm", quietly = TRUE) & requireNamespace("AER", quietly = TRUE) if (pkgs_available) withAutoprint({ # first differencing data(Fatalities, package = "AER") Fatalities <- within(Fatalities, { frate <- 10000 * fatal / pop drinkagec <- cut(drinkage, breaks = 18:22, include.lowest = TRUE, right = FALSE) drinkagec <- relevel(drinkagec, ref = 4) }) plm_FD <- plm(frate ~ beertax + drinkagec + miles + unemp + log(income), data = Fatalities, index = c("state", "year"), model = "fd") vcovHC(plm_FD, method="arellano", type = "sss", cluster = "group") vcovCR(plm_FD, type = "CR1S") vcovCR(plm_FD, type = "CR2") })
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from a
rma.mv object.
## S3 method for class 'rma.mv' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)## S3 method for class 'rma.mv' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Optional expression or vector indicating which observations belong to the same cluster. If not specified, will be set to the factor in the random-effects structure with the fewest distinct levels. Caveat emptor: the function does not check that the random effects are nested. |
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Optional logical indicating whether the weights used in
fitting the model are inverse-variance. If not specified, |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
pkgs_available <- requireNamespace("metafor", quietly = TRUE) & requireNamespace("metadat", quietly = TRUE) if (pkgs_available) withAutoprint({ library(metafor) data(dat.assink2016, package = "metadat") mfor_fit <- rma.mv(yi ~ year + deltype, V = vi, random = ~ 1 | study / esid, data = dat.assink2016) mfor_fit mfor_CR2 <- vcovCR(mfor_fit, type = "CR2") mfor_CR2 coef_test(mfor_fit, vcov = mfor_CR2, test = c("Satterthwaite", "saddlepoint")) Wald_test(mfor_fit, constraints = constrain_zero(3:4), vcov = mfor_CR2) })pkgs_available <- requireNamespace("metafor", quietly = TRUE) & requireNamespace("metadat", quietly = TRUE) if (pkgs_available) withAutoprint({ library(metafor) data(dat.assink2016, package = "metadat") mfor_fit <- rma.mv(yi ~ year + deltype, V = vi, random = ~ 1 | study / esid, data = dat.assink2016) mfor_fit mfor_CR2 <- vcovCR(mfor_fit, type = "CR2") mfor_CR2 coef_test(mfor_fit, vcov = mfor_CR2, test = c("Satterthwaite", "saddlepoint")) Wald_test(mfor_fit, constraints = constrain_zero(3:4), vcov = mfor_CR2) })
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from a
rma.uni object.
## S3 method for class 'rma.uni' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)## S3 method for class 'rma.uni' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Expression or vector indicating which observations
belong to the same cluster. Required for |
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Optional logical indicating whether the weights used in
fitting the model are inverse-variance. If not specified, |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
pkgs_available <- requireNamespace("metafor", quietly = TRUE) & requireNamespace("metadat", quietly = TRUE) if (pkgs_available) withAutoprint({ library(metafor) data(dat.assink2016, package = "metadat") mfor_fit <- rma.uni(yi ~ year + deltype, vi = vi, data = dat.assink2016) mfor_fit mfor_CR2 <- vcovCR(mfor_fit, type = "CR2", cluster = dat.assink2016$study) mfor_CR2 coef_test(mfor_fit, vcov = mfor_CR2, test = c("Satterthwaite", "saddlepoint")) Wald_test(mfor_fit, constraints = constrain_zero(2:4), vcov = mfor_CR2) })pkgs_available <- requireNamespace("metafor", quietly = TRUE) & requireNamespace("metadat", quietly = TRUE) if (pkgs_available) withAutoprint({ library(metafor) data(dat.assink2016, package = "metadat") mfor_fit <- rma.uni(yi ~ year + deltype, vi = vi, data = dat.assink2016) mfor_fit mfor_CR2 <- vcovCR(mfor_fit, type = "CR2", cluster = dat.assink2016$study) mfor_CR2 coef_test(mfor_fit, vcov = mfor_CR2, test = c("Satterthwaite", "saddlepoint")) Wald_test(mfor_fit, constraints = constrain_zero(2:4), vcov = mfor_CR2) })
vcovCR returns a sandwich estimate of the variance-covariance matrix
of a set of regression coefficient estimates from a
robu object.
## S3 method for class 'robu' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)## S3 method for class 'robu' vcovCR(obj, cluster, type, target, inverse_var, form = "sandwich", ...)
obj |
Fitted model for which to calculate the variance-covariance matrix |
cluster |
Optional expression or vector indicating which observations
belong to the same cluster. If not specified, will be set to the
|
type |
Character string specifying which small-sample adjustment should
be used, with available options |
target |
Optional matrix or vector describing the working
variance-covariance model used to calculate the |
inverse_var |
Optional logical indicating whether the weights used in
fitting the model are inverse-variance. If not specified, |
form |
Controls the form of the returned matrix. The default
|
... |
Additional arguments available for some classes of objects. |
An object of class c("vcovCR","clubSandwich"), which consists
of a matrix of the estimated variance of and covariances between the
regression coefficient estimates.
if (requireNamespace("robumeta", quietly = TRUE)) withAutoprint({ library(robumeta) data(hierdat) robu_fit <- robu(effectsize ~ binge + followup + sreport + age, data = hierdat, studynum = studyid, var.eff.size = var, modelweights = "HIER") robu_fit robu_CR2 <- vcovCR(robu_fit, type = "CR2") robu_CR2 coef_test(robu_fit, vcov = robu_CR2, test = c("Satterthwaite", "saddlepoint")) Wald_test(robu_fit, constraints = constrain_zero(c(2,4)), vcov = robu_CR2) Wald_test(robu_fit, constraints = constrain_zero(2:5), vcov = robu_CR2) })if (requireNamespace("robumeta", quietly = TRUE)) withAutoprint({ library(robumeta) data(hierdat) robu_fit <- robu(effectsize ~ binge + followup + sreport + age, data = hierdat, studynum = studyid, var.eff.size = var, modelweights = "HIER") robu_fit robu_CR2 <- vcovCR(robu_fit, type = "CR2") robu_CR2 coef_test(robu_fit, vcov = robu_CR2, test = c("Satterthwaite", "saddlepoint")) Wald_test(robu_fit, constraints = constrain_zero(c(2,4)), vcov = robu_CR2) Wald_test(robu_fit, constraints = constrain_zero(2:5), vcov = robu_CR2) })
Wald_test reports Wald-type tests of linear contrasts from a fitted
linear regression model, using a sandwich estimator for the
variance-covariance matrix and a small sample correction for the p-value.
Several different small-sample corrections are available.
Wald_test( obj, constraints, vcov, null_constant = 0, test = "HTZ", tidy = FALSE, adjustment_method = "none", ... )Wald_test( obj, constraints, vcov, null_constant = 0, test = "HTZ", tidy = FALSE, adjustment_method = "none", ... )
obj |
Fitted model for which to calculate Wald tests. |
constraints |
constraint or list of multiple constraints to test. See details and examples. |
vcov |
Variance covariance matrix estimated using |
null_constant |
vector of null values or list of such vectors for each
set of constraints to test. For a single |
test |
Character vector specifying which small-sample correction(s) to
calculate. The following corrections are available: |
tidy |
Logical value controlling whether to tidy the test results. If
|
adjustment_method |
A character string indicating a multiple comparisons
correction to apply to p-values in instances where multiple tests are run.
Possible options are from |
... |
Further arguments passed to |
Constraints can be specified directly as q X p matrices or
indirectly through constrain_equal,
constrain_zero, or constrain_pairwise. By
default, each constraint will be tested against the null hypothesis that it
equal to a zero vector. Non-zero values for null-hypotheses can be
specified using the null_constant argument.
A list of test results.
vcovCR, constrain_equal,
constrain_zero, constrain_pairwise
if (requireNamespace("carData", quietly = TRUE)) withAutoprint({ data(Duncan, package = "carData") Duncan$cluster <- sample(LETTERS[1:8], size = nrow(Duncan), replace = TRUE) Duncan_fit <- lm(prestige ~ 0 + type + income + type:income + type:education, data=Duncan) # Note that type:income terms are interactions because main effect of income is included # but type:education terms are separate slopes for each unique level of type # Test equality of intercepts Wald_test(Duncan_fit, constraints = constrain_equal(1:3), vcov = "CR2", cluster = Duncan$cluster) # Test equality of type-by-education slopes Wald_test(Duncan_fit, constraints = constrain_equal(":education", reg_ex = TRUE), vcov = "CR2", cluster = Duncan$cluster) # Pairwise comparisons of type-by-education slopes Wald_test(Duncan_fit, constraints = constrain_pairwise(":education", reg_ex = TRUE), vcov = "CR2", cluster = Duncan$cluster) # Test type-by-income interactions Wald_test(Duncan_fit, constraints = constrain_zero(":income", reg_ex = TRUE), vcov = "CR2", cluster = Duncan$cluster) # Pairwise comparisons of type-by-income interactions Wald_test(Duncan_fit, constraints = constrain_pairwise(":income", reg_ex = TRUE, with_zero = TRUE), vcov = "CR2", cluster = Duncan$cluster) # Pairwise comparisons of type-by-education slopes, with two tests and # multiple comparisons p-value adjustment Wald_test(Duncan_fit, constraints = constrain_pairwise(":education", reg_ex = TRUE), vcov = "CR2", cluster = Duncan$cluster, test = c("HTZ","chi-sq"), adjustment_method = "holm") })if (requireNamespace("carData", quietly = TRUE)) withAutoprint({ data(Duncan, package = "carData") Duncan$cluster <- sample(LETTERS[1:8], size = nrow(Duncan), replace = TRUE) Duncan_fit <- lm(prestige ~ 0 + type + income + type:income + type:education, data=Duncan) # Note that type:income terms are interactions because main effect of income is included # but type:education terms are separate slopes for each unique level of type # Test equality of intercepts Wald_test(Duncan_fit, constraints = constrain_equal(1:3), vcov = "CR2", cluster = Duncan$cluster) # Test equality of type-by-education slopes Wald_test(Duncan_fit, constraints = constrain_equal(":education", reg_ex = TRUE), vcov = "CR2", cluster = Duncan$cluster) # Pairwise comparisons of type-by-education slopes Wald_test(Duncan_fit, constraints = constrain_pairwise(":education", reg_ex = TRUE), vcov = "CR2", cluster = Duncan$cluster) # Test type-by-income interactions Wald_test(Duncan_fit, constraints = constrain_zero(":income", reg_ex = TRUE), vcov = "CR2", cluster = Duncan$cluster) # Pairwise comparisons of type-by-income interactions Wald_test(Duncan_fit, constraints = constrain_pairwise(":income", reg_ex = TRUE, with_zero = TRUE), vcov = "CR2", cluster = Duncan$cluster) # Pairwise comparisons of type-by-education slopes, with two tests and # multiple comparisons p-value adjustment Wald_test(Duncan_fit, constraints = constrain_pairwise(":education", reg_ex = TRUE), vcov = "CR2", cluster = Duncan$cluster, test = c("HTZ","chi-sq"), adjustment_method = "holm") })