Function to obtain cluster-robust tests and confidence intervals (also known as robust variance estimation) of the model coefficients for objects of class "rma".

robust(x, cluster, ...)

# S3 method for class 'rma.uni'
robust(x, cluster, adjust=TRUE, clubSandwich=FALSE, digits, ...)
# S3 method for class 'rma.mv'
robust(x, cluster, adjust=TRUE, clubSandwich=FALSE, digits, ...)

Arguments

x

an object of class "rma.uni" or "rma.mv".

cluster

vector to specify the clustering variable to use for constructing the sandwich estimator of the variance-covariance matrix.

adjust

logical to specify whether a small-sample correction should be applied to the variance-covariance matrix.

clubSandwich

logical to specify whether the clubSandwich package should be used to obtain the cluster-robust tests and confidence intervals.

digits

optional integer to specify the number of decimal places to which the printed results should be rounded. If unspecified, the default is to take the value from the object.

...

other arguments.

Details

The function constructs a cluster-robust estimate of the variance-covariance matrix of the model coefficients based on a sandwich-type estimator and then computes tests and confidence intervals of the model coefficients. This function will often be part of a general workflow for meta-analyses involving complex dependency structures as described here.

By default, tests of individual coefficients and confidence intervals are based on a t-distribution with \(n-p\) degrees of freedom, while the omnibus test uses an F-distribution with \(m\) and \(n-p\) degrees of freedom, where \(n\) is the number of clusters, \(p\) denotes the total number of model coefficients (including the intercept if it is present), and \(m\) denotes the number of coefficients tested by the omnibus test. This is sometimes called the ‘residual’ method for approximating the (denominator) degrees of freedom.

When adjust=TRUE (the default), the cluster-robust estimate of the variance-covariance matrix is multiplied by the factor \(n/(n-p)\), which serves as a small-sample adjustment that tends to improve the performance of the method when the number of clusters is small. This is sometimes called the ‘CR1’ adjustment/estimator (in contrast to ‘CR0’ when adjust=FALSE).

For an even better small-sample adjustment, one can set clubSandwich=TRUE in which case the clubSandwich package is used to obtain the cluster-robust tests and confidence intervals. The variance-covariance matrix of the model coefficients is then estimated using the ‘bias-reduced linearization’ adjustment proposed by Bell and McCaffrey (2002) and further developed in Tipton (2015) and Pustejovsky and Tipton (2018). This is sometimes called the ‘CR2’ adjustment/estimator. The degrees of freedom of the t-tests are then estimated using a Satterthwaite approximation. F-tests are then based on an approximate Hotelling's T-squared reference distribution, with denominator degrees of freedom estimated using a method by Zhang (2012, 2013), as further described in Tipton and Pustejovky (2015).

Value

An object of class "robust.rma". The object is a list containing the following components:

beta

estimated coefficients of the model.

se

robust standard errors of the coefficients.

zval

test statistics of the coefficients.

pval

corresponding p-values.

ci.lb

lower bound of the confidence intervals for the coefficients.

ci.ub

upper bound of the confidence intervals for the coefficients.

vb

robust variance-covariance matrix of the estimated coefficients.

QM

test statistic of the omnibus test of moderators.

QMp

corresponding p-value.

...

some additional elements/values.

The results are formatted and printed with the print.rma.uni and print.rma.mv functions (depending on the type of model).

Predicted/fitted values based on "robust.rma" objects can be obtained with the predict function. Tests for sets of model coefficients or linear combinations thereof can be obtained with the anova function.

Note

The variable specified via cluster is assumed to be of the same length as the data originally passed to the rma.uni or rma.mv functions (and if the data argument was used in the original model fit, then the variable will be searched for within this data frame first). Any subsetting and removal of studies with missing values that was applied during the model fitting is also automatically applied to the variable specified via the cluster argument.

The idea of the robust (sandwich-type) estimator for models with unspecified heteroscedasticity can be traced back to Eicker (1967), Huber (1967), and White (1980, 1984). Hence, the method in general is often referred to as the Eicker-Huber-White method. Some small-sample improvements to the method are described by MacKinnon and White (1985). The extension to the cluster-robust estimator can be found in Froot (1989) and Williams (2000), which is also related to the GEE approach by Liang and Zeger (1986). Cameron and Miller (2015) provide an extensive overview of cluster-robust methods. Sidik and Jonkman (2005, 2006) introduced robust methods in the meta-analytic context for standard random/mixed-effects models. The use of cluster-robust methods for multivariate/multilevel meta-analytic models was introduced by Hedges, Tipton, and Johnson (2010).

References

Bell, R. M., & McCaffry, D. F. (2002). Bias reduction in standard errors for linear regression with multi-stage samples. Survey Methodology, 28(2), 169–181. https://www150.statcan.gc.ca/n1/en/catalogue/12-001-X20020029058

Cameron, A. C., & Miller, D. L. (2015). A practitioner's guide to cluster-robust inference. Journal of Human Resources, 50(2), 317–372. https://doi.org/10.3368/jhr.50.2.317

Eicker, F. (1967). Limit theorems for regressions with unequal and dependent errors. In L. M. LeCam & J. Neyman (Eds.), Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability (pp. 59–82). Berkeley: University of California Press.

Froot, K. A. (1989). Consistent covariance matrix estimation with cross-sectional dependence and heteroskedasticity in financial data. Journal of Financial and Quantitative Analysis, 24(3), 333–355. https://doi.org/10.2307/2330815

Hedges, L. V., Tipton, E., & Johnson, M. C. (2010). Robust variance estimation in meta-regression with dependent effect size estimates. Research Synthesis Methods, 1(1), 39–65. https://doi.org/10.1002/jrsm.5

Huber, P. (1967). The behavior of maximum-likelihood estimates under nonstandard conditions. In L. M. LeCam & J. Neyman (Eds.), Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability (pp. 221–233). University of California Press.

Liang, K. Y., & Zeger, S. L. (1986). Longitudinal data analysis using generalized linear models. Biometrika, 73(1), 13–22. https://doi.org/10.1093/biomet/73.1.13

MacKinnon, J. G., & White, H. (1985). Some heteroskedasticity-consistent covariance matrix estimators with improved finite sample properties. Journal of Econometrics, 29(3), 305–325. https://doi.org/10.1016/0304-4076(85)90158-7

Tipton, E. (2015). Small sample adjustments for robust variance estimation with meta-regression. Psychological Methods, 20(3), 375–393. https://doi.org/10.1037/met0000011

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. https://doi.org/10.3102/1076998615606099

Sidik, K., & Jonkman, J. N. (2005). A note on variance estimation in random effects meta-regression. Journal of Biopharmaceutical Statistics, 15(5), 823–838. https://doi.org/10.1081/BIP-200067915

Sidik, K., & Jonkman, J. N. (2006). Robust variance estimation for random effects meta-analysis. Computational Statistics & Data Analysis, 50(12), 3681–3701. https://doi.org/10.1016/j.csda.2005.07.019

White, H. (1980). A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica, 48(4), 817–838. https://doi.org/10.2307/1912934

White, H. (1984). Asymptotic theory for econometricians. Orlando, FL: Academic Press.

Williams, R. L. (2000). A note on robust variance estimation for cluster-correlated data. Biometrics, 56(2), 645–646. https://doi.org/10.1111/j.0006-341x.2000.00645.x

Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3), 1–48. https://doi.org/10.18637/jss.v036.i03

Zhang, J.-T. (2012). An approximate Hotelling T2-test for heteroscedastic one-way MANOVA. Open Journal of Statistics, 2(1), 1–11. https://doi.org/10.4236/ojs.2012.21001

Zhang, J.-T. (2013). Tests of linear hypotheses in the ANOVA under heteroscedasticity. International Journal of Advanced Statistics and Probability, 1, 9–24. https://doi.org/10.14419/ijasp.v1i2.908

See also

rma.uni and rma.mv for functions to fit models for which cluster-robust tests and confidence intervals can be obtained.

Examples

############################################################################

### copy data from Bangert-Drowns et al. (2004) into 'dat'
dat <- dat.bangertdrowns2004

### fit random-effects model
res <- rma(yi, vi, data=dat)
res
#> 
#> Random-Effects Model (k = 48; tau^2 estimator: REML)
#> 
#> tau^2 (estimated amount of total heterogeneity): 0.0499 (SE = 0.0197)
#> tau (square root of estimated tau^2 value):      0.2235
#> I^2 (total heterogeneity / total variability):   58.37%
#> H^2 (total variability / sampling variability):  2.40
#> 
#> Test for Heterogeneity:
#> Q(df = 47) = 107.1061, p-val < .0001
#> 
#> Model Results:
#> 
#> estimate      se    zval    pval   ci.lb   ci.ub      
#>   0.2219  0.0460  4.8209  <.0001  0.1317  0.3122  *** 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 

### use cluster-robust inference methods
robust(res, cluster=id)
#> 
#> Random-Effects Model (k = 48; tau^2 estimator: REML)
#> 
#> tau^2 (estimated amount of total heterogeneity): 0.0499 (SE = 0.0197)
#> tau (square root of estimated tau^2 value):      0.2235
#> I^2 (total heterogeneity / total variability):   58.37%
#> H^2 (total variability / sampling variability):  2.40
#> 
#> Test for Heterogeneity:
#> Q(df = 47) = 107.1061, p-val < .0001
#> 
#> Number of estimates:   48
#> Number of clusters:    48
#> Estimates per cluster: 1
#> 
#> Model Results:
#> 
#> estimate      se¹    tval¹  df¹    pval¹   ci.lb¹   ci.ub¹      
#>   0.2219  0.0460   4.8283   47   <.0001   0.1295   0.3144   *** 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> 1) results based on cluster-robust inference (var-cov estimator: CR1,
#>    approx t-test and confidence interval, df: residual method)
#> 

### use methods from the clubSandwich package
robust(res, cluster=id, clubSandwich=TRUE)
#> 
#> Random-Effects Model (k = 48; tau^2 estimator: REML)
#> 
#> tau^2 (estimated amount of total heterogeneity): 0.0499 (SE = 0.0197)
#> tau (square root of estimated tau^2 value):      0.2235
#> I^2 (total heterogeneity / total variability):   58.37%
#> H^2 (total variability / sampling variability):  2.40
#> 
#> Test for Heterogeneity:
#> Q(df = 47) = 107.1061, p-val < .0001
#> 
#> Number of estimates:   48
#> Number of clusters:    48
#> Estimates per cluster: 1
#> 
#> Model Results:
#> 
#> estimate      se¹    tval¹     df¹    pval¹   ci.lb¹   ci.ub¹      
#>   0.2219  0.0460   4.8280   41.04   <.0001   0.1291   0.3148   *** 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> 1) results based on cluster-robust inference (var-cov estimator: CR2,
#>    approx t-test and confidence interval, df: Satterthwaite approx)
#> 

### fit meta-regression model
res <- rma(yi, vi, mods = ~ length, data=dat)
#> Warning: 2 studies with NAs omitted from model fitting.
res
#> 
#> Mixed-Effects Model (k = 46; tau^2 estimator: REML)
#> 
#> tau^2 (estimated amount of residual heterogeneity):     0.0441 (SE = 0.0188)
#> tau (square root of estimated tau^2 value):             0.2100
#> I^2 (residual heterogeneity / unaccounted variability): 55.26%
#> H^2 (unaccounted variability / sampling variability):   2.24
#> R^2 (amount of heterogeneity accounted for):            5.08%
#> 
#> Test for Residual Heterogeneity:
#> QE(df = 44) = 96.2810, p-val < .0001
#> 
#> Test of Moderators (coefficient 2):
#> QM(df = 1) = 4.2266, p-val = 0.0398
#> 
#> Model Results:
#> 
#>          estimate      se    zval    pval    ci.lb   ci.ub    
#> intrcpt    0.0692  0.0825  0.8384  0.4018  -0.0925  0.2309    
#> length     0.0149  0.0073  2.0559  0.0398   0.0007  0.0292  * 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 

### use cluster-robust inference methods
robust(res, cluster=id)
#> 
#> Mixed-Effects Model (k = 46; tau^2 estimator: REML)
#> 
#> tau^2 (estimated amount of residual heterogeneity):     0.0441 (SE = 0.0188)
#> tau (square root of estimated tau^2 value):             0.2100
#> I^2 (residual heterogeneity / unaccounted variability): 55.26%
#> H^2 (unaccounted variability / sampling variability):   2.24
#> R^2 (amount of heterogeneity accounted for):            5.08%
#> 
#> Test for Residual Heterogeneity:
#> QE(df = 44) = 96.2810, p-val < .0001
#> 
#> Number of estimates:   46
#> Number of clusters:    46
#> Estimates per cluster: 1
#> 
#> Test of Moderators (coefficient 2):¹
#> F(df1 = 1, df2 = 44) = 6.4843, p-val = 0.0145
#> 
#> Model Results:
#> 
#>          estimate      se¹    tval¹  df¹    pval¹    ci.lb¹   ci.ub¹    
#> intrcpt    0.0692  0.0651   1.0628   44   0.2937   -0.0620   0.2004     
#> length     0.0149  0.0059   2.5464   44   0.0145    0.0031   0.0268   * 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> 1) results based on cluster-robust inference (var-cov estimator: CR1,
#>    approx t/F-tests and confidence intervals, df: residual method)
#> 

### use methods from the clubSandwich package
robust(res, cluster=id, clubSandwich=TRUE)
#> 
#> Mixed-Effects Model (k = 46; tau^2 estimator: REML)
#> 
#> tau^2 (estimated amount of residual heterogeneity):     0.0441 (SE = 0.0188)
#> tau (square root of estimated tau^2 value):             0.2100
#> I^2 (residual heterogeneity / unaccounted variability): 55.26%
#> H^2 (unaccounted variability / sampling variability):   2.24
#> R^2 (amount of heterogeneity accounted for):            5.08%
#> 
#> Test for Residual Heterogeneity:
#> QE(df = 44) = 96.2810, p-val < .0001
#> 
#> Number of estimates:   46
#> Number of clusters:    46
#> Estimates per cluster: 1
#> 
#> Test of Moderators (coefficient 2):¹
#> F(df1 = 1, df2 = 19.16) = 6.3961, p-val = 0.0204
#> 
#> Model Results:
#> 
#>          estimate      se¹    tval¹     df¹    pval¹    ci.lb¹   ci.ub¹    
#> intrcpt    0.0692  0.0655   1.0560   16.51   0.3062   -0.0694   0.2077     
#> length     0.0149  0.0059   2.5290   19.16   0.0204    0.0026   0.0273   * 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> 1) results based on cluster-robust inference (var-cov estimator: CR2,
#>    approx t/F-tests and confidence intervals, df: Satterthwaite approx)
#> 

############################################################################

### copy data from Konstantopoulos (2011) into 'dat'
dat <- dat.konstantopoulos2011

### fit multilevel random-effects model
res <- rma.mv(yi, vi, random = ~ 1 | district/school, data=dat)
res
#> 
#> Multivariate Meta-Analysis Model (k = 56; method: REML)
#> 
#> Variance Components:
#> 
#>             estim    sqrt  nlvls  fixed           factor 
#> sigma^2.1  0.0651  0.2551     11     no         district 
#> sigma^2.2  0.0327  0.1809     56     no  district/school 
#> 
#> Test for Heterogeneity:
#> Q(df = 55) = 578.8640, p-val < .0001
#> 
#> Model Results:
#> 
#> estimate      se    zval    pval   ci.lb   ci.ub    
#>   0.1847  0.0846  2.1845  0.0289  0.0190  0.3504  * 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 

### use cluster-robust inference methods
robust(res, cluster=district)
#> 
#> Multivariate Meta-Analysis Model (k = 56; method: REML)
#> 
#> Variance Components:
#> 
#>             estim    sqrt  nlvls  fixed           factor 
#> sigma^2.1  0.0651  0.2551     11     no         district 
#> sigma^2.2  0.0327  0.1809     56     no  district/school 
#> 
#> Test for Heterogeneity:
#> Q(df = 55) = 578.8640, p-val < .0001
#> 
#> Number of estimates:   56
#> Number of clusters:    11
#> Estimates per cluster: 3-11 (mean: 5.09, median: 4)
#> 
#> Model Results:
#> 
#> estimate      se¹    tval¹  df¹    pval¹    ci.lb¹   ci.ub¹    
#>   0.1847  0.0845   2.1859   10   0.0537   -0.0036   0.3730   . 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> 1) results based on cluster-robust inference (var-cov estimator: CR1,
#>    approx t-test and confidence interval, df: residual method)
#> 

### use methods from the clubSandwich package
robust(res, cluster=district, clubSandwich=TRUE)
#> 
#> Multivariate Meta-Analysis Model (k = 56; method: REML)
#> 
#> Variance Components:
#> 
#>             estim    sqrt  nlvls  fixed           factor 
#> sigma^2.1  0.0651  0.2551     11     no         district 
#> sigma^2.2  0.0327  0.1809     56     no  district/school 
#> 
#> Test for Heterogeneity:
#> Q(df = 55) = 578.8640, p-val < .0001
#> 
#> Number of estimates:   56
#> Number of clusters:    11
#> Estimates per cluster: 3-11 (mean: 5.09, median: 4)
#> 
#> Model Results:
#> 
#> estimate      se¹    tval¹    df¹    pval¹    ci.lb¹   ci.ub¹    
#>   0.1847  0.0845   2.1870   9.86   0.0540   -0.0038   0.3733   . 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> 1) results based on cluster-robust inference (var-cov estimator: CR2,
#>    approx t-test and confidence interval, df: Satterthwaite approx)
#> 

############################################################################

### copy data from Berkey et al. (1998) into 'dat'
dat <- dat.berkey1998

### variables v1i and v2i correspond to the 2x2 var-cov matrices of the studies;
### so use these variables to construct the V matrix (note: since v1i and v2i are
### var-cov matrices and not correlation matrices, set vi=1 for all rows)
V <- vcalc(vi=1, cluster=author, rvars=c(v1i, v2i), data=dat)

### fit multivariate model
res <- rma.mv(yi, V, mods = ~ 0 + outcome, random = ~ outcome | trial, struct="UN", data=dat)
res
#> 
#> Multivariate Meta-Analysis Model (k = 10; method: REML)
#> 
#> Variance Components:
#> 
#> outer factor: trial   (nlvls = 5)
#> inner factor: outcome (nlvls = 2)
#> 
#>             estim    sqrt  k.lvl  fixed  level 
#> tau^2.1    0.0327  0.1807      5     no     AL 
#> tau^2.2    0.0117  0.1083      5     no     PD 
#> 
#>     rho.AL  rho.PD    AL  PD 
#> AL       1             -   5 
#> PD  0.6088       1    no   - 
#> 
#> Test for Residual Heterogeneity:
#> QE(df = 8) = 128.2267, p-val < .0001
#> 
#> Test of Moderators (coefficients 1:2):
#> QM(df = 2) = 108.8616, p-val < .0001
#> 
#> Model Results:
#> 
#>            estimate      se     zval    pval    ci.lb    ci.ub      
#> outcomeAL   -0.3392  0.0879  -3.8589  0.0001  -0.5115  -0.1669  *** 
#> outcomePD    0.3534  0.0588   6.0057  <.0001   0.2381   0.4688  *** 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 

### use cluster-robust inference methods
robust(res, cluster=trial)
#> 
#> Multivariate Meta-Analysis Model (k = 10; method: REML)
#> 
#> Variance Components:
#> 
#> outer factor: trial   (nlvls = 5)
#> inner factor: outcome (nlvls = 2)
#> 
#>             estim    sqrt  k.lvl  fixed  level 
#> tau^2.1    0.0327  0.1807      5     no     AL 
#> tau^2.2    0.0117  0.1083      5     no     PD 
#> 
#>     rho.AL  rho.PD    AL  PD 
#> AL       1             -   5 
#> PD  0.6088       1    no   - 
#> 
#> Test for Residual Heterogeneity:
#> QE(df = 8) = 128.2267, p-val < .0001
#> 
#> Number of estimates:   10
#> Number of clusters:    5
#> Estimates per cluster: 2
#> 
#> Test of Moderators (coefficients 1:2):¹
#> F(df1 = 2, df2 = 3) = 41.6138, p-val = 0.0065
#> 
#> Model Results:
#> 
#>            estimate      se¹     tval¹  df¹    pval¹    ci.lb¹    ci.ub¹    
#> outcomeAL   -0.3392  0.1010   -3.3597    3   0.0437   -0.6605   -0.0179   * 
#> outcomePD    0.3534  0.0675    5.2338    3   0.0136    0.1385    0.5683   * 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> 1) results based on cluster-robust inference (var-cov estimator: CR1,
#>    approx t/F-tests and confidence intervals, df: residual method)
#> 

### use methods from the clubSandwich package
robust(res, cluster=trial, clubSandwich=TRUE)
#> 
#> Multivariate Meta-Analysis Model (k = 10; method: REML)
#> 
#> Variance Components:
#> 
#> outer factor: trial   (nlvls = 5)
#> inner factor: outcome (nlvls = 2)
#> 
#>             estim    sqrt  k.lvl  fixed  level 
#> tau^2.1    0.0327  0.1807      5     no     AL 
#> tau^2.2    0.0117  0.1083      5     no     PD 
#> 
#>     rho.AL  rho.PD    AL  PD 
#> AL       1             -   5 
#> PD  0.6088       1    no   - 
#> 
#> Test for Residual Heterogeneity:
#> QE(df = 8) = 128.2267, p-val < .0001
#> 
#> Number of estimates:   10
#> Number of clusters:    5
#> Estimates per cluster: 2
#> 
#> Test of Moderators (coefficients 1:2):¹
#> F(df1 = 2, df2 = 2.76) = 40.7911, p-val = 0.0089
#> 
#> Model Results:
#> 
#>            estimate      se¹     tval¹    df¹    pval¹    ci.lb¹    ci.ub¹     
#> outcomeAL   -0.3392  0.0892   -3.8031   3.81   0.0208   -0.5919   -0.0866    * 
#> outcomePD    0.3534  0.0582    6.0776   3.79   0.0044    0.1884    0.5185   ** 
#> 
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> 1) results based on cluster-robust inference (var-cov estimator: CR2,
#>    approx t/F-tests and confidence intervals, df: Satterthwaite approx)
#> 

############################################################################