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 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
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(mrate ~ legal + beertaxa, data = MV_deaths,
- plm_unweighted effect = "twoways", index = c("state","year"))
- coef_test(plm_unweighted, vcov = "CR1", cluster = "individual", test = "naive-t")
library(plm)
+<- plm(mrate ~ legal + beertaxa, data = MV_deaths,
+ plm_unweighted 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(mrate ~ 0 + legal + beertaxa + factor(state) + factor(year),
- lm_weighted weights = pop, data = MV_deaths)
- coef_test(lm_weighted, vcov = "CR1",
-cluster = MV_deaths$state, test = "naive-t")[1:2,]
<- lm(mrate ~ 0 + legal + beertaxa + factor(state) + factor(year),
+ lm_weighted 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(mrate ~ 0 + legal + beertaxa + year, data = MV_deaths,
- plm_random effect = "individual", index = c("state","year"),
- model = "random")
- coef_test(plm_random, vcov = "CR1", test = "naive-t")[1:2,]
<- plm(mrate ~ 0 + legal + beertaxa + year, data = MV_deaths,
+ plm_random 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):
-<- within(MV_deaths, {
- MV_deaths <- legal - tapply(legal, state, mean)[factor(state)]
- legal_cent <- beertaxa - tapply(beertaxa, state, mean)[factor(state)]
- beer_cent
- })
-<- plm(mrate ~ 0 + legal + beertaxa + legal_cent + beer_cent + factor(year),
- plm_Hausman data = MV_deaths,
- effect = "individual", index = c("state","year"),
- model = "random")
- coef_test(plm_Hausman, vcov = "CR2", test = "Satterthwaite")[1:4,]
<- within(MV_deaths, {
+ MV_deaths <- legal - tapply(legal, state, mean)[factor(state)]
+ legal_cent <- beertaxa - tapply(beertaxa, state, mean)[factor(state)]
+ beer_cent
+ })
+<- plm(mrate ~ 0 + legal + beertaxa + legal_cent + beer_cent + factor(year),
+ plm_Hausman 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")