tests/testthat/test_equal_to_sandwich.R

In baruuum/jars: Just Another Robust Standard errors package

context("Check that results agree with sandwich package")

# generate data
n = 300
x1 = runif(n)
x2 = rnorm(n)
z = sample(1:6, n, replace = T)
y = 2 + runif(1, -1, 1) * x1 + runif(1, -1, 1) * x2 + rnorm(n)

# fit model
m = lm(y ~ x1 + x2)

for (a in paste0("HC", 0:4)) {
    
    test_that(paste0("Het: estimators agree on ", a), {
        
            expect_equal(
                jars::.robustSE(model.matrix(m), resid(m), type = a),
                sandwich::vcovHC(m, type = a),
                check.attributes = F
            )
        
            expect_equal(
                vcov(summary(m, robust = TRUE, type = a)),
                sandwich::vcovHC(m, type = a),
                check.attributes = F
            )
        }
    )
}


for (a in paste0("HC", 0:1)) {
    
    test_that(paste0("Clust: estimators agree on ", a), {
        
            expect_equal(
                jars::.clusterSE(model.matrix(m), resid(m), clust = z, type = a),
                sandwich::vcovCL(m, cluster = z, type = a),
                check.attributes = F
            )
        
            expect_equal(
                vcov(summary(m, cluster = z, type = a)),
                sandwich::vcovCL(m, cluster = z, type = a),
                check.attributes = F
            )
        }
    )
}



context("Behavior When Missing Values are Present")

# generate some missing values
dat = data.frame(y, x1, x2)
n_miss = 20L
which_var = sample(1:3, n_miss, replace = T)
which_row = sample.int(n, n_miss, replace = F)
for (i in 1:20) dat[which_row[i], which_var[i]] = NA

# refit model
m = lm(y ~ x1 + x2, data = dat)

for (a in paste0("HC", 0:4)) {
    
    test_that(paste0("Het: estimators agree on ", 
                     a, " in the presence of NA"), {
        
            expect_equal(
                vcov(summary(m, robust = TRUE, type = a)),
                sandwich::vcovHC(m, type = a),
                check.attributes = F
            )

        }
    )
}

for (a in paste0("HC", 0:1)) {
    
    test_that(paste0("Clust: estimators agree on ", 
                      a, " in the presence of NA (in vars)"), {
        
            expect_equal(
                vcov(summary(m, cluster = z, type = a)),
                sandwich::vcovCL(m, cluster = z, type = a),
                check.attributes = F
            )

        }
    )
}



# add missing values to cluster var
z[sample.int(n, 5, replace = F)] = NA

for (a in paste0("HC", 0:1)) {
    
    test_that(paste0("Clust: estimators agree on ", 
                      a, " in the presence of NA (in vars & clust)"), {
        
            expect_false(
                isTRUE(all.equal(
                    vcov(summary(m, cluster = z, type = a)),
                    sandwich::vcovCL(m, cluster = z, type = a),
                    check.attributes = F
                ))
            )
                          
            expect_warning(
                vcov(summary(m, cluster = z, type = a))
            )

        }
    )
}

# recode missing values in clusters as own cluster
z[is.na(z)] = 99

for (a in paste0("HC", 0:1)) {
    
    test_that(
        paste0("Clust: estimators agree on ", 
               a, " in the presence of NA (in vars & clust)"), {
        
            expect_equal(
                    vcov(summary(m, cluster = z, type = a)),
                    sandwich::vcovCL(m, cluster = z, type = a),
                    check.attributes = F
            )

        }
    )
}