tests/testthat/test_predict_ci.R

In baruuum/btoolbox: random toolbox

context("Test errors in predict_ci")

# obs
n = 30000

# coefficients
true_beta = c(.8, .25)
true_tau = c(-.5, 2)

# predictors (dummy & continuous)
x1 = sample(c(0, 1), n, replace = T)
x2 = runif(n, -2, 2)

# linear predictor
xb = cbind(x1, x2) %*% true_beta

# generate latent outcome
ystar = as.numeric(xb + rlogis(n))

# generate obseved outcome
y = ifelse(
    ystar < true_tau[1], 1,
        ifelse(ystar >= true_tau[1] & ystar < true_tau[2],
               2, 3)
)

N = rpois(1L, 15)
pred_dat = data.frame(
    x1 = sample.int(2, N, replace = T) - 1L, 
    x2 = rnorm(N)
)

test_that("function works for lm/glm models", {
    
    y_bin = ifelse(y > 1, 1, 0)
    
    f_list = list(binomial("logit"),
                  binomial("probit"),
                  quasibinomial(link = "logit"),
                  Gamma("inverse"),
                  Gamma("log"),
                  poisson("log"),
                  quasipoisson("log"))
    
    for (i in 1:(length(f_list) + 1L)) {
        
        if (i <= 3) {
        
            f = glm(y_bin ~ x1 + x2, family = f_list[[i]])
            
        } else if (i > length(f_list))  {
            
            f = lm(y ~ x1 + x2)
            
        } else {
            
            f = glm(y ~ x1 + x2, family = f_list[[i]])
            
        }
        
        pred =  predict_ci(f, pred_dat, level = .5)
        
        expect_equal(nrow(pred), nrow(pred_dat))
        expect_equal(ncol(pred), 2L)
        
        pred_point = predict(f, pred_dat, type = "response")
        
        expect_true(
            all(pred_point >= pred[,1] & pred_point <= pred[,2])
        )
        
        # get Sigma
        Sigma = vcov(f)
        
        # run ci with sigma
        pred_sig =  predict_ci(f, pred_dat, Sigma, level = .5)
        
        expect_equal(pred, pred_sig)
        
        # get het-sigma
        Sigma_het = sandwich::vcovHC(f, type = "HC1")
        
        pred_het =  predict_ci(
            f, pred_dat, Sigma_het, level = .5)
        
        expect_true(!isTRUE(all.equal(pred, pred_het)))
        
    }

})


test_that("function works for polr/multinom models", {
    
    for (i in 1:2) {
        
        if (i == 1) {
        
            f = MASS::polr(factor(y) ~ x1 + x2, Hess = T)
            
        } else {
            
            f = nnet::multinom(factor(y) ~ x1 + x2, Hess = T)
            
        }
        
        pred =  predict_ci(f, pred_dat, level = .5, n_sim = 50000)
        
        expect_equal(dim(pred)[1], nrow(pred_dat))
        expect_equal(dim(pred)[2], 2L)
        expect_equal(dim(pred)[3], length(unique(y)))
        
        pred_point = predict(f, pred_dat, type = "p")
        
        for (j in seq.int(ncol(pred_point))) {
            expect_true(
                all(pred_point[, j] >= pred[, 1, j] & 
                        pred_point[, j] <= pred[, 2, j])
            )
        }
        
        # get Sigma
        Sigma = vcov(f)
        
        # run ci with sigma
        pred_sig =  predict_ci(f, pred_dat, Sigma, level = .5, n_sim = 50000)
        
        expect_true(all(abs(pred - pred_sig) < .001))
        
        # # get het-sigma (cannot implement robust)
        # Sigma_het = sandwich::vcovHC(f, type = "HC1")
        # 
        # pred_het =  predict_ci(
        #     f, pred_dat, Sigma_het, level = .5)
        # 
        # expect_true(!isTRUE(all.equal(pred, pred_het)))
        # 
    }

})