Binary files /tmp/tmp2nyyhvj_/NDQfUx3Xa8/r-cran-clubsandwich-0.5.7/build/partial.rdb and /tmp/tmp2nyyhvj_/L7P4owpT8u/r-cran-clubsandwich-0.5.8/build/partial.rdb differ diff -Nru r-cran-clubsandwich-0.5.7/debian/changelog r-cran-clubsandwich-0.5.8/debian/changelog --- r-cran-clubsandwich-0.5.7/debian/changelog 2022-06-22 09:50:11.000000000 +0000 +++ r-cran-clubsandwich-0.5.8/debian/changelog 2022-08-19 12:42:14.000000000 +0000 @@ -1,3 +1,10 @@ +r-cran-clubsandwich (0.5.8-1) unstable; urgency=medium + + * Team upload. + * New upstream version 0.5.8 + + -- Nilesh Patra Fri, 19 Aug 2022 18:12:14 +0530 + r-cran-clubsandwich (0.5.7-1) unstable; urgency=medium * New upstream version diff -Nru r-cran-clubsandwich-0.5.7/DESCRIPTION r-cran-clubsandwich-0.5.8/DESCRIPTION --- r-cran-clubsandwich-0.5.7/DESCRIPTION 2022-06-15 09:20:16.000000000 +0000 +++ r-cran-clubsandwich-0.5.8/DESCRIPTION 2022-08-15 07:50:01.000000000 +0000 @@ -1,7 +1,7 @@ Package: clubSandwich Title: Cluster-Robust (Sandwich) Variance Estimators with Small-Sample Corrections -Version: 0.5.7 +Version: 0.5.8 Authors@R: person("James", "Pustejovsky", email = "jepusto@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0003-0591-9465")) Description: Provides several cluster-robust variance estimators (i.e., sandwich estimators) for ordinary and weighted least squares linear regression @@ -15,9 +15,10 @@ 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(), ivreg() (from package 'AER'), plm() (from package 'plm'), gls() - and lme() (from 'nlme'), lmer() (from `lme4`), robu() (from 'robumeta'), and - rma.uni() and rma.mv() (from 'metafor'). + objects, glm(), 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'), and rma.uni() + and rma.mv() (from 'metafor'). URL: http://jepusto.github.io/clubSandwich/ BugReports: https://github.com/jepusto/clubSandwich/issues Depends: R (>= 3.0.0) @@ -27,13 +28,13 @@ Imports: stats, sandwich Suggests: Formula, knitr, carData, geepack, metafor, metadat, robumeta, nlme, mlmRev, AER, plm (>= 1.6-4), Matrix, lme4, zoo, testthat, - rmarkdown, covr + rmarkdown, covr, ivreg RoxygenNote: 7.1.2 Encoding: UTF-8 Language: en-US NeedsCompilation: no -Packaged: 2022-06-14 16:42:53 UTC; jepus +Packaged: 2022-08-13 14:34:57 UTC; jepusto Author: James Pustejovsky [aut, cre] () Maintainer: James Pustejovsky Repository: CRAN -Date/Publication: 2022-06-15 09:20:16 UTC +Date/Publication: 2022-08-15 07:50:01 UTC diff -Nru r-cran-clubsandwich-0.5.7/inst/doc/meta-analysis-with-CRVE.html r-cran-clubsandwich-0.5.8/inst/doc/meta-analysis-with-CRVE.html --- r-cran-clubsandwich-0.5.7/inst/doc/meta-analysis-with-CRVE.html 2022-06-14 16:42:47.000000000 +0000 +++ r-cran-clubsandwich-0.5.8/inst/doc/meta-analysis-with-CRVE.html 2022-08-13 14:34:48.000000000 +0000 @@ -12,7 +12,7 @@ - + Meta-analysis with cluster-robust variance estimation @@ -334,7 +334,7 @@

Meta-analysis with cluster-robust variance estimation

James E. Pustejovsky

-

2022-06-14

+

2022-08-13

diff -Nru r-cran-clubsandwich-0.5.7/inst/doc/panel-data-CRVE.html r-cran-clubsandwich-0.5.8/inst/doc/panel-data-CRVE.html --- r-cran-clubsandwich-0.5.7/inst/doc/panel-data-CRVE.html 2022-06-14 16:42:50.000000000 +0000 +++ r-cran-clubsandwich-0.5.8/inst/doc/panel-data-CRVE.html 2022-08-13 14:34:52.000000000 +0000 @@ -12,7 +12,7 @@ - + Cluster-robust standard errors and hypothesis tests in panel data models @@ -334,7 +334,7 @@

Cluster-robust standard errors and hypothesis tests in panel data models

James E. Pustejovsky

-

2022-06-14

+

2022-08-13

@@ -470,15 +470,14 @@ fixed effects before estimating the effect of legal. The clubSandwich package works with fitted plm models too:

-
library(plm)
-
## Warning: package 'plm' was built under R version 4.1.3
-
plm_unweighted <- plm(mrate ~ legal + beertaxa, data = MV_deaths, 
-                      effect = "twoways", index = c("state","year"))
-coef_test(plm_unweighted, vcov = "CR1", cluster = "individual", test = "naive-t")
+
library(plm)
+plm_unweighted <- plm(mrate ~ legal + beertaxa, data = MV_deaths, 
+                      effect = "twoways", index = c("state","year"))
+coef_test(plm_unweighted, vcov = "CR1", cluster = "individual", test = "naive-t")
##     Coef. Estimate   SE t-stat d.f. (naive-t) p-val (naive-t) Sig.
 ##     legal     7.59 2.44  3.108             49         0.00313   **
 ##  beertaxa     3.82 5.14  0.743             49         0.46128
-
coef_test(plm_unweighted, vcov = "CR2", cluster = "individual", test = "Satterthwaite")
+
coef_test(plm_unweighted, vcov = "CR2", cluster = "individual", test = "Satterthwaite")
##     Coef. Estimate   SE t-stat d.f. (Satt) p-val (Satt) Sig.
 ##     legal     7.59 2.51  3.019       24.58      0.00583   **
 ##  beertaxa     3.82 5.27  0.725        5.77      0.49663
@@ -490,15 +489,15 @@ bit if the model is estimated using population weights. We go back to fitting in lm with dummies for all the fixed effects because plm does not handle weighted least squares.

-
lm_weighted <- lm(mrate ~ 0 + legal + beertaxa + factor(state) + factor(year), 
-                  weights = pop, data = MV_deaths)
-coef_test(lm_weighted, vcov = "CR1", 
-          cluster = MV_deaths$state, test = "naive-t")[1:2,]
+
lm_weighted <- lm(mrate ~ 0 + legal + beertaxa + factor(state) + factor(year), 
+                  weights = pop, data = MV_deaths)
+coef_test(lm_weighted, vcov = "CR1", 
+          cluster = MV_deaths$state, test = "naive-t")[1:2,]
##     Coef. Estimate   SE t-stat d.f. (naive-t) p-val (naive-t) Sig.
 ##     legal     7.78 2.01   3.87             49          <0.001  ***
 ##  beertaxa    11.16 4.20   2.66             49          0.0106    *
-
coef_test(lm_weighted, vcov = "CR2", 
-          cluster = MV_deaths$state, test = "Satterthwaite")[1:2,]
+
coef_test(lm_weighted, vcov = "CR2", 
+          cluster = MV_deaths$state, test = "Satterthwaite")[1:2,]
##     Coef. Estimate   SE t-stat d.f. (Satt) p-val (Satt) Sig.
 ##     legal     7.78 2.13   3.64        8.52      0.00588   **
 ##  beertaxa    11.16 4.37   2.55        6.85      0.03854    *
@@ -517,9 +516,9 @@ standardized rate and so a better assumption might be that the error variances are inversely proportional to population size. The following code uses this alternate working model:

-
coef_test(lm_weighted, vcov = "CR2", 
-          cluster = MV_deaths$state, target = 1 / MV_deaths$pop, 
-          test = "Satterthwaite")[1:2,]
+
coef_test(lm_weighted, vcov = "CR2", 
+          cluster = MV_deaths$state, target = 1 / MV_deaths$pop, 
+          test = "Satterthwaite")[1:2,]
##     Coef. Estimate   SE t-stat d.f. (Satt) p-val (Satt) Sig.
 ##     legal     7.78 2.03   3.83       12.64      0.00221   **
 ##  beertaxa    11.16 4.17   2.68        5.06      0.04333    *
@@ -533,14 +532,14 @@ with the regressors, then a more efficient way to estimate \(\gamma,\delta\) is by weighted least squares, with weights based on a random effects model. We still treat the year effects as fixed.

-
plm_random <- plm(mrate ~ 0 + legal + beertaxa + year, data = MV_deaths, 
-                  effect = "individual", index = c("state","year"),
-                  model = "random")
-coef_test(plm_random, vcov = "CR1", test = "naive-t")[1:2,]
+
plm_random <- plm(mrate ~ 0 + legal + beertaxa + year, data = MV_deaths, 
+                  effect = "individual", index = c("state","year"),
+                  model = "random")
+coef_test(plm_random, vcov = "CR1", test = "naive-t")[1:2,]
##     Coef. Estimate   SE t-stat d.f. (naive-t) p-val (naive-t) Sig.
 ##     legal     7.31 2.39  3.054             49         0.00364   **
 ##  beertaxa     3.37 5.11  0.661             49         0.51202
-
coef_test(plm_random, vcov = "CR2", test = "Satterthwaite")[1:2,]
+
coef_test(plm_random, vcov = "CR2", test = "Satterthwaite")[1:2,]
##     Coef. Estimate   SE t-stat d.f. (Satt) p-val (Satt) Sig.
 ##     legal     7.31 2.46  2.966       25.18      0.00652   **
 ##  beertaxa     3.37 5.22  0.647        5.78      0.54258
@@ -576,16 +575,16 @@ effects.

For efficiency, we estimate this specification using weighted least squares (although OLS would be valid too):

-
MV_deaths <- within(MV_deaths, {
-  legal_cent <- legal - tapply(legal, state, mean)[factor(state)]
-  beer_cent <- beertaxa - tapply(beertaxa, state, mean)[factor(state)]
-})
-
-plm_Hausman <- plm(mrate ~ 0 + legal + beertaxa + legal_cent + beer_cent + factor(year), 
-                   data = MV_deaths,
-                   effect = "individual", index = c("state","year"),
-                   model = "random")
-coef_test(plm_Hausman, vcov = "CR2", test = "Satterthwaite")[1:4,]
+
MV_deaths <- within(MV_deaths, {
+  legal_cent <- legal - tapply(legal, state, mean)[factor(state)]
+  beer_cent <- beertaxa - tapply(beertaxa, state, mean)[factor(state)]
+})
+
+plm_Hausman <- plm(mrate ~ 0 + legal + beertaxa + legal_cent + beer_cent + factor(year), 
+                   data = MV_deaths,
+                   effect = "individual", index = c("state","year"),
+                   model = "random")
+coef_test(plm_Hausman, vcov = "CR2", test = "Satterthwaite")[1:4,]
##       Coef. Estimate   SE  t-stat d.f. (Satt) p-val (Satt) Sig.
 ##       legal   -9.180 7.62 -1.2042       24.94       0.2398     
 ##    beertaxa    3.395 9.40  0.3613        6.44       0.7295     
@@ -597,18 +596,18 @@
 robust Wald statistic, then use a \(\chi^2_2\) reference distribution (or
 equivalently, compare a re-scaled Wald statistic to an \(F(2,\infty)\) distribution). The
 Wald_test function reports the latter version:

-
Wald_test(plm_Hausman, 
-          constraints = constrain_zero(c("legal_cent","beer_cent")), 
-          vcov = "CR1", test = "chi-sq")
+
Wald_test(plm_Hausman, 
+          constraints = constrain_zero(c("legal_cent","beer_cent")), 
+          vcov = "CR1", test = "chi-sq")
##    test Fstat df_num df_denom  p_val sig
 ##  chi-sq  2.93      2      Inf 0.0534   .

The test is just shy of significance at the 5% level. If we instead use the CR2 variance estimator and our newly proposed approximate F-test (which is the default in Wald_test), then we get:

-
Wald_test(plm_Hausman, 
-          constraints = constrain_zero(c("legal_cent","beer_cent")), 
-          vcov = "CR2")
+
Wald_test(plm_Hausman, 
+          constraints = constrain_zero(c("legal_cent","beer_cent")), 
+          vcov = "CR2")
##  test Fstat df_num df_denom p_val sig
 ##   HTZ  2.56      2     11.7  0.12

The low degrees of freedom of the test indicate that we’re definitely diff -Nru r-cran-clubsandwich-0.5.7/inst/doc/Wald-tests-in-clubSandwich.html r-cran-clubsandwich-0.5.8/inst/doc/Wald-tests-in-clubSandwich.html --- r-cran-clubsandwich-0.5.7/inst/doc/Wald-tests-in-clubSandwich.html 2022-06-14 16:42:41.000000000 +0000 +++ r-cran-clubsandwich-0.5.8/inst/doc/Wald-tests-in-clubSandwich.html 2022-08-13 14:34:37.000000000 +0000 @@ -12,7 +12,7 @@ - + Wald tests of multiple-constraint null hypotheses @@ -356,27 +356,26 @@

Wald tests of multiple-constraint null hypotheses

James E. Pustejovsky

-

2022-06-14

+

2022-08-13

diff -Nru r-cran-clubsandwich-0.5.7/man/vcovCR.ivreg.Rd r-cran-clubsandwich-0.5.8/man/vcovCR.ivreg.Rd --- r-cran-clubsandwich-0.5.7/man/vcovCR.ivreg.Rd 2022-06-12 18:43:57.000000000 +0000 +++ r-cran-clubsandwich-0.5.8/man/vcovCR.ivreg.Rd 2022-08-12 19:46:50.000000000 +0000 @@ -50,11 +50,18 @@ } \description{ \code{vcovCR} returns a sandwich estimate of the variance-covariance matrix -of a set of regression coefficient estimates from an \code{\link[AER]{ivreg}} object. +of a set of regression coefficient estimates from an ivreg object fitted +from the \CRANpkg{AER} package or the \CRANpkg{ivreg} package. +} +\details{ +For any "ivreg" objects fitted via the \code{\link[ivreg]{ivreg}} + function from the \CRANpkg{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"). } \examples{ - if (requireNamespace("AER", quietly = TRUE)) withAutoprint({ library(AER) @@ -65,11 +72,29 @@ tdiff <- (taxs - tax)/cpi }) - iv_fit <- ivreg(log(packs) ~ log(rprice) + log(rincome) | + iv_fit_AER <- AER::ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), data = Cigs) - vcovCR(iv_fit, cluster = Cigs$state, type = "CR2") - coef_test(iv_fit, vcov = "CR2", cluster = Cigs$state) + 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) }) } diff -Nru r-cran-clubsandwich-0.5.7/MD5 r-cran-clubsandwich-0.5.8/MD5 --- r-cran-clubsandwich-0.5.7/MD5 2022-06-15 09:20:16.000000000 +0000 +++ r-cran-clubsandwich-0.5.8/MD5 2022-08-15 07:50:01.000000000 +0000 @@ -1,6 +1,6 @@ -723cd7eca5b305c69cf7f031ba05aba4 *DESCRIPTION +17b0a78c7c4774ad8385ee53539c148a *DESCRIPTION 07ff211393c01ab2b3534ddd7a601d52 *NAMESPACE -d8bbaa928430b94928701fb477869d97 *NEWS.md +b9cbba1f8657902b4c58c2e5d0c9aec0 *NEWS.md a1c4db8cc7d1e4f2209e81d9049bccac *R/CR-adjustments.R 91d2efb468fdd6227c91cf2300b0044e *R/S3-methods.R 35afb3ddff980418f586fb27433088ca *R/Wald_test.R @@ -11,7 +11,7 @@ bb0d8db3acb8652cbbca2efb0cdfbf39 *R/get_arrays.R 3a3fe1e4f49c85dad4d8e8623414cda1 *R/glm.R 3bdd747acc2b0bbbed0a26ea5f30371f *R/gls.R -8d19ac6d19cc9d12e136963bcf316802 *R/ivreg.R +fb5e175f375601fd969792b3cb8c129b *R/ivreg.R 6a425ebb1f8b54dcd4195ae7b4bf70c6 *R/lm.R 878431f8d11b487b99e73efc9aaa9432 *R/lme.R 236c086171e1e53b30e331b1a72547ad *R/lmer.R @@ -22,7 +22,7 @@ 11d98e5c8cf564f1a2d394303ecfe221 *R/rma-uni.R 337e2a1285271f501d6af62af63b9ed5 *R/robu.R efa27f93c711b2fb47c4e3350ca8d6bf *R/utilities.R -9ad66b75c26e484cc7a3991fb847c8da *build/partial.rdb +fbeb53e2d6d5a3c2e2e3e1da5a252b55 *build/partial.rdb 06bf516326af02fb398fa248a0371a5a *build/vignette.rds 8e550f05dda31b5c40d1817bae9a4363 *data/AchievementAwardsRCT.RData 34360e013336c721ad97d9c3658db4df *data/MortalityRates.RData @@ -30,13 +30,13 @@ 6402f15efe8b3298b8dea6d9811a42ce *data/dropoutPrevention.RData d0902ccacf3defd478c6075e6c00720d *inst/doc/Wald-tests-in-clubSandwich.R 57cb6087c04840ccb9e69e0499cd05c9 *inst/doc/Wald-tests-in-clubSandwich.Rmd -d1364451d07dd96a03d5ac4aeac40d00 *inst/doc/Wald-tests-in-clubSandwich.html +d99d73b5396310ddaf1f42a419d70ddb *inst/doc/Wald-tests-in-clubSandwich.html f5d5ffe95e040995200ee6a8fc756210 *inst/doc/meta-analysis-with-CRVE.R ebc446904160ac4add09b6ecf4e9a603 *inst/doc/meta-analysis-with-CRVE.Rmd -54e643734dd3ce43ace980f760a15ee9 *inst/doc/meta-analysis-with-CRVE.html +2716f1679fbffa58f87b9fb70a48ee04 *inst/doc/meta-analysis-with-CRVE.html fbd774bb057c0dbac28800d79de7a8da *inst/doc/panel-data-CRVE.R 4694bc1c8daeb75adc066f6398ac4d0c *inst/doc/panel-data-CRVE.Rmd -73062da08d62172d0657120f985c9754 *inst/doc/panel-data-CRVE.html +0e7f9095aa6c8871e90a61c1ecabe16b *inst/doc/panel-data-CRVE.html cf2f124f9a215d635cbf9ef33c1c74a6 *man/AchievementAwardsRCT.Rd 9dcfded510ec3a852686250e20c258e6 *man/MortalityRates.Rd 7304421949562f2c392131aaafb1afe0 *man/SATcoaching.Rd @@ -52,7 +52,7 @@ 198f468de2915ab1157789f00487cd23 *man/vcovCR.Rd b89fc3af7660f17a83ebb90877f76ace *man/vcovCR.glm.Rd ef6e66f7286ce0bf4093e79580610575 *man/vcovCR.gls.Rd -af208398bac2c5bb6d2fe65267067e03 *man/vcovCR.ivreg.Rd +641da2ab49f0fa699e13c90318508f6e *man/vcovCR.ivreg.Rd 56e10ed54e4b8cfedf2aad822c90a861 *man/vcovCR.lm.Rd 5a3090b35b0fe3d26d48b7be39194c33 *man/vcovCR.lme.Rd a6508917e5d0fa4e1916d6fe975bbfc6 *man/vcovCR.lmerMod.Rd @@ -62,6 +62,7 @@ b3fb1637bee90baeb2755c88f5dd2895 *man/vcovCR.rma.uni.Rd 461b9e7729343401861bdfa1f9deebf3 *man/vcovCR.robu.Rd 77ee716352046c1d86deb383d1d5497a *tests/testthat.R +7b4feabc94745c032a3b14e1ced26d5c *tests/testthat/test_AER_ivreg.R 420b29172eb62c3e952d5f690fb0ed7e *tests/testthat/test_Wald.R 496dade7481c704545b42e21d53040ff *tests/testthat/test_coef.R 57264754aa780ff807794104686098c3 *tests/testthat/test_conf_int.R @@ -71,7 +72,7 @@ 31083daff7e3904291808f2051f8b244 *tests/testthat/test_ignore_absorption.R d544ad56c93e1f10a7e00c09521f364a *tests/testthat/test_impute_covariance_matrix.R f31d0900fcf58534bc397aa43473f560 *tests/testthat/test_intercept_formulas.R -9afb94344689953f303b6795a4e78721 *tests/testthat/test_ivreg.R +4f654c9b380ea95765496c051607a742 *tests/testthat/test_ivreg_ivreg.R af12f4a8e33cda24ae068e8f88f3cd44 *tests/testthat/test_linear_contrast.R e1463e069b6705b8eb2aebff355d58fa *tests/testthat/test_lm.R 73f82f16b1cac2e57a23a362f1e7f41c *tests/testthat/test_lme_2level.R @@ -80,7 +81,7 @@ 14b742f504dc42d71f27be05d19b0b9f *tests/testthat/test_lmerMod.R 0032bb715b62139c7294f6bb88bd8df1 *tests/testthat/test_mlm.R 617b13bf16fdb40eb0b06457d1e22474 *tests/testthat/test_plm-ID-variables.R -b3b447b54ac32d73740d429d8f906353 *tests/testthat/test_plm-first-differences.R +588b5b0d192ac5155fb007744d367ee4 *tests/testthat/test_plm-first-differences.R 37add75e3bf443e3fb06038e6a014293 *tests/testthat/test_plm-fixed-effects.R 56da486857e8d4504e8aa8c1fe7937da *tests/testthat/test_plm-random-effects.R 4597e0d7dfbedc264c0ea336bc5229be *tests/testthat/test_plm-unbalanced-fixed-effects.R diff -Nru r-cran-clubsandwich-0.5.7/NEWS.md r-cran-clubsandwich-0.5.8/NEWS.md --- r-cran-clubsandwich-0.5.7/NEWS.md 2022-06-14 13:57:15.000000000 +0000 +++ r-cran-clubsandwich-0.5.8/NEWS.md 2022-08-12 20:47:03.000000000 +0000 @@ -1,3 +1,8 @@ +# clubSandwich 0.5.8 + +* Added support for `ivreg::ivreg` objects when estimated by ordinary least squares (support for objects estimated by 2SM and 2SMM is not yet implemented). +* Updated unit tests for `plm::plm()` when `method = "FD"` to account for bug fixes in version 2.6-2 of plm. + # clubSandwich 0.5.7 * Fixed bug in methods for multi-variate multi-level models estimated with lme(). @@ -20,8 +25,8 @@ ## New features -* `Wald_test()` gains an option for test = "Naive-Fp", which uses denominator degrees of freedom equal to the number of clusters minus the number of coefficients in the fitted model. -* `coef_test()` and `conf_int()` gain an option for test = "naive-tp", which uses denominator degrees of freedom equal to the number of clusters minus the number of coefficients in the fitted model. +* `Wald_test()` gains an option for `test = "Naive-Fp"`, which uses denominator degrees of freedom equal to the number of clusters minus the number of coefficients in the fitted model. +* `coef_test()` and `conf_int()` gain an option for `test = "naive-tp"`, which uses denominator degrees of freedom equal to the number of clusters minus the number of coefficients in the fitted model. ## Minor improvements and bug fixes diff -Nru r-cran-clubsandwich-0.5.7/R/ivreg.R r-cran-clubsandwich-0.5.8/R/ivreg.R --- r-cran-clubsandwich-0.5.7/R/ivreg.R 2022-06-12 18:41:28.000000000 +0000 +++ r-cran-clubsandwich-0.5.8/R/ivreg.R 2022-08-12 19:46:50.000000000 +0000 @@ -5,7 +5,8 @@ #' Cluster-robust variance-covariance matrix for an ivreg object. #' #' \code{vcovCR} returns a sandwich estimate of the variance-covariance matrix -#' of a set of regression coefficient estimates from an \code{\link[AER]{ivreg}} object. +#' of a set of regression coefficient estimates from an ivreg object fitted +#' from the \CRANpkg{AER} package or the \CRANpkg{ivreg} package. #' #' @param cluster Expression or vector indicating which observations belong to #' the same cluster. Required for \code{ivreg} objects. @@ -15,6 +16,12 @@ #' diagonal. If not specified, the target is taken to be an identity matrix. #' @param inverse_var Not used for \code{ivreg} objects. #' @inheritParams vcovCR +#' +#' @details For any "ivreg" objects fitted via the \code{\link[ivreg]{ivreg}} +#' function from the \CRANpkg{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"). #' #' @return An object of class \code{c("vcovCR","clubSandwich")}, which consists #' of a matrix of the estimated variance of and covariances between the @@ -24,7 +31,6 @@ #' #' @examples #' -#' #' if (requireNamespace("AER", quietly = TRUE)) withAutoprint({ #' #' library(AER) @@ -35,11 +41,29 @@ #' tdiff <- (taxs - tax)/cpi #' }) #' -#' iv_fit <- ivreg(log(packs) ~ log(rprice) + log(rincome) | +#' iv_fit_AER <- AER::ivreg(log(packs) ~ log(rprice) + log(rincome) | #' log(rincome) + tdiff + I(tax/cpi), data = Cigs) -#' vcovCR(iv_fit, cluster = Cigs$state, type = "CR2") -#' coef_test(iv_fit, vcov = "CR2", cluster = Cigs$state) +#' 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) #' }) #' #' @export @@ -47,6 +71,7 @@ vcovCR.ivreg <- function(obj, cluster, type, target = NULL, inverse_var = FALSE, form = "sandwich", ...) { if (missing(cluster)) stop("You must specify a clustering variable.") if (inverse_var != FALSE) stop("Unfortunately, the inverse_var option is not available for ivreg models.") + if (!is.null(obj$method) && obj$method %in% c("M", "MM")) stop("clubSandwich does not currently support ivreg models estimated using method = 'M' or method = 'MM'.") vcov_CR(obj, cluster = cluster, type = type, target = target, inverse_var = inverse_var, form = form) } diff -Nru r-cran-clubsandwich-0.5.7/tests/testthat/test_AER_ivreg.R r-cran-clubsandwich-0.5.8/tests/testthat/test_AER_ivreg.R --- r-cran-clubsandwich-0.5.7/tests/testthat/test_AER_ivreg.R 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-clubsandwich-0.5.8/tests/testthat/test_AER_ivreg.R 2022-08-12 19:46:50.000000000 +0000 @@ -0,0 +1,239 @@ +context("AER::ivreg objects") +set.seed(20190513) + +skip_if_not_installed("zoo") +skip_if_not_installed("AER") + +library(zoo, quietly=TRUE) +library(AER, quietly=TRUE) + +data("CigarettesSW", package = "AER") + +Cigs <- within(CigarettesSW, { + rprice <- price/cpi + rincome <- income/population/cpi + tdiff <- (taxs - tax)/cpi +}) + +CR_types <- paste0("CR",0:4) + +obj_un <- ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), + data = Cigs) +obj_wt <- ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), + data = Cigs, + weights = population) + +X <- model.matrix(obj_wt, component = "regressors") +Z <- model.matrix(obj_wt, component = "instruments") +y <- log(CigarettesSW$packs) +w <- weights(obj_wt) + +test_that("Basic calculations from ivreg agree for unweighted model.", { + XZ <- model.matrix(obj_un, component = "projected") + ZtZ_inv <- chol2inv(chol(t(Z) %*% Z)) + XZ_check <- Z %*% ZtZ_inv %*% t(Z) %*% X + + expect_equal(XZ, XZ_check, check.attributes=FALSE) + expect_equal(coef(obj_un), lm.fit(XZ, y)$coefficients) + expect_equal (bread(obj_un), chol2inv(chol(t(XZ) %*% XZ)) * nobs(obj_un), check.attributes=FALSE) + + hii <- diag(X %*% chol2inv(chol(t(XZ) %*% XZ)) %*% t(XZ)) + expect_equal(hatvalues(obj_un), hii) + + r <- as.vector(y - X %*% coef(obj_un)) + expect_equal(r, as.vector(residuals_CS(obj_un))) +}) + +test_that("Basic calculations from ivreg agree for weighted model.", { + XZ <- model.matrix(obj_wt, component = "projected") + ZwZ_inv <- chol2inv(chol(t(Z) %*% (w * Z))) + XZ_check <- Z %*% ZwZ_inv %*% t(Z) %*% (w * X) + + expect_equal(XZ, XZ_check, check.attributes=FALSE) + expect_equal(coef(obj_wt), lm.wfit(XZ, y, w)$coefficients) + expect_equal(bread(obj_wt), chol2inv(chol(t(XZ) %*% (w * XZ))) * nobs(obj_wt), check.attributes=FALSE) + + hii <- diag(X%*% chol2inv(chol(t(XZ) %*% (w * XZ))) %*% t(w * XZ)) + expect_false(all(hatvalues(obj_wt) == hii)) # does not agree because hatvalues doesn't work with weighting + + r <- as.vector(y - X %*% coef(obj_wt)) + expect_equal(r, as.vector(residuals_CS(obj_wt))) +}) + +test_that("bread works", { + + expect_true(check_bread(obj_un, cluster = Cigs$state, y = log(Cigs$packs))) + tsls_vcov <- bread(obj_un) * summary(obj_un)$sigma^2 / v_scale(obj_un) + expect_equal(vcov(obj_un), tsls_vcov) + + expect_true(check_bread(obj_wt, cluster = Cigs$state, y = log(Cigs$packs))) + wtsls_vcov <- bread(obj_wt) * summary(obj_wt)$sigma^2 / v_scale(obj_wt) + expect_equal(vcov(obj_wt), wtsls_vcov) +}) + + +test_that("vcovCR options don't matter for CR0", { + expect_error(vcovCR(obj_un, type = "CR0")) + CR0 <- vcovCR(obj_un, cluster = Cigs$state, type = "CR0") + expect_output(print(CR0)) + attr(CR0, "target") <- NULL + attr(CR0, "inverse_var") <- NULL + CR0_A <- vcovCR(obj_un, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population) + attr(CR0_A, "target") <- NULL + attr(CR0_A, "inverse_var") <- NULL + expect_identical(CR0_A, CR0) + CR0_B <- vcovCR(obj_un, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population, inverse_var = FALSE) + attr(CR0_B, "target") <- NULL + attr(CR0_B, "inverse_var") <- NULL + expect_identical(CR0_A, CR0) + + expect_error(vcovCR(obj_un, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population, inverse_var = TRUE)) + + wCR0 <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR0") + attr(wCR0, "target") <- NULL + attr(wCR0, "inverse_var") <- NULL + wCR0_A <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population) + attr(wCR0_A, "target") <- NULL + attr(wCR0_A, "inverse_var") <- NULL + expect_identical(wCR0_A, wCR0) + wCR0_B <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population, inverse_var = FALSE) + attr(wCR0_B, "target") <- NULL + attr(wCR0_B, "inverse_var") <- NULL + expect_identical(wCR0_B, wCR0) + + expect_error(vcovCR(obj_wt, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population, inverse_var = TRUE)) +}) + +test_that("vcovCR options work for CR2", { + + CR2_iv <- vcovCR(obj_un, cluster = Cigs$state, type = "CR2") + expect_equal(vcovCR(obj_un, cluster = Cigs$state, type = "CR2", target = rep(1, nobs(obj_un))), CR2_iv) + expect_false(identical(vcovCR(obj_un, cluster = Cigs$state, type = "CR2", target = 1 / Cigs$population), CR2_iv)) + + wCR2_id <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR2") + expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR2", inverse_var = FALSE), wCR2_id) + expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR2", target = rep(1, nobs(obj_un))), wCR2_id) + expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR2", target = rep(1, nobs(obj_un)), inverse_var = FALSE), wCR2_id) + +}) + +test_that("vcovCR options work for CR4", { + + CR4_not <- vcovCR(obj_un, cluster = Cigs$state, type = "CR4") + expect_identical(vcovCR(obj_un, cluster = Cigs$state, type = "CR4", target = rep(1, nobs(obj_un))), CR4_not) + expect_identical(vcovCR(obj_un, cluster = Cigs$state, type = "CR4", target = rep(1, nobs(obj_un)), inverse_var = FALSE), CR4_not) + expect_false(identical(vcovCR(obj_un, cluster = Cigs$state, type = "CR4", target = 1 / Cigs$population), CR4_not)) + + wCR4_id <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR4") + expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR4", inverse_var = FALSE), wCR4_id) + expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR4", target = rep(1, nobs(obj_wt))), wCR4_id) + expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR4", target = rep(1, nobs(obj_wt)), inverse_var = FALSE), wCR4_id) + +}) + + +test_that("CR2 is target-unbiased", { + expect_true(check_CR(obj_un, vcov = "CR2", cluster = Cigs$state)) + expect_true(check_CR(obj_wt, vcov = "CR2", cluster = Cigs$state)) +}) + +test_that("CR4 is target-unbiased", { + skip("Need to understand target-unbiasedness for ivreg objects.") + expect_true(check_CR(obj_un, vcov = "CR4", cluster = Cigs$state)) + expect_true(check_CR(obj_wt, vcov = "CR4", cluster = Cigs$state)) +}) + +test_that("vcovCR is equivalent to vcovHC (with HC0 or HC1) when clusters are all of size 1", { + library(sandwich, quietly=TRUE) + CR0 <- vcovCR(obj_un, cluster = 1:nobs(obj_un), type = "CR0") + expect_equal(vcovHC(obj_un, type = "HC0"), as.matrix(CR0)) + CR1 <- vcovCR(obj_un, cluster = 1:nobs(obj_un), type = "CR1S") + expect_equal(vcovHC(obj_un, type = "HC1"), as.matrix(CR1)) + CR2 <- vcovCR(obj_un, cluster = 1:nobs(obj_un), type = "CR2") + expect_false(all(vcovHC(obj_un, type = "HC2") == as.matrix(CR2))) + CR3 <- vcovCR(obj_un, cluster = 1:nobs(obj_un), type = "CR3") + expect_false(all(vcovHC(obj_un, type = "HC3") == as.matrix(CR3))) +}) + +test_that("Order doesn't matter.",{ + + check_sort_order(obj_wt, Cigs, "state") + +}) + +test_that("clubSandwich works with dropped observations", { + dat_miss <- Cigs + dat_miss$rincome[sample.int(nrow(Cigs), size = round(nrow(Cigs) / 10))] <- NA + iv_dropped <- update(obj_un, data = dat_miss) + dat_complete <- subset(dat_miss, !is.na(rincome)) + iv_complete <- update(obj_un, data = dat_complete) + + CR_drop <- lapply(CR_types, function(x) vcovCR(iv_dropped, cluster = dat_miss$state, type = x)) + CR_complete <- lapply(CR_types, function(x) vcovCR(iv_complete, cluster = dat_complete$state, type = x)) + expect_equal(CR_drop, CR_complete) + + test_drop <- lapply(CR_types, function(x) coef_test(iv_dropped, vcov = x, cluster = dat_miss$state, test = "All", p_values = FALSE)) + test_complete <- lapply(CR_types, function(x) coef_test(iv_complete, vcov = x, cluster = dat_complete$state, test = "All", p_values = FALSE)) + expect_equal(test_drop, test_complete) +}) + + + +test_that("weight scale doesn't matter", { + + iv_fit_w <- update(obj_un, weights = rep(4, nobs(obj_un))) + + unweighted_fit <- lapply(CR_types, function(x) vcovCR(obj_un, cluster = Cigs$state, type = x)) + weighted_fit <- lapply(CR_types, function(x) vcovCR(iv_fit_w, cluster = Cigs$state, type = x)) + + expect_equal(lapply(unweighted_fit, as.matrix), + lapply(weighted_fit, as.matrix), + tol = 5 * 10^-7) + + target <- 1 + rpois(nrow(Cigs), lambda = 8) + unweighted_fit <- lapply(CR_types, function(x) vcovCR(obj_un, cluster = Cigs$state, type = x, target = target)) + weighted_fit <- lapply(CR_types, function(x) vcovCR(iv_fit_w, cluster = Cigs$state, type = x, target = target * 15)) + + expect_equal(lapply(unweighted_fit, as.matrix), + lapply(weighted_fit, as.matrix), + tol = 5 * 10^-7) + +}) + +test_that("clubSandwich works with weights of zero.", { + + n_Cigs <- nrow(Cigs) + Cigs$awt <- rpois(n_Cigs, lambda = 1.4) + table(Cigs$awt) + + iv_full <- update(obj_un, weights = awt) + Cigs_sub <- subset(Cigs, awt > 0) + iv_sub <- update(iv_full, data = Cigs_sub) + + CR_full <- lapply(CR_types, function(x) vcovCR(iv_full, cluster = Cigs$state, type = x)) + CR_sub <- lapply(CR_types, function(x) vcovCR(iv_sub, cluster = Cigs_sub$state, type = x)) + expect_equal(CR_full, CR_sub, check.attributes = FALSE) + + test_full <- lapply(CR_types, function(x) coef_test(iv_full, vcov = x, cluster = Cigs$state, test = c("z","naive-t","Satterthwaite"), p_values = TRUE)) + test_sub <- lapply(CR_types, function(x) coef_test(iv_sub, vcov = x, cluster = Cigs_sub$state, test = c("z","naive-t","Satterthwaite"), p_values = TRUE)) + expect_equal(test_full, test_sub, check.attributes = FALSE) + + dat_miss <- Cigs + miss_indicator <- sample.int(n_Cigs, size = round(n_Cigs / 10)) + dat_miss$rincome[miss_indicator] <- NA + with(dat_miss, table(awt, is.na(rincome))) + + iv_dropped <- update(iv_full, data = dat_miss) + dat_complete <- subset(dat_miss, !is.na(rincome)) + iv_complete <- update(iv_full, data = dat_complete) + + CR_drop <- lapply(CR_types, function(x) vcovCR(iv_dropped, cluster = dat_miss$state, type = x)) + CR_complete <- lapply(CR_types, function(x) vcovCR(iv_complete, cluster = dat_complete$state, type = x)) + expect_equal(CR_drop, CR_complete) + + test_drop <- lapply(CR_types, function(x) coef_test(iv_dropped, vcov = x, cluster = dat_miss$state, test = "All", p_values = FALSE)) + test_complete <- lapply(CR_types, function(x) coef_test(iv_complete, vcov = x, cluster = dat_complete$state, test = "All", p_values = FALSE)) + expect_equal(test_drop, test_complete) + +}) + diff -Nru r-cran-clubsandwich-0.5.7/tests/testthat/test_ivreg_ivreg.R r-cran-clubsandwich-0.5.8/tests/testthat/test_ivreg_ivreg.R --- r-cran-clubsandwich-0.5.7/tests/testthat/test_ivreg_ivreg.R 1970-01-01 00:00:00.000000000 +0000 +++ r-cran-clubsandwich-0.5.8/tests/testthat/test_ivreg_ivreg.R 2022-08-12 19:46:50.000000000 +0000 @@ -0,0 +1,269 @@ +############################# + +context("ivreg::ivreg objects") +set.seed(20190513) + +skip_if_not_installed("ivreg") + +library(ivreg, quietly=TRUE) + +data("CigarettesSW", package = "AER") + +Cigs <- within(CigarettesSW, { + rprice <- price/cpi + rincome <- income/population/cpi + tdiff <- (taxs - tax)/cpi +}) + +CR_types <- paste0("CR",0:4) + +obj_un <- ivreg::ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), + data = Cigs) +obj_wt <- ivreg::ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), + data = Cigs, + weights = population) + +X <- model.matrix(obj_wt, component = "regressors") +Z <- model.matrix(obj_wt, component = "instruments") +y <- log(CigarettesSW$packs) +w <- weights(obj_wt) + +test_that("Basic calculations from ivreg agree for unweighted model.", { + XZ <- model.matrix(obj_un, component = "projected") + ZtZ_inv <- chol2inv(chol(t(Z) %*% Z)) + XZ_check <- Z %*% ZtZ_inv %*% t(Z) %*% X + + expect_equal(XZ, XZ_check, check.attributes=FALSE) + expect_equal(coef(obj_un), lm.fit(XZ, y)$coefficients) + expect_equal (bread(obj_un), chol2inv(chol(t(XZ) %*% XZ)) * nobs(obj_un), check.attributes=FALSE) + + hii <- diag(XZ %*% chol2inv(chol(t(XZ) %*% XZ)) %*% t(XZ)) + expect_equal(hatvalues(obj_un, type = "stage2"), hii) + + r <- as.vector(y - X %*% coef(obj_un)) + expect_equal(r, as.vector(residuals_CS(obj_un))) +}) + +test_that("Basic calculations from ivreg agree for weighted model.", { + XZ <- model.matrix(obj_wt, component = "projected") + ZwZ_inv <- chol2inv(chol(t(Z) %*% (w * Z))) + XZ_check <- Z %*% ZwZ_inv %*% t(Z) %*% (w * X) + + expect_equal(XZ, XZ_check, check.attributes=FALSE) + expect_equal(coef(obj_wt), lm.wfit(XZ, y, w)$coefficients) + expect_equal(bread(obj_wt), chol2inv(chol(t(XZ) %*% (w * XZ))) * nobs(obj_wt), check.attributes=FALSE) + + hii <- diag(X%*% chol2inv(chol(t(XZ) %*% (w * XZ))) %*% t(w * XZ)) + expect_false(all(hatvalues(obj_wt) == hii)) # does not agree because hatvalues doesn't work with weighting + + r <- as.vector(y - X %*% coef(obj_wt)) + expect_equal(r, as.vector(residuals_CS(obj_wt))) +}) + +test_that("bread works", { + + expect_true(check_bread(obj_un, cluster = Cigs$state, y = log(Cigs$packs))) + tsls_vcov <- bread(obj_un) * summary(obj_un)$sigma^2 / v_scale(obj_un) + expect_equal(vcov(obj_un), tsls_vcov) + + expect_true(check_bread(obj_wt, cluster = Cigs$state, y = log(Cigs$packs))) + wtsls_vcov <- bread(obj_wt) * summary(obj_wt)$sigma^2 / v_scale(obj_wt) + expect_equal(vcov(obj_wt), wtsls_vcov) +}) + + +test_that("vcovCR options don't matter for CR0", { + expect_error(vcovCR(obj_un, type = "CR0")) + CR0 <- vcovCR(obj_un, cluster = Cigs$state, type = "CR0") + expect_output(print(CR0)) + attr(CR0, "target") <- NULL + attr(CR0, "inverse_var") <- NULL + CR0_A <- vcovCR(obj_un, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population) + attr(CR0_A, "target") <- NULL + attr(CR0_A, "inverse_var") <- NULL + expect_identical(CR0_A, CR0) + CR0_B <- vcovCR(obj_un, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population, inverse_var = FALSE) + attr(CR0_B, "target") <- NULL + attr(CR0_B, "inverse_var") <- NULL + expect_identical(CR0_A, CR0) + + expect_error(vcovCR(obj_un, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population, inverse_var = TRUE)) + + wCR0 <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR0") + attr(wCR0, "target") <- NULL + attr(wCR0, "inverse_var") <- NULL + wCR0_A <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population) + attr(wCR0_A, "target") <- NULL + attr(wCR0_A, "inverse_var") <- NULL + expect_identical(wCR0_A, wCR0) + wCR0_B <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population, inverse_var = FALSE) + attr(wCR0_B, "target") <- NULL + attr(wCR0_B, "inverse_var") <- NULL + expect_identical(wCR0_B, wCR0) + + expect_error(vcovCR(obj_wt, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population, inverse_var = TRUE)) +}) + +test_that("vcovCR options work for CR2", { + + CR2_iv <- vcovCR(obj_un, cluster = Cigs$state, type = "CR2") + expect_equal(vcovCR(obj_un, cluster = Cigs$state, type = "CR2", target = rep(1, nobs(obj_un))), CR2_iv) + expect_false(identical(vcovCR(obj_un, cluster = Cigs$state, type = "CR2", target = 1 / Cigs$population), CR2_iv)) + + wCR2_id <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR2") + expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR2", inverse_var = FALSE), wCR2_id) + expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR2", target = rep(1, nobs(obj_un))), wCR2_id) + expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR2", target = rep(1, nobs(obj_un)), inverse_var = FALSE), wCR2_id) + +}) + +test_that("vcovCR options work for CR4", { + + CR4_not <- vcovCR(obj_un, cluster = Cigs$state, type = "CR4") + expect_identical(vcovCR(obj_un, cluster = Cigs$state, type = "CR4", target = rep(1, nobs(obj_un))), CR4_not) + expect_identical(vcovCR(obj_un, cluster = Cigs$state, type = "CR4", target = rep(1, nobs(obj_un)), inverse_var = FALSE), CR4_not) + expect_false(identical(vcovCR(obj_un, cluster = Cigs$state, type = "CR4", target = 1 / Cigs$population), CR4_not)) + + wCR4_id <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR4") + expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR4", inverse_var = FALSE), wCR4_id) + expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR4", target = rep(1, nobs(obj_wt))), wCR4_id) + expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR4", target = rep(1, nobs(obj_wt)), inverse_var = FALSE), wCR4_id) + +}) + + +test_that("CR2 is target-unbiased", { + expect_true(check_CR(obj_un, vcov = "CR2", cluster = Cigs$state)) + expect_true(check_CR(obj_wt, vcov = "CR2", cluster = Cigs$state)) +}) + +test_that("CR4 is target-unbiased", { + skip("Need to understand target-unbiasedness for ivreg objects.") + expect_true(check_CR(obj_un, vcov = "CR4", cluster = Cigs$state)) + expect_true(check_CR(obj_wt, vcov = "CR4", cluster = Cigs$state)) +}) + +test_that("vcovCR is equivalent to vcovHC (with HC0 or HC1) when clusters are all of size 1", { + library(sandwich, quietly=TRUE) + CR0 <- vcovCR(obj_un, cluster = 1:nobs(obj_un), type = "CR0") + expect_equal(vcovHC(obj_un, type = "HC0"), as.matrix(CR0)) + CR1 <- vcovCR(obj_un, cluster = 1:nobs(obj_un), type = "CR1S") + expect_equal(vcovHC(obj_un, type = "HC1"), as.matrix(CR1)) + CR2 <- vcovCR(obj_un, cluster = 1:nobs(obj_un), type = "CR2") + expect_false(all(vcovHC(obj_un, type = "HC2") == as.matrix(CR2))) + CR3 <- vcovCR(obj_un, cluster = 1:nobs(obj_un), type = "CR3") + expect_false(all(vcovHC(obj_un, type = "HC3") == as.matrix(CR3))) +}) + +test_that("Order doesn't matter.",{ + + check_sort_order(obj_wt, Cigs, "state") + +}) +# Clustering variable must have length equal to the number of rows in the data used to fit obj. + +test_that("clubSandwich works with dropped observations", { + dat_miss <- Cigs + dat_miss$rincome[sample.int(nrow(Cigs), size = round(nrow(Cigs) / 10))] <- NA + iv_dropped <- update(obj_un, data = dat_miss) + dat_complete <- subset(dat_miss, !is.na(rincome)) + iv_complete <- update(obj_un, data = dat_complete) + + CR_drop <- lapply(CR_types, function(x) vcovCR(iv_dropped, cluster = dat_miss$state, type = x)) + CR_complete <- lapply(CR_types, function(x) vcovCR(iv_complete, cluster = dat_complete$state, type = x)) + expect_equal(CR_drop, CR_complete) + + test_drop <- lapply(CR_types, function(x) coef_test(iv_dropped, vcov = x, cluster = dat_miss$state, test = "All", p_values = FALSE)) + test_complete <- lapply(CR_types, function(x) coef_test(iv_complete, vcov = x, cluster = dat_complete$state, test = "All", p_values = FALSE)) + expect_equal(test_drop, test_complete) +}) + + + +test_that("weight scale doesn't matter", { + + iv_fit_w <- update(obj_un, weights = rep(4, nobs(obj_un))) + + unweighted_fit <- lapply(CR_types, function(x) vcovCR(obj_un, cluster = Cigs$state, type = x)) + weighted_fit <- lapply(CR_types, function(x) vcovCR(iv_fit_w, cluster = Cigs$state, type = x)) + + expect_equal(lapply(unweighted_fit, as.matrix), + lapply(weighted_fit, as.matrix), + tol = 5 * 10^-7) + + target <- 1 + rpois(nrow(Cigs), lambda = 8) + unweighted_fit <- lapply(CR_types, function(x) vcovCR(obj_un, cluster = Cigs$state, type = x, target = target)) + weighted_fit <- lapply(CR_types, function(x) vcovCR(iv_fit_w, cluster = Cigs$state, type = x, target = target * 15)) + + expect_equal(lapply(unweighted_fit, as.matrix), + lapply(weighted_fit, as.matrix), + tol = 5 * 10^-7) + +}) + +test_that("clubSandwich works with weights of zero.", { + + n_Cigs <- nrow(Cigs) + Cigs$awt <- rpois(n_Cigs, lambda = 1.4) + table(Cigs$awt) + + iv_full <- update(obj_un, weights = awt) + Cigs_sub <- subset(Cigs, awt > 0) + iv_sub <- update(iv_full, data = Cigs_sub) + + CR_full <- lapply(CR_types, function(x) vcovCR(iv_full, cluster = Cigs$state, type = x)) + CR_sub <- lapply(CR_types, function(x) vcovCR(iv_sub, cluster = Cigs_sub$state, type = x)) + expect_equal(CR_full, CR_sub, check.attributes = FALSE) + + test_full <- lapply(CR_types, function(x) coef_test(iv_full, vcov = x, cluster = Cigs$state, test = c("z","naive-t","Satterthwaite"), p_values = TRUE)) + test_sub <- lapply(CR_types, function(x) coef_test(iv_sub, vcov = x, cluster = Cigs_sub$state, test = c("z","naive-t","Satterthwaite"), p_values = TRUE)) + expect_equal(test_full, test_sub, check.attributes = FALSE) + + dat_miss <- Cigs + miss_indicator <- sample.int(n_Cigs, size = round(n_Cigs / 10)) + dat_miss$rincome[miss_indicator] <- NA + with(dat_miss, table(awt, is.na(rincome))) + + iv_dropped <- update(iv_full, data = dat_miss) + dat_complete <- subset(dat_miss, !is.na(rincome)) + iv_complete <- update(iv_full, data = dat_complete) + + CR_drop <- lapply(CR_types, function(x) vcovCR(iv_dropped, cluster = dat_miss$state, type = x)) + CR_complete <- lapply(CR_types, function(x) vcovCR(iv_complete, cluster = dat_complete$state, type = x)) + expect_equal(CR_drop, CR_complete) + + test_drop <- lapply(CR_types, function(x) coef_test(iv_dropped, vcov = x, cluster = dat_miss$state, test = "All", p_values = FALSE)) + test_complete <- lapply(CR_types, function(x) coef_test(iv_complete, vcov = x, cluster = dat_complete$state, test = "All", p_values = FALSE)) + expect_equal(test_drop, test_complete) + +}) + + +#------------------------------------------------------------------------------- +# Other estimation methods + +ols_un <- ivreg::ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), + data = Cigs, method = "OLS") +ols_wt <- ivreg::ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), + data = Cigs, + weights = population, method = "OLS") + +mom_un <- ivreg::ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), + data = Cigs, method = "M") +mom_wt <- ivreg::ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), + data = Cigs, + weights = population, method = "M") + +rob_un <- ivreg::ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), + data = Cigs, method = "MM") +rob_wt <- ivreg::ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), + data = Cigs, + weights = population, method = "MM") + +test_that("vcovCR does not currently support ivreg models estimated using method = 'M' or method = 'MM'", { + expect_error(vcovCR(mom_un, cluster = Cigs$state, type = "CR2")) + expect_error(vcovCR(mom_wt, cluster = Cigs$state, type = "CR2")) + expect_error(vcovCR(rob_un, cluster = Cigs$state, type = "CR2")) + expect_error(vcovCR(rob_wt, cluster = Cigs$state, type = "CR2")) +}) + diff -Nru r-cran-clubsandwich-0.5.7/tests/testthat/test_ivreg.R r-cran-clubsandwich-0.5.8/tests/testthat/test_ivreg.R --- r-cran-clubsandwich-0.5.7/tests/testthat/test_ivreg.R 2022-06-11 01:49:36.000000000 +0000 +++ r-cran-clubsandwich-0.5.8/tests/testthat/test_ivreg.R 1970-01-01 00:00:00.000000000 +0000 @@ -1,239 +0,0 @@ -context("ivreg objects") -set.seed(20190513) - -skip_if_not_installed("zoo") -skip_if_not_installed("AER") - -library(zoo, quietly=TRUE) -library(AER, quietly=TRUE) - -data("CigarettesSW", package = "AER") - -Cigs <- within(CigarettesSW, { - rprice <- price/cpi - rincome <- income/population/cpi - tdiff <- (taxs - tax)/cpi -}) - -CR_types <- paste0("CR",0:4) - -obj_un <- ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), - data = Cigs) -obj_wt <- ivreg(log(packs) ~ log(rprice) + log(rincome) | log(rincome) + tdiff + I(tax/cpi), - data = Cigs, - weights = population) - -X <- model.matrix(obj_wt, component = "regressors") -Z <- model.matrix(obj_wt, component = "instruments") -y <- log(CigarettesSW$packs) -w <- weights(obj_wt) - -test_that("Basic calculations from ivreg agree for unweighted model.", { - XZ <- model.matrix(obj_un, component = "projected") - ZtZ_inv <- chol2inv(chol(t(Z) %*% Z)) - XZ_check <- Z %*% ZtZ_inv %*% t(Z) %*% X - - expect_equal(XZ, XZ_check, check.attributes=FALSE) - expect_equal(coef(obj_un), lm.fit(XZ, y)$coefficients) - expect_equal (bread(obj_un), chol2inv(chol(t(XZ) %*% XZ)) * nobs(obj_un), check.attributes=FALSE) - - hii <- diag(X %*% chol2inv(chol(t(XZ) %*% XZ)) %*% t(XZ)) - expect_equal(hatvalues(obj_un), hii) - - r <- as.vector(y - X %*% coef(obj_un)) - expect_equal(r, as.vector(residuals_CS(obj_un))) -}) - -test_that("Basic calculations from ivreg agree for weighted model.", { - XZ <- model.matrix(obj_wt, component = "projected") - ZwZ_inv <- chol2inv(chol(t(Z) %*% (w * Z))) - XZ_check <- Z %*% ZwZ_inv %*% t(Z) %*% (w * X) - - expect_equal(XZ, XZ_check, check.attributes=FALSE) - expect_equal(coef(obj_wt), lm.wfit(XZ, y, w)$coefficients) - expect_equal(bread(obj_wt), chol2inv(chol(t(XZ) %*% (w * XZ))) * nobs(obj_wt), check.attributes=FALSE) - - hii <- diag(X%*% chol2inv(chol(t(XZ) %*% (w * XZ))) %*% t(w * XZ)) - expect_false(all(hatvalues(obj_wt) == hii)) # does not agree because hatvalues doesn't work with weighting - - r <- as.vector(y - X %*% coef(obj_wt)) - expect_equal(r, as.vector(residuals_CS(obj_wt))) -}) - -test_that("bread works", { - - expect_true(check_bread(obj_un, cluster = Cigs$state, y = log(Cigs$packs))) - tsls_vcov <- bread(obj_un) * summary(obj_un)$sigma^2 / v_scale(obj_un) - expect_equal(vcov(obj_un), tsls_vcov) - - expect_true(check_bread(obj_wt, cluster = Cigs$state, y = log(Cigs$packs))) - wtsls_vcov <- bread(obj_wt) * summary(obj_wt)$sigma^2 / v_scale(obj_wt) - expect_equal(vcov(obj_wt), wtsls_vcov) -}) - - -test_that("vcovCR options don't matter for CR0", { - expect_error(vcovCR(obj_un, type = "CR0")) - CR0 <- vcovCR(obj_un, cluster = Cigs$state, type = "CR0") - expect_output(print(CR0)) - attr(CR0, "target") <- NULL - attr(CR0, "inverse_var") <- NULL - CR0_A <- vcovCR(obj_un, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population) - attr(CR0_A, "target") <- NULL - attr(CR0_A, "inverse_var") <- NULL - expect_identical(CR0_A, CR0) - CR0_B <- vcovCR(obj_un, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population, inverse_var = FALSE) - attr(CR0_B, "target") <- NULL - attr(CR0_B, "inverse_var") <- NULL - expect_identical(CR0_A, CR0) - - expect_error(vcovCR(obj_un, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population, inverse_var = TRUE)) - - wCR0 <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR0") - attr(wCR0, "target") <- NULL - attr(wCR0, "inverse_var") <- NULL - wCR0_A <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population) - attr(wCR0_A, "target") <- NULL - attr(wCR0_A, "inverse_var") <- NULL - expect_identical(wCR0_A, wCR0) - wCR0_B <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population, inverse_var = FALSE) - attr(wCR0_B, "target") <- NULL - attr(wCR0_B, "inverse_var") <- NULL - expect_identical(wCR0_B, wCR0) - - expect_error(vcovCR(obj_wt, cluster = Cigs$state, type = "CR0", target = 1 / Cigs$population, inverse_var = TRUE)) -}) - -test_that("vcovCR options work for CR2", { - - CR2_iv <- vcovCR(obj_un, cluster = Cigs$state, type = "CR2") - expect_equal(vcovCR(obj_un, cluster = Cigs$state, type = "CR2", target = rep(1, nobs(obj_un))), CR2_iv) - expect_false(identical(vcovCR(obj_un, cluster = Cigs$state, type = "CR2", target = 1 / Cigs$population), CR2_iv)) - - wCR2_id <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR2") - expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR2", inverse_var = FALSE), wCR2_id) - expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR2", target = rep(1, nobs(obj_un))), wCR2_id) - expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR2", target = rep(1, nobs(obj_un)), inverse_var = FALSE), wCR2_id) - -}) - -test_that("vcovCR options work for CR4", { - - CR4_not <- vcovCR(obj_un, cluster = Cigs$state, type = "CR4") - expect_identical(vcovCR(obj_un, cluster = Cigs$state, type = "CR4", target = rep(1, nobs(obj_un))), CR4_not) - expect_identical(vcovCR(obj_un, cluster = Cigs$state, type = "CR4", target = rep(1, nobs(obj_un)), inverse_var = FALSE), CR4_not) - expect_false(identical(vcovCR(obj_un, cluster = Cigs$state, type = "CR4", target = 1 / Cigs$population), CR4_not)) - - wCR4_id <- vcovCR(obj_wt, cluster = Cigs$state, type = "CR4") - expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR4", inverse_var = FALSE), wCR4_id) - expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR4", target = rep(1, nobs(obj_wt))), wCR4_id) - expect_identical(vcovCR(obj_wt, cluster = Cigs$state, type = "CR4", target = rep(1, nobs(obj_wt)), inverse_var = FALSE), wCR4_id) - -}) - - -test_that("CR2 is target-unbiased", { - expect_true(check_CR(obj_un, vcov = "CR2", cluster = Cigs$state)) - expect_true(check_CR(obj_wt, vcov = "CR2", cluster = Cigs$state)) -}) - -test_that("CR4 is target-unbiased", { - skip("Need to understand target-unbiasedness for ivreg objects.") - expect_true(check_CR(obj_un, vcov = "CR4", cluster = Cigs$state)) - expect_true(check_CR(obj_wt, vcov = "CR4", cluster = Cigs$state)) -}) - -test_that("vcovCR is equivalent to vcovHC (with HC0 or HC1) when clusters are all of size 1", { - library(sandwich, quietly=TRUE) - CR0 <- vcovCR(obj_un, cluster = 1:nobs(obj_un), type = "CR0") - expect_equal(vcovHC(obj_un, type = "HC0"), as.matrix(CR0)) - CR1 <- vcovCR(obj_un, cluster = 1:nobs(obj_un), type = "CR1S") - expect_equal(vcovHC(obj_un, type = "HC1"), as.matrix(CR1)) - CR2 <- vcovCR(obj_un, cluster = 1:nobs(obj_un), type = "CR2") - expect_false(all(vcovHC(obj_un, type = "HC2") == as.matrix(CR2))) - CR3 <- vcovCR(obj_un, cluster = 1:nobs(obj_un), type = "CR3") - expect_false(all(vcovHC(obj_un, type = "HC3") == as.matrix(CR3))) -}) - -test_that("Order doesn't matter.",{ - - check_sort_order(obj_wt, Cigs, "state") - -}) - -test_that("clubSandwich works with dropped observations", { - dat_miss <- Cigs - dat_miss$rincome[sample.int(nrow(Cigs), size = round(nrow(Cigs) / 10))] <- NA - iv_dropped <- update(obj_un, data = dat_miss) - dat_complete <- subset(dat_miss, !is.na(rincome)) - iv_complete <- update(obj_un, data = dat_complete) - - CR_drop <- lapply(CR_types, function(x) vcovCR(iv_dropped, cluster = dat_miss$state, type = x)) - CR_complete <- lapply(CR_types, function(x) vcovCR(iv_complete, cluster = dat_complete$state, type = x)) - expect_equal(CR_drop, CR_complete) - - test_drop <- lapply(CR_types, function(x) coef_test(iv_dropped, vcov = x, cluster = dat_miss$state, test = "All", p_values = FALSE)) - test_complete <- lapply(CR_types, function(x) coef_test(iv_complete, vcov = x, cluster = dat_complete$state, test = "All", p_values = FALSE)) - expect_equal(test_drop, test_complete) -}) - - - -test_that("weight scale doesn't matter", { - - iv_fit_w <- update(obj_un, weights = rep(4, nobs(obj_un))) - - unweighted_fit <- lapply(CR_types, function(x) vcovCR(obj_un, cluster = Cigs$state, type = x)) - weighted_fit <- lapply(CR_types, function(x) vcovCR(iv_fit_w, cluster = Cigs$state, type = x)) - - expect_equal(lapply(unweighted_fit, as.matrix), - lapply(weighted_fit, as.matrix), - tol = 5 * 10^-7) - - target <- 1 + rpois(nrow(Cigs), lambda = 8) - unweighted_fit <- lapply(CR_types, function(x) vcovCR(obj_un, cluster = Cigs$state, type = x, target = target)) - weighted_fit <- lapply(CR_types, function(x) vcovCR(iv_fit_w, cluster = Cigs$state, type = x, target = target * 15)) - - expect_equal(lapply(unweighted_fit, as.matrix), - lapply(weighted_fit, as.matrix), - tol = 5 * 10^-7) - -}) - -test_that("clubSandwich works with weights of zero.", { - - n_Cigs <- nrow(Cigs) - Cigs$awt <- rpois(n_Cigs, lambda = 1.4) - table(Cigs$awt) - - iv_full <- update(obj_un, weights = awt) - Cigs_sub <- subset(Cigs, awt > 0) - iv_sub <- update(iv_full, data = Cigs_sub) - - CR_full <- lapply(CR_types, function(x) vcovCR(iv_full, cluster = Cigs$state, type = x)) - CR_sub <- lapply(CR_types, function(x) vcovCR(iv_sub, cluster = Cigs_sub$state, type = x)) - expect_equal(CR_full, CR_sub, check.attributes = FALSE) - - test_full <- lapply(CR_types, function(x) coef_test(iv_full, vcov = x, cluster = Cigs$state, test = c("z","naive-t","Satterthwaite"), p_values = TRUE)) - test_sub <- lapply(CR_types, function(x) coef_test(iv_sub, vcov = x, cluster = Cigs_sub$state, test = c("z","naive-t","Satterthwaite"), p_values = TRUE)) - expect_equal(test_full, test_sub, check.attributes = FALSE) - - dat_miss <- Cigs - miss_indicator <- sample.int(n_Cigs, size = round(n_Cigs / 10)) - dat_miss$rincome[miss_indicator] <- NA - with(dat_miss, table(awt, is.na(rincome))) - - iv_dropped <- update(iv_full, data = dat_miss) - dat_complete <- subset(dat_miss, !is.na(rincome)) - iv_complete <- update(iv_full, data = dat_complete) - - CR_drop <- lapply(CR_types, function(x) vcovCR(iv_dropped, cluster = dat_miss$state, type = x)) - CR_complete <- lapply(CR_types, function(x) vcovCR(iv_complete, cluster = dat_complete$state, type = x)) - expect_equal(CR_drop, CR_complete) - - test_drop <- lapply(CR_types, function(x) coef_test(iv_dropped, vcov = x, cluster = dat_miss$state, test = "All", p_values = FALSE)) - test_complete <- lapply(CR_types, function(x) coef_test(iv_complete, vcov = x, cluster = dat_complete$state, test = "All", p_values = FALSE)) - expect_equal(test_drop, test_complete) - -}) - diff -Nru r-cran-clubsandwich-0.5.7/tests/testthat/test_plm-first-differences.R r-cran-clubsandwich-0.5.8/tests/testthat/test_plm-first-differences.R --- r-cran-clubsandwich-0.5.7/tests/testthat/test_plm-first-differences.R 2022-06-11 01:49:36.000000000 +0000 +++ r-cran-clubsandwich-0.5.8/tests/testthat/test_plm-first-differences.R 2022-08-12 20:43:15.000000000 +0000 @@ -34,6 +34,11 @@ as.matrix(vcovCR(plm_FD, type = "CR0")), check.attributes = FALSE) expect_equal(vcovHC(plm_FD, method="arellano", type = "sss", cluster = "group"), as.matrix(vcovCR(plm_FD, type = "CR1S")), check.attributes = FALSE) +}) + +test_that("CR0 and CR1S agree with arellano vcov for versions <= 2.6-1", { + + skip_if(packageVersion("plm") > "2.6-1") X <- model_matrix(plm_FD) e <- residuals(plm_FD) @@ -55,6 +60,16 @@ vcovHC(plm_FD, method="arellano", type = "HC0", cluster = "time"), check.attributes = FALSE) }) +test_that("CR0 and CR1S agree with arellano vcov for versions > 2.6-1", { + + skip_if_not_installed("plm", minimum_version = "2.6-2") + + expect_equal(vcovHC(plm_FD, method="arellano", type = "HC0", cluster = "time"), + as.matrix(vcovCR(plm_FD, type = "CR0", cluster = "time")), check.attributes = FALSE) + expect_equal(vcovHC(plm_FD, method="arellano", type = "sss", cluster = "time"), + as.matrix(vcovCR(plm_FD, type = "CR1S", cluster = "time")), check.attributes = FALSE) +}) + test_that("vcovCR options work for CR2", { CR2_iv <- vcovCR(plm_FD, type = "CR2")