Nothing
tests/testthat/test-resi.R
In RESI: Robust Effect Size Index (RESI) Estimation
library(car)
#library(splines)
library(sandwich)
#library(survival)
library(nlme)
#library(lme4)
library(tibble)
#library(geepack)
data = RESI::insurance
mod = glm(charges ~ region * age + bmi + sex, data = data)
mod.r = glm(charges ~ bmi +sex, data = data)
mod.lm = lm(charges ~ region * age + bmi + sex, data = data)
mod.log = glm(smoker ~ age + region, data = data, family = "binomial")
mod.glm.int = glm(smoker ~ 1, data = data, family = "binomial")
mod.lm.int = lm(charges ~ 1, data = data)
# splines
if(requireNamespace("splines")){
mod.s = glm(charges ~ region * splines::ns(age, df=3) + bmi + sex, data = data)
mod.s.r = glm(charges ~ region * splines::ns(age, df=3), data = data)
mod.lm.s = lm(charges ~ region*splines::ns(age, df = 3) + bmi + sex, data = data)
mod.lm.s.r = lm(charges ~ region * splines::ns(age, df=3), data = data)
}
# missing data models
data.nas = data
data.nas[12:17, 'bmi'] = NA
mod.na = lm(charges ~ region * age + bmi + sex, data = data.nas)
mod.log.na = glm(smoker ~ age + region + bmi, data = data.nas, family = "binomial")
#nls
x <- seq(0, 100, 1)
y<-((runif(1, 10, 20)*x)/(runif(1, 0, 10) + x)) +
rnorm(101, 0, 1)
data.nls = data.frame(X = x, Y = y)
mod.nls = nls(Y~a * X/(b + X), data = data.nls, start = list(a = 1, b = 2))
#survival
if(requireNamespace("survival")){
data.surv = survival::lung
mod.surv = survival::survreg(survival::Surv(time, status) ~ age + sex + ph.karno, data=data.surv,
dist='weibull', robust = TRUE)
mod.surv.red = survival::survreg(survival::Surv(time, status) ~ age, data=data.surv[which(!(1:nrow(data.surv)%in% mod.surv$na.action)),],
dist='weibull', robust = TRUE)
mod.surv.nr = survival::survreg(survival::Surv(time, status) ~ age + sex + ph.karno, data=data.surv,
dist='weibull', robust = FALSE)
mod.surv.int = survival::survreg(survival::Surv(time, status) ~ 1, data=data.surv,
dist='weibull', robust = TRUE)
mod.coxph = survival::coxph(survival::Surv(time, status) ~ age + sex + wt.loss, data=survival::lung, robust = TRUE)
mod.coxph.red = survival::coxph(survival::Surv(time, status) ~ age, data=survival::lung, robust = TRUE)
}
# hurdle
if(requireNamespace("pscl")){
data.hurdle = pscl::bioChemists
mod.hurdle = pscl::hurdle(art ~ fem + mar + kid5 + phd + ment | fem + mar + kid5 + phd + ment, data = data.hurdle)
mod.hurdle.r = pscl::hurdle(art ~ kid5 + phd | kid5 + phd, data = data.hurdle)
mod.hurdle.int = pscl::hurdle(art ~ 1, data = data.hurdle)
# zeroinfl
mod.zinf = pscl::zeroinfl(art ~ fem + mar + kid5 + phd + ment | fem + mar + kid5 + phd + ment, data = data.hurdle)
mod.zinf.r = pscl::zeroinfl(art ~ kid5 + phd | kid5 + phd, data = data.hurdle)
mod.zinf.int = pscl::zeroinfl(art ~ 1, data = data.hurdle)
}
data.gee = RESI::depression
# gee and lme
if(requireNamespace("gee")){
mod.gee = gee::gee(depression ~ diagnose + drug*time, data = data.gee,
id = id, family = binomial, corstr = "independence")
mod.gee.int = gee::gee(depression ~ 1, data = data.gee,
id = id, family = binomial, corstr = "independence")
}
if(requireNamespace("geepack")){
library(geepack)
mod.geeglm = geepack::geeglm(depression ~ diagnose + drug*time, data = data.gee,
id = id, family = binomial, corstr = "independence")
mod.geeglm.r = geepack::geeglm(depression ~ diagnose, data = data.gee, id = id,
family = binomial, corstr = "independence")
mod.geeglm.int = geepack::geeglm(depression ~ 1, data = data.gee, id = id,
family = binomial, corstr = "independence")
}
if(requireNamespace("glmtoolbox")){
#library(glmtoolbox)
mod.glmgee = glmtoolbox::glmgee(depression ~ diagnose + drug*time, data = data.gee,
id = id, family = binomial, corstr = "independence")
mod.glmgee.r = glmtoolbox::glmgee(depression ~ diagnose, data = data.gee, id = id,
family = binomial, corstr = "independence")
mod.glmgee.int = glmtoolbox::glmgee(depression ~ 1, data = data.gee, id = id,
family = binomial, corstr = "independence")
}
mod.lme = nlme::lme(distance ~ age + Sex, data = nlme::Orthodont, random = ~ 1)
mod.lme.int = nlme::lme(distance ~ 1, data = nlme::Orthodont, random = ~ 1)
if(requireNamespace("lme4")){
mod.lmerMod = lme4::lmer(Reaction ~ Days + (Days | Subject), lme4::sleepstudy)
data.lmer = nlme::Orthodont
data.lmer$nsex <- as.numeric(data.lmer$Sex=="Male")
data.lmer$nsexage <- with(data.lmer, nsex*age)
mod.lmerMod2 = lme4::lmer(distance ~ age + (age|Subject) + (0+nsex|Subject) +
(0 + nsexage|Subject), data=data.lmer)
mod.lmerMod.int = lme4::lmer(Reaction ~ 1 + (1 | Subject), lme4::sleepstudy)
}
# tibbles
if(requireNamespace("tibble")){
tib.dat <- tibble::as_tibble(data)
mod.glm.tib <- glm(charges ~ region * age + bmi + sex, data = tib.dat)
mod.lm.tib <- lm(charges ~ region * age + bmi + sex, data = tib.dat)
tib.nls <- tibble::as_tibble(data.nls)
mod.nls.tib <- nls(Y~a * X/(b + X), data = tib.nls, start = list(a = 1, b = 2))
if(requireNamespace("survival")){
tib.surv <- tibble::as_tibble(data.surv)
mod.surv.tib <- survival::survreg(survival::Surv(time, status) ~ age + sex + ph.karno, data=tib.surv,
dist='weibull', robust = TRUE)
mod.coxph.tib <- survival::coxph(survival::Surv(time, status) ~ age + sex + wt.loss,
data=tibble::as_tibble(survival::lung), robust = TRUE)
}
if(requireNamespace("pscl")){
tib.hurdle <- tibble::as_tibble(data.hurdle)
mod.hurdle.tib <- pscl::hurdle(art ~ fem + mar + kid5 + phd + ment | fem + mar +
kid5 + phd + ment, data = tib.hurdle)
mod.zinf.tib <- pscl::zeroinfl(art ~ fem + mar + kid5 + phd + ment | fem + mar +
kid5 + phd + ment, data = tib.hurdle)
}
tib.gee <- tibble::as_tibble(data.gee)
if(requireNamespace("gee")){
mod.gee.tib <- gee::gee(depression ~ diagnose + drug*time, data = tib.gee,
id = id, family = binomial, corstr = "independence")
}
if(requireNamespace("geepack")){
mod.geeglm.tib <- geepack::geeglm(depression ~ diagnose + drug*time,
data = tib.gee, id = id, family = binomial,
corstr = "independence")
}
if(requireNamespace("glmtoolbox")){
mod.glmgee.tib <- glmtoolbox::glmgee(depression ~ diagnose + drug*time,
data = tib.gee, id = id, family = binomial,
corstr = "independence")
}
if(requireNamespace("lme4")){
lmerModtibdat = tibble::as_tibble(lme4::sleepstudy)
mod.lmerMod.tib <- lme4::lmer(Reaction ~ Days + (Days | Subject),
lmerModtibdat)
}}
# glmmTMB
if(requireNamespace("glmmTMB")){
model_glmmTMB <- glmmTMB::glmmTMB(Reaction ~ Days + (1 | Subject), data = lme4::sleepstudy, family = gaussian())
}
# Models with transformed responses -- used to test that bootstrap (resi_stat)
# can update() the reduced model on bootstrap data that has the raw columns
# (e.g. 'charges') but NOT the pre-computed transformed column ('log10(charges)').
mod.lm.log10 <- lm(log10(charges) ~ region + age + bmi + sex, data = data)
mod.glm.I <- glm(I(charges > 15000) ~ region + age + bmi + sex,
data = data, family = binomial())
## tests
test_that("Specifying non-allowed vcov produces warning",{
expect_warning(resi(mod.nls, data = data.nls, nboot = 1, vcovfunc = sandwich::vcovHC),
"Sandwich vcov function not applicable for nls model type, vcovfunc set to regtools::nlshc")
if(requireNamespace("survival")){
expect_warning(resi(mod.surv, data = data.surv, nboot = 1, vcovfunc = sandwich::vcovHC),
"vcovfunc argument ignored for survreg objects")
expect_warning(resi(mod.coxph, data = data.surv, nboot = 1, vcovfunc = sandwich::vcovHC),
"vcovfunc argument ignored for coxph objects")}
})
test_that("boot.method = 'bayes' only works for lm and nls", {
expect_true(resi(mod.lm, nboot = 1, boot.method = "bayes")$boot.method == "bayes")
expect_true(resi(mod.nls, data = data.nls, nboot = 1, boot.method = "bayes")$boot.method == "bayes")
expect_error(resi(mod, nboot = 1, boot.method = "bayes"))
if(requireNamespace("survival")){
expect_error(resi(mod.surv, data = data.surv, nboot = 1, boot.method = "bayes"))
expect_error(resi(mod.coxph, data = data.surv, nboot = 1, boot.method = "bayes"))}
if(requireNamespace("pscl")){
expect_error(resi(mod.hurdle, nboot = 1, boot.method = "bayes"))
expect_error(resi(mod.zinf, nboot = 1, boot.method = "bayes"))}
})
test_that("data is needed for certain models", {
# mod.s now uses ci.method='qf' by default, which succeeds without data;
# use ci.method='boot' to test that bootstrapping fails without data
expect_error(resi(mod.s, nboot = 1, ci.method = "boot"))
# mod.lm.s now uses ci.method='qf' by default, which succeeds without data;
# use ci.method='boot' to test that bootstrapping fails without data
expect_error(resi(mod.lm.s, nboot = 1, ci.method = "boot"))
expect_error(resi(mod.nls, nboot = 1))
if(requireNamespace("survival")){
expect_error(resi(mod.surv, nboot = 1))
expect_error(resi(mod.coxph, nboot = 1))}
if(requireNamespace("gee")){
expect_error(resi(mod.gee))}
})
# using different boot.results storage now
# test_that("boot.results stores correctly",{
# expect_equal(colnames(resi(mod.log, nboot = 2, store.boot = TRUE)$boot.results$t), c("Overall", "(Intercept)", "age", "regionnorthwest", "regionsoutheast", "regionsouthwest", "age", "region"))
# })
test_that("resi produces the correct estimates", {
expect_equal(unname(resi(mod, nboot = 1)$estimates), c(0.35982590, -0.06733537, -0.02670248, -0.03341748, -0.00246893, 0.14980504, 0.15238719,
0.05839381, 0.01667149, 0.03610640, -0.01497171, 0.03655364, 0.29533571, 0.14991488,
0.05159896, 0.01617303), tolerance = 1e-07)
expect_equal(unname(resi(mod.lm, nboot = 1)$estimates), c(0.3595407548, -0.0672973369,
-0.0266873974, -0.0333986021,
-0.0024675352, 0.1497204142,
0.1523011075, 0.0583608197,
0.0166620711, 0.0360860037,
-0.0149632506, 0.0364798507,
0.2951113263, 0.1497981921,
0.0515491712, 0.0160560332), tolerance = 1e-07)
expect_equal(unname(resi(mod.s, nboot = 1, data = data)$estimates[c(1, 6, 11, 24)]),
c(0.35634034, 0.06143486, 0.00748577, 0.05339006), tolerance = 1e-07)
if(requireNamespace("survival")){
expect_equal(unname(resi(mod.surv, nboot = 1, data = data.surv)$estimates),
c(0.194002944, 0.515785774, -0.080800275, 0.200549139, 0.105353954, -0.276392759,
0.046080344, 0.189247644, 0.081817903), tolerance = 1e-07)
expect_equal(unname(resi(mod.coxph, nboot = 1, data = data.surv)$estimates),
c(0.224354226, 0.136138740, -0.213293162, 0.008641209, 0.117732151,
0.202042262, 0.000000000))
}
if(requireNamespace("pscl")){
expect_equal(unname(resi(mod.hurdle, nboot = 1)$estimates[c(1, 2, 5, 8, 10)]),
c(0.28063439, 0.12770489, -0.06898884, 0.02686305, 0.06039791),
tolerance = 1e-07)
expect_equal(unname(resi(mod.zinf, nboot = 1)$estimates[c(1, 2, 5, 8, 10)]),
c(0.23725614, 0.11892936, -0.07014934, -0.03481673, -0.03461235),
tolerance = 1e-05)
}
if(requireNamespace("gee")){
expect_equal(unname(resi(mod.gee, nboot = 10, data = data.gee)$coefficients[,'L-RESI']),
c(-0.02585617, -0.48811210, 0.01414410, 0.56177861, -0.29398258),
tolerance = 1e-07)
expect_equal(unname(resi(mod.gee, nboot = 10, data = data.gee)$coefficients[,'CS-RESI']),
c(-0.007006305,-0.121182453,0.003769394,0.119402133,-0.069223594),
tolerance = 1e-07)
}
if(requireNamespace("geepack")){
expect_equal(unname(resi(mod.geeglm, nboot = 10, data = data.gee)$coefficients[,'L-RESI']),
c(-0.02585617, -0.48811210, 0.01414410, 0.56177861, -0.29398258),
tolerance = 1e-07)
expect_equal(unname(resi(mod.geeglm, nboot = 10, data = data.gee)$coefficients[,'CS-RESI']),
c(-0.007006305,-0.121182453,0.003769394,0.119402133,-0.069223594),
tolerance = 1e-07)
}
# if(requireNamespace("glmtoolbox")){
# for some reason test is saying mod.glmgee not found, can't resolve. manually tested
# # tests don't seem to like passing data argument. I think glmtoolbox deals with environments strangely
# mod.glmgee = glmtoolbox::glmgee(depression ~ diagnose + drug*time, data = data.gee,
# id = id, family = binomial, corstr = "independence")
# expect_equal(unname(resi(mod.glmgee, nboot = 10)$coefficients[,'L-RESI']),
# c(-0.02585617, -0.48811210, 0.01414410, 0.56177861, -0.29398258),
# tolerance = 1e-07)
# expect_equal(unname(resi(mod.glmgee, nboot = 10)$coefficients[,'CS-RESI']),
# c(-0.007006305,-0.121182453,0.003769394,0.119402133,-0.069223594),
# tolerance = 1e-07)
#
# # # SNV: temporary pasted in for now
# # library(glmtoolbox)
# # data(spruce)
# # m1 = glmgee(size ~ poly(days,4) + treat, id=tree, family=Gamma(log), data=spruces)
# # debug(resi.glmgee)
# # glmtoolboxRESI = resi(m1, nboot=10, data=spruces)
# # glmtoolboxRESI = resi_pe(m1, nboot=10, data=spruces)
# # ggmod = geeglm(size ~ poly(days,4) + treat, id=tree, family=Gamma(log), data=spruces)
# # geeglmRESI = resi_pe(ggmod)
# }
expect_equal(unname(suppressWarnings(resi(mod.lme, nboot = 10)$coefficients[,'RESI'])),
c(3.659090, 1.739166, 0.512371), tolerance = 1e-07)
if(requireNamespace("lme4")){
expect_equal(unname(suppressWarnings(resi(mod.lmerMod, nboot = 10)$coefficients[,'L-RESI'])),
c(8.434942, 1.533073), tolerance = 1e-07)}
if(requireNamespace("glmmTMB")){
expect_equal(unname(suppressWarnings(resi(model_glmmTMB, nboot = 10)$coefficients[,'RESI'])),
c( 6.22903, 3.06823), tolerance = 1e-05)
}
})
test_that("RESI estimates are in between the confidence limits", {
resi.obj = resi(mod, nboot = 500, store.boot = TRUE)
expect_true(all(resi.obj$coefficients$RESI >= resi.obj$coefficients$`2.5%`) & all(resi.obj$coefficients$RESI <= resi.obj$coefficients$`97.5%`))
expect_true(all(resi.obj$anova$RESI >= resi.obj$anova$`2.5%`) & all(resi.obj$anova$RESI <= resi.obj$anova$`97.5%`))
An.obj = car::Anova(resi.obj, alpha = 0.01)
expect_true(all(An.obj$RESI >= An.obj$`0.05%`) & all(An.obj$RESI <= An.obj$`99.5%`))
# ci.method='boot' needed: resi.lm now defaults to 'qf' which skips bootstrap
# and store.boot=TRUE would have no effect, breaking Anova.resi at a new alpha
resi.obj = resi(mod.lm, nboot = 500, store.boot = TRUE, ci.method = "boot")
expect_true(all(resi.obj$coefficients$RESI >= resi.obj$coefficients$`2.5%`) & all(resi.obj$coefficients$RESI <= resi.obj$coefficients$`97.5%`))
expect_true(all(resi.obj$anova[1:5, "RESI"] >= resi.obj$anova[1:5, "2.5%"]) & all(resi.obj$anova[1:5, "RESI"] <= resi.obj$anova[1:5, "RESI"]))
An.obj = car::Anova(resi.obj, alpha = 0.01)
expect_true(all(An.obj$RESI[1:5] >= An.obj$`0.05%`[1:5]) & all(An.obj$RESI[1:5] <= An.obj$`99.5%`[1:5]))
resi.obj = resi(mod.nls, data = data.nls)
expect_true(all(resi.obj$coefficients$RESI >= resi.obj$coefficients$`2.5%`) & all(resi.obj$coefficients$RESI <= resi.obj$coefficients$`97.5%`))
if(requireNamespace("survival")){
resi.obj = resi(mod.coxph, data = data.surv, nboot = 500)
expect_true(all(resi.obj$coefficients$RESI >= resi.obj$coefficients$`2.5%`) & all(resi.obj$coefficients$RESI <= resi.obj$coefficients$`97.5%`))
expect_true(all(resi.obj$anova$RESI >= resi.obj$anova$`2.5%`) & all(resi.obj$anova$RESI <= resi.obj$anova$`97.5%`))}
# commented for time
# resi.obj = resi(mod.hurdle, nboot = 500)
# expect_true(all(resi.obj$coefficients$RESI >= resi.obj$coefficients$`2.5%`) & all(resi.obj$coefficients$RESI <= resi.obj$coefficients$`97.5%`))
if(requireNamespace("geepack")){
resi.obj = resi(mod.geeglm, nboot = 500)
expect_true(all(resi.obj$coefficients$`L-RESI` >= resi.obj$coefficients$`L 2.5%`) & all(resi.obj$coefficients$RESI <= resi.obj$coefficients$`L 97.5%`))
expect_true(all(resi.obj$coefficients$`CS-RESI` >= resi.obj$coefficients$`CS 2.5%`) & all(resi.obj$coefficients$RESI <= resi.obj$coefficients$`CS 97.5%`))
expect_true(all(resi.obj$anova$`L-RESI` >= resi.obj$anova$`L 2.5%`) & all(resi.obj$anova$RESI <= resi.obj$anova$`L 97.5%`))
expect_true(all(resi.obj$anova$`CS-RESI` >= resi.obj$anova$`CS 2.5%`) & all(resi.obj$anova$RESI <= resi.obj$anova$`CS 97.5%`))}
})
# Regression test: bootstrap must succeed for models with transformed responses.
# Bug: resi.default previously fit the intercept-only reduced model using a
# backtick-quoted column name (e.g. `log10(charges)` ~ 1) from model.frame,
# which broke resi_stat's update(model.reduced, data = boot.data) because
# boot.data has 'charges' but not a pre-computed 'log10(charges)' column.
# Every replicate returned NA -> "Bootstrapping failed" error.
test_that("bootstrap succeeds for models with transformed responses", {
set.seed(1)
# lm with log-transformed response
r.lm <- resi(mod.lm.log10, data = data, nboot = 20, ci.method = "boot")
expect_s3_class(r.lm, "resi")
expect_false(anyNA(r.lm$coefficients$RESI))
expect_true(all(r.lm$coefficients$RESI >= r.lm$coefficients$`2.5%`))
expect_true(all(r.lm$coefficients$RESI <= r.lm$coefficients$`97.5%`))
# glm with I() binary response
r.glm <- suppressWarnings(
resi(mod.glm.I, data = data, nboot = 20, ci.method = "boot")
)
expect_s3_class(r.glm, "resi")
expect_false(anyNA(r.glm$coefficients$RESI))
expect_true(all(r.glm$coefficients$RESI >= r.glm$coefficients$`2.5%`))
expect_true(all(r.glm$coefficients$RESI <= r.glm$coefficients$`97.5%`))
})
if(requireNamespace("lme4")){
test_that("lmerMod: L-RESI and CS-RESI point estimates are correct", {
pe = suppressWarnings(resi_pe(mod.lmerMod))
expect_equal(unname(pe$coefficients[,'L-RESI']),
c(8.434942, 1.533073), tolerance = 1e-07)
expect_equal(unname(pe$coefficients[,'CS-RESI']),
c(3.6680844, 0.6185898), tolerance = 1e-06)
expect_equal(unname(pe$anova[,'L-RESI']), 1.533073, tolerance = 1e-06)
expect_equal(unname(pe$anova[,'CS-RESI']), 0.5719247, tolerance = 1e-06)
})
if(requireNamespace("geepack")){
test_that("geeglm (exchangeable, positive rho) L-RESI matches lmerMod (CR0) L-RESI", {
sdata = lme4::sleepstudy
sdata$subjid = as.integer(sdata$Subject)
mod_lmer_ri = lme4::lmer(Reaction ~ Days + (1 | Subject), sdata)
mod_gee_exch = geepack::geeglm(Reaction ~ Days, data = sdata, id = subjid,
family = gaussian, corstr = "exchangeable")
pe_lmer = suppressWarnings(resi_pe(mod_lmer_ri, vcov.args = list(type = "CR0")))
pe_gee = resi_pe(mod_gee_exch)
# for positive exchangeable correlation (alpha ~ 0.58), GEE and CR0 sandwich
# give the same L-RESI up to finite-sample differences
expect_equal(unname(pe_lmer$coefficients[,'L-RESI']),
unname(pe_gee$coefficients[,'L-RESI']), tolerance = 0.05)
# CS-RESI is identical: both compute the same weighted independence lm + HC0
expect_equal(unname(pe_lmer$coefficients[,'CS-RESI']),
unname(pe_gee$coefficients[,'CS-RESI']), tolerance = 1e-07)
})}}
if(requireNamespace("gee") & requireNamespace("geepack")){
test_that("vcov methods same for gee and geeglm",{
expect_equal(unname(mod.gee$robust.variance), mod.geeglm$geese$vbeta)
modg1 = glm(formula = formula(mod.gee), family = mod.gee$family, data = data.gee,
contrasts = mod.gee$contrasts)
modg2 = glm(formula = formula(mod.geeglm), family = mod.geeglm$family, data = data.gee,
contrasts = mod.geeglm$contrasts)
expect_equal(sandwich::vcovHC(modg1, type = "HC0"), sandwich::vcovHC(modg2, type = "HC0"))
})
test_that("boot.results same (approx) for gee and geeglm",{
set.seed(123)
resi.obj = resi(mod.geeglm, nboot = 5, store.boot = T, anova = F)
set.seed(123)
resi.obj2 = resi(mod.gee, data = data.gee, nboot = 5, store.boot = T)
expect_equal(resi.obj$boot.results$t[,-1], resi.obj2$boot.results$t, tolerance = 1e-09)
})}
if(requireNamespace("glmtoolbox") & requireNamespace("geepack")){
test_that("vcov and estimates are the same for glmgee and geeglm",{
expect_equal(unname(mod.glmgee$R), mod.geeglm$geese$vbeta)
})
test_that("boot.results same (approx) for gee and geeglm",{
set.seed(123)
resi.obj = resi(mod.geeglm, nboot = 5, store.boot = T, anova = F)
set.seed(123)
resi.obj2 = resi(mod.gee, data = data.gee, nboot = 5, store.boot = T)
expect_equal(resi.obj$boot.results$t[,-1], resi.obj2$boot.results$t, tolerance = 1e-09)
})}
test_that("unbiased = FALSE returns same abs. RESI as Chi-sq/F",{
resi.obj = resi(mod, nboot = 1, unbiased = FALSE)
expect_equal(abs(resi.obj$coefficients[6:7, 'RESI']), resi.obj$anova[3:4, 'RESI'], tolerance = 1e-03)
if(requireNamespace("survival")){
resi.obj = resi(mod.surv, data = data.surv, nboot = 1, unbiased = FALSE)
expect_equal(abs(resi.obj$coefficients[2:4, 'RESI']), resi.obj$anova[1:3, 'RESI'], tolerance = 1e-07)
resi.obj = resi(mod.coxph, data = data.surv, nboot = 1, unbiased = FALSE)
expect_equal(abs(resi.obj$coefficients[1:3, 'RESI']), resi.obj$anova[1:3, 'RESI'], tolerance = 1e-07)}
resi.obj = resi(mod.lm, nboot = 1, unbiased = FALSE)
expect_equal(abs(resi.obj$coefficients[6:7,'RESI']), resi.obj$anova[3:4,'RESI'], tolerance = 1e-07)
})
test_that("Using a matrix for data works", {
expect_true(class(resi(mod.nls, nboot = 10, data = as.matrix(data.nls))) == "resi")
})
test_that("Specifying a reduced model only changes overall output", {
resi.obj = resi(mod, nboot = 1)
resi.objr = resi(mod, model.reduced = mod.r, nboot = 1)
expect_equal(resi.obj$coefficients$RESI, resi.objr$coefficients$RESI)
expect_equal(resi.obj$anova$RESI, resi.objr$anova$RESI)
expect_false(resi.obj$overall$RESI[2] == resi.objr$overall$RESI[2])
resi.obj = resi(mod.lm.s, data = data, nboot = 1)
resi.objr = resi(mod.lm.s, model.reduced = mod.lm.s.r, data = data, nboot = 1)
expect_equal(resi.obj$coefficients$RESI, resi.objr$coefficients$RESI)
expect_equal(resi.obj$anova$RESI, resi.objr$anova$RESI)
expect_false(resi.obj$overall$RESI[2] == resi.objr$overall$RESI[2])
if(requireNamespace("survival")){
resi.obj = resi(mod.surv, data = data.surv, nboot = 1)
resi.objr = resi(mod.surv, model.reduced = mod.surv.red, data = data.surv[which(!(1:nrow(data.surv)%in% mod.surv$na.action)),], nboot = 10)
expect_false(resi.obj$overall$RESI[2] == resi.objr$overall$RESI[2])}
if(requireNamespace("pscl")){
resi.obj = resi(mod.hurdle, nboot = 1)
resi.objr = resi(mod.hurdle, model.reduced = mod.hurdle.r, nboot = 1)
expect_equal(resi.obj$coefficients$RESI, resi.objr$coefficients$RESI)
expect_equal(resi.obj$anova$RESI, resi.objr$anova$RESI)
expect_false(resi.obj$overall$RESI[2] == resi.objr$overall$RESI[2])}
})
test_that("specifying additional vcov args works",{
expect_false(any(resi(mod, nboot = 1, vcov.args = list(type = "HC0"))$coefficients[,'RESI'] == resi(mod, nboot = 1)$coefficients[,'RESI']))
})
test_that("specifying additional Anova args works",{
expect_false(all(resi(mod, nboot = 1, Anova.args = list(type = "3"))$anova[2:6,'RESI'] == resi(mod, nboot = 1)$anova[,'RESI']))
if(requireNamespace("survival")){
expect_false(length(resi(mod.surv, data = data.surv, nboot = 1, Anova.args = list(type = "3"))$anova$RESI) == length(resi(mod.surv, data = data.surv, nboot = 1)$anova$RESI))}
})
test_that("vcovfunc = vcov changes naive.var to TRUE",{
expect_true(resi(mod, vcovfunc = vcov, nboot = 10)$naive.var == TRUE)
})
if(requireNamespace("tibble")){
test_that("tibbles work", {
expect_equal(unname(resi(mod.glm.tib, nboot = 1)$estimates), c(0.35982590, -0.06733537, -0.02670248, -0.03341748, -0.00246893, 0.14980504, 0.15238719,
0.05839381, 0.01667149, 0.03610640, -0.01497171, 0.03655364, 0.29533571, 0.14991488,
0.05159896, 0.01617303), tolerance = 1e-07)
expect_equal(unname(resi(mod.lm.tib, nboot = 1)$estimates), c(0.3595407548, -0.0672973369,
-0.0266873974, -0.0333986021,
-0.0024675352, 0.1497204142,
0.1523011075, 0.0583608197,
0.0166620711, 0.0360860037,
-0.0149632506, 0.0364798507,
0.2951113263, 0.1497981921,
0.0515491712, 0.0160560332))
if(requireNamespace("survival")){
expect_equal(unname(resi(mod.surv.tib, nboot = 1, data = data.surv)$estimates), c(0.194002944, 0.515785774, -0.080800275, 0.200549139, 0.105353954, -0.276392759,
0.046080344, 0.189247644, 0.081817903), tolerance = 1e-07)
expect_equal(unname(resi(mod.coxph.tib, nboot = 1, data = data.surv)$estimates), c(0.224354226, 0.136138740, -0.213293162, 0.008641209, 0.117732151, 0.202042262, 0.000000000))
}
if(requireNamespace("pscl")){
expect_equal(unname(resi(mod.hurdle.tib, nboot = 1)$estimates[c(1, 2, 5, 8, 10)]), c(0.28063439, 0.12770489, -0.06898884, 0.02686305, 0.06039791), tolerance = 1e-07)
expect_equal(unname(resi(mod.zinf.tib, nboot = 1)$estimates[c(1, 2, 5, 8, 10)]), c(0.23725614, 0.11892936, -0.07014934, -0.03481673, -0.03461235), tolerance = 1e-05)
}
if(requireNamespace("gee")){
expect_equal(unname(resi(mod.gee.tib, nboot = 10, data = data.gee)$coefficients[,'L-RESI']), c(-0.02585617, -0.48811210, 0.01414410, 0.56177861, -0.29398258), tolerance = 1e-07)
}
if(requireNamespace("geepack")){
expect_equal(unname(resi(mod.geeglm.tib, nboot = 10, data = data.gee)$coefficients[,'L-RESI']), c(-0.02585617, -0.48811210, 0.01414410, 0.56177861, -0.29398258), tolerance = 1e-07)
}
# if(requireNamespace("glmtoolbox")){
# expect_equal(unname(resi(mod.glmgee.tib, nboot = 10)$coefficients[,'L-RESI']), c(-0.02585617, -0.48811210, 0.01414410, 0.56177861, -0.29398258), tolerance = 1e-07)
# }
if(requireNamespace("lme4")){
expect_equal(unname(suppressWarnings(resi(mod.lmerMod.tib, nboot = 10))$coefficients[,'L-RESI']),c(8.434942, 1.533073), tolerance = 1e-07)
}
})}
test_that("intercept-only full model works and doesn't have overall element", {
expect_true(is.null(resi(mod.glm.int, nboot = 1, Anova.args = list(type = 3))$overall))
expect_true(is.null(resi(mod.lm.int, nboot = 1, Anova.args = list(type = 3))$overall))
if(requireNamespace("survival")){
expect_true(is.null(resi(mod.surv.int, data = data.surv, nboot = 1, Anova.args = list(type = 3))$overall))
}
if(requireNamespace("pscl")){
expect_true(is.null(resi(mod.hurdle.int, nboot = 1, Anova.args = list(type = 3))$overall))
expect_true(is.null(resi(mod.zinf.int, nboot = 1, Anova.args = list(type = 3))$overall))
}
# to be implemented
# if(requireNamespace("geepack")){
# expect_true(is.null(resi(mod.geeglm.int, nboot = 1)$overall))
# }
})
test_that("same model.full and model.reduced produces error",{
expect_error(resi(model.full = mod, model.reduced = mod, nboot = 1))
})
# important because some resi_pe methods use waldtest and some use wald.test
test_that("wald.test and waldtest still consistent",{
expect_equal(unname(lmtest::waldtest(mod.lm, vcov = vcovHC, test = "Chisq")$Chisq[2]),
unname(aod::wald.test(vcovHC(mod.lm), coef(mod.lm),
Terms = 2:nrow(vcovHC(mod.lm)))$result$chi2["chi2"]))
})
# issue from github
test_that("overall = FALSE doesn't produce an error", {
expect_equal(unname(resi(mod.lm, nboot = 1, overall = F, coefficients = F)$estimates),
c(0.03647985, 0.29511133, 0.14979819, 0.05154917, 0.01605603), tolerance = 1e-05)
})
# robustbase: lmrob and glmrob
if (requireNamespace("robustbase", quietly = TRUE)) {
set.seed(7)
n_rob <- 150
df_rob <- data.frame(
x1 = rnorm(n_rob),
x2 = factor(sample(c("a", "b", "c"), n_rob, replace = TRUE)),
y = rnorm(n_rob)
)
df_rob$y <- df_rob$y + 2 * df_rob$x1 +
ifelse(df_rob$x2 == "b", 1.5, ifelse(df_rob$x2 == "c", -1, 0))
df_rob$y[1:5] <- df_rob$y[1:5] + 15 # outliers
mod.lmrob <- robustbase::lmrob(y ~ x1 + x2, data = df_rob)
mod.lmrob.r <- robustbase::lmrob(y ~ x1, data = df_rob)
mod.lmrob.int <- robustbase::lmrob(y ~ 1, data = df_rob)
set.seed(7)
n_grob <- 150
df_glmrob <- data.frame(
x1 = rnorm(n_grob),
x2 = factor(sample(c("a", "b", "c"), n_grob, replace = TRUE))
)
lp_grob <- 1.5 * df_glmrob$x1 +
ifelse(df_glmrob$x2 == "b", 1, -0.5)
df_glmrob$y <- rbinom(n_grob, 1, plogis(lp_grob))
mod.glmrob <- suppressWarnings(robustbase::glmrob(y ~ x1 + x2, family = binomial, data = df_glmrob))
mod.glmrob.r <- suppressWarnings(robustbase::glmrob(y ~ x1, family = binomial, data = df_glmrob))
mod.glmrob.int <- suppressWarnings(robustbase::glmrob(y ~ 1, family = binomial, data = df_glmrob))
test_that("lmrob: resi_pe produces correct estimates", {
pe <- suppressWarnings(resi_pe(mod.lmrob))
expect_equal(unname(pe$estimates),
c(2.1317677, 0.0321117, 1.7228418, 0.5869813, -0.2935287,
1.7298297, 0.9128559),
tolerance = 1e-06)
expect_equal(unname(pe$anova$RESI), c(1.7298297, 0.9128559), tolerance = 1e-06)
expect_equal(unname(pe$coefficients$RESI),
c(0.0321117, 1.7228418, 0.5869813, -0.2935287), tolerance = 1e-06)
})
test_that("lmrob: resi() returns RESI within CI bounds", {
set.seed(1)
ro <- suppressWarnings(resi(mod.lmrob, nboot = 100))
expect_true(all(ro$coefficients$RESI >= ro$coefficients$`2.5%`))
expect_true(all(ro$coefficients$RESI <= ro$coefficients$`97.5%`))
expect_true(all(ro$anova$RESI >= ro$anova$`2.5%`))
expect_true(all(ro$anova$RESI <= ro$anova$`97.5%`))
})
test_that("lmrob: specifying a reduced model only changes overall output", {
ro <- suppressWarnings(resi(mod.lmrob, nboot = 1))
ror <- suppressWarnings(resi(mod.lmrob, model.reduced = mod.lmrob.r, nboot = 1))
expect_equal(ro$coefficients$RESI, ror$coefficients$RESI)
expect_equal(ro$anova$RESI, ror$anova$RESI)
expect_false(ro$overall$RESI[2] == ror$overall$RESI[2])
})
test_that("lmrob: intercept-only model has no overall element", {
expect_true(is.null(suppressWarnings(resi(mod.lmrob.int, nboot = 1, anova = FALSE))$overall))
})
test_that("glmrob: resi_pe produces correct estimates", {
pe2 <- suppressWarnings(resi_pe(mod.glmrob))
expect_equal(unname(pe2$estimates),
c(0.4340252, -0.1257804, 0.4254354, 0.2549712, 0.0386988,
0.4175268, 0.2551174),
tolerance = 1e-06)
expect_equal(unname(pe2$anova$RESI), c(0.4175268, 0.2551174), tolerance = 1e-06)
})
test_that("glmrob: resi() returns RESI within CI bounds", {
set.seed(1)
ro2 <- suppressWarnings(resi(mod.glmrob, data = df_glmrob, nboot = 100))
expect_true(all(ro2$coefficients$RESI >= ro2$coefficients$`2.5%`))
expect_true(all(ro2$coefficients$RESI <= ro2$coefficients$`97.5%`))
expect_true(all(ro2$anova$RESI >= ro2$anova$`2.5%`))
expect_true(all(ro2$anova$RESI <= ro2$anova$`97.5%`))
})
test_that("glmrob: vcovHC redirects to stats::vcov with warning", {
expect_warning(
resi_pe(mod.glmrob, vcovfunc = sandwich::vcovHC),
"sandwich::vcovHC is not supported for glmrob"
)
})
test_that("glmrob: intercept-only model has no overall element", {
expect_true(is.null(suppressWarnings(resi(mod.glmrob.int, data = df_glmrob, nboot = 1, anova = FALSE))$overall))
})
}
Try the RESI package in your browser
Any scripts or data that you put into this service are public.
RESI documentation built on July 25, 2026, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(car)
#library(splines)
library(sandwich)
#library(survival)
library(nlme)
#library(lme4)
library(tibble)
#library(geepack)
data = RESI::insurance
mod = glm(charges ~ region * age + bmi + sex, data = data)
mod.r = glm(charges ~ bmi +sex, data = data)
mod.lm = lm(charges ~ region * age + bmi + sex, data = data)
mod.log = glm(smoker ~ age + region, data = data, family = "binomial")
mod.glm.int = glm(smoker ~ 1, data = data, family = "binomial")
mod.lm.int = lm(charges ~ 1, data = data)
# splines
if(requireNamespace("splines")){
mod.s = glm(charges ~ region * splines::ns(age, df=3) + bmi + sex, data = data)
mod.s.r = glm(charges ~ region * splines::ns(age, df=3), data = data)
mod.lm.s = lm(charges ~ region*splines::ns(age, df = 3) + bmi + sex, data = data)
mod.lm.s.r = lm(charges ~ region * splines::ns(age, df=3), data = data)
}
# missing data models
data.nas = data
data.nas[12:17, 'bmi'] = NA
mod.na = lm(charges ~ region * age + bmi + sex, data = data.nas)
mod.log.na = glm(smoker ~ age + region + bmi, data = data.nas, family = "binomial")
#nls
x <- seq(0, 100, 1)
y<-((runif(1, 10, 20)*x)/(runif(1, 0, 10) + x)) +
rnorm(101, 0, 1)
data.nls = data.frame(X = x, Y = y)
mod.nls = nls(Y~a * X/(b + X), data = data.nls, start = list(a = 1, b = 2))
#survival
if(requireNamespace("survival")){
data.surv = survival::lung
mod.surv = survival::survreg(survival::Surv(time, status) ~ age + sex + ph.karno, data=data.surv,
dist='weibull', robust = TRUE)
mod.surv.red = survival::survreg(survival::Surv(time, status) ~ age, data=data.surv[which(!(1:nrow(data.surv)%in% mod.surv$na.action)),],
dist='weibull', robust = TRUE)
mod.surv.nr = survival::survreg(survival::Surv(time, status) ~ age + sex + ph.karno, data=data.surv,
dist='weibull', robust = FALSE)
mod.surv.int = survival::survreg(survival::Surv(time, status) ~ 1, data=data.surv,
dist='weibull', robust = TRUE)
mod.coxph = survival::coxph(survival::Surv(time, status) ~ age + sex + wt.loss, data=survival::lung, robust = TRUE)
mod.coxph.red = survival::coxph(survival::Surv(time, status) ~ age, data=survival::lung, robust = TRUE)
}
# hurdle
if(requireNamespace("pscl")){
data.hurdle = pscl::bioChemists
mod.hurdle = pscl::hurdle(art ~ fem + mar + kid5 + phd + ment | fem + mar + kid5 + phd + ment, data = data.hurdle)
mod.hurdle.r = pscl::hurdle(art ~ kid5 + phd | kid5 + phd, data = data.hurdle)
mod.hurdle.int = pscl::hurdle(art ~ 1, data = data.hurdle)
# zeroinfl
mod.zinf = pscl::zeroinfl(art ~ fem + mar + kid5 + phd + ment | fem + mar + kid5 + phd + ment, data = data.hurdle)
mod.zinf.r = pscl::zeroinfl(art ~ kid5 + phd | kid5 + phd, data = data.hurdle)
mod.zinf.int = pscl::zeroinfl(art ~ 1, data = data.hurdle)
}
data.gee = RESI::depression
# gee and lme
if(requireNamespace("gee")){
mod.gee = gee::gee(depression ~ diagnose + drug*time, data = data.gee,
id = id, family = binomial, corstr = "independence")
mod.gee.int = gee::gee(depression ~ 1, data = data.gee,
id = id, family = binomial, corstr = "independence")
}
if(requireNamespace("geepack")){
library(geepack)
mod.geeglm = geepack::geeglm(depression ~ diagnose + drug*time, data = data.gee,
id = id, family = binomial, corstr = "independence")
mod.geeglm.r = geepack::geeglm(depression ~ diagnose, data = data.gee, id = id,
family = binomial, corstr = "independence")
mod.geeglm.int = geepack::geeglm(depression ~ 1, data = data.gee, id = id,
family = binomial, corstr = "independence")
}
if(requireNamespace("glmtoolbox")){
#library(glmtoolbox)
mod.glmgee = glmtoolbox::glmgee(depression ~ diagnose + drug*time, data = data.gee,
id = id, family = binomial, corstr = "independence")
mod.glmgee.r = glmtoolbox::glmgee(depression ~ diagnose, data = data.gee, id = id,
family = binomial, corstr = "independence")
mod.glmgee.int = glmtoolbox::glmgee(depression ~ 1, data = data.gee, id = id,
family = binomial, corstr = "independence")
}
mod.lme = nlme::lme(distance ~ age + Sex, data = nlme::Orthodont, random = ~ 1)
mod.lme.int = nlme::lme(distance ~ 1, data = nlme::Orthodont, random = ~ 1)
if(requireNamespace("lme4")){
mod.lmerMod = lme4::lmer(Reaction ~ Days + (Days | Subject), lme4::sleepstudy)
data.lmer = nlme::Orthodont
data.lmer$nsex <- as.numeric(data.lmer$Sex=="Male")
data.lmer$nsexage <- with(data.lmer, nsex*age)
mod.lmerMod2 = lme4::lmer(distance ~ age + (age|Subject) + (0+nsex|Subject) +
(0 + nsexage|Subject), data=data.lmer)
mod.lmerMod.int = lme4::lmer(Reaction ~ 1 + (1 | Subject), lme4::sleepstudy)
}
# tibbles
if(requireNamespace("tibble")){
tib.dat <- tibble::as_tibble(data)
mod.glm.tib <- glm(charges ~ region * age + bmi + sex, data = tib.dat)
mod.lm.tib <- lm(charges ~ region * age + bmi + sex, data = tib.dat)
tib.nls <- tibble::as_tibble(data.nls)
mod.nls.tib <- nls(Y~a * X/(b + X), data = tib.nls, start = list(a = 1, b = 2))
if(requireNamespace("survival")){
tib.surv <- tibble::as_tibble(data.surv)
mod.surv.tib <- survival::survreg(survival::Surv(time, status) ~ age + sex + ph.karno, data=tib.surv,
dist='weibull', robust = TRUE)
mod.coxph.tib <- survival::coxph(survival::Surv(time, status) ~ age + sex + wt.loss,
data=tibble::as_tibble(survival::lung), robust = TRUE)
}
if(requireNamespace("pscl")){
tib.hurdle <- tibble::as_tibble(data.hurdle)
mod.hurdle.tib <- pscl::hurdle(art ~ fem + mar + kid5 + phd + ment | fem + mar +
kid5 + phd + ment, data = tib.hurdle)
mod.zinf.tib <- pscl::zeroinfl(art ~ fem + mar + kid5 + phd + ment | fem + mar +
kid5 + phd + ment, data = tib.hurdle)
}
tib.gee <- tibble::as_tibble(data.gee)
if(requireNamespace("gee")){
mod.gee.tib <- gee::gee(depression ~ diagnose + drug*time, data = tib.gee,
id = id, family = binomial, corstr = "independence")
}
if(requireNamespace("geepack")){
mod.geeglm.tib <- geepack::geeglm(depression ~ diagnose + drug*time,
data = tib.gee, id = id, family = binomial,
corstr = "independence")
}
if(requireNamespace("glmtoolbox")){
mod.glmgee.tib <- glmtoolbox::glmgee(depression ~ diagnose + drug*time,
data = tib.gee, id = id, family = binomial,
corstr = "independence")
}
if(requireNamespace("lme4")){
lmerModtibdat = tibble::as_tibble(lme4::sleepstudy)
mod.lmerMod.tib <- lme4::lmer(Reaction ~ Days + (Days | Subject),
lmerModtibdat)
}}
# glmmTMB
if(requireNamespace("glmmTMB")){
model_glmmTMB <- glmmTMB::glmmTMB(Reaction ~ Days + (1 | Subject), data = lme4::sleepstudy, family = gaussian())
}
# Models with transformed responses -- used to test that bootstrap (resi_stat)
# can update() the reduced model on bootstrap data that has the raw columns
# (e.g. 'charges') but NOT the pre-computed transformed column ('log10(charges)').
mod.lm.log10 <- lm(log10(charges) ~ region + age + bmi + sex, data = data)
mod.glm.I <- glm(I(charges > 15000) ~ region + age + bmi + sex,
data = data, family = binomial())
## tests
test_that("Specifying non-allowed vcov produces warning",{
expect_warning(resi(mod.nls, data = data.nls, nboot = 1, vcovfunc = sandwich::vcovHC),
"Sandwich vcov function not applicable for nls model type, vcovfunc set to regtools::nlshc")
if(requireNamespace("survival")){
expect_warning(resi(mod.surv, data = data.surv, nboot = 1, vcovfunc = sandwich::vcovHC),
"vcovfunc argument ignored for survreg objects")
expect_warning(resi(mod.coxph, data = data.surv, nboot = 1, vcovfunc = sandwich::vcovHC),
"vcovfunc argument ignored for coxph objects")}
})
test_that("boot.method = 'bayes' only works for lm and nls", {
expect_true(resi(mod.lm, nboot = 1, boot.method = "bayes")$boot.method == "bayes")
expect_true(resi(mod.nls, data = data.nls, nboot = 1, boot.method = "bayes")$boot.method == "bayes")
expect_error(resi(mod, nboot = 1, boot.method = "bayes"))
if(requireNamespace("survival")){
expect_error(resi(mod.surv, data = data.surv, nboot = 1, boot.method = "bayes"))
expect_error(resi(mod.coxph, data = data.surv, nboot = 1, boot.method = "bayes"))}
if(requireNamespace("pscl")){
expect_error(resi(mod.hurdle, nboot = 1, boot.method = "bayes"))
expect_error(resi(mod.zinf, nboot = 1, boot.method = "bayes"))}
})
test_that("data is needed for certain models", {
# mod.s now uses ci.method='qf' by default, which succeeds without data;
# use ci.method='boot' to test that bootstrapping fails without data
expect_error(resi(mod.s, nboot = 1, ci.method = "boot"))
# mod.lm.s now uses ci.method='qf' by default, which succeeds without data;
# use ci.method='boot' to test that bootstrapping fails without data
expect_error(resi(mod.lm.s, nboot = 1, ci.method = "boot"))
expect_error(resi(mod.nls, nboot = 1))
if(requireNamespace("survival")){
expect_error(resi(mod.surv, nboot = 1))
expect_error(resi(mod.coxph, nboot = 1))}
if(requireNamespace("gee")){
expect_error(resi(mod.gee))}
})
# using different boot.results storage now
# test_that("boot.results stores correctly",{
# expect_equal(colnames(resi(mod.log, nboot = 2, store.boot = TRUE)$boot.results$t), c("Overall", "(Intercept)", "age", "regionnorthwest", "regionsoutheast", "regionsouthwest", "age", "region"))
# })
test_that("resi produces the correct estimates", {
expect_equal(unname(resi(mod, nboot = 1)$estimates), c(0.35982590, -0.06733537, -0.02670248, -0.03341748, -0.00246893, 0.14980504, 0.15238719,
0.05839381, 0.01667149, 0.03610640, -0.01497171, 0.03655364, 0.29533571, 0.14991488,
0.05159896, 0.01617303), tolerance = 1e-07)
expect_equal(unname(resi(mod.lm, nboot = 1)$estimates), c(0.3595407548, -0.0672973369,
-0.0266873974, -0.0333986021,
-0.0024675352, 0.1497204142,
0.1523011075, 0.0583608197,
0.0166620711, 0.0360860037,
-0.0149632506, 0.0364798507,
0.2951113263, 0.1497981921,
0.0515491712, 0.0160560332), tolerance = 1e-07)
expect_equal(unname(resi(mod.s, nboot = 1, data = data)$estimates[c(1, 6, 11, 24)]),
c(0.35634034, 0.06143486, 0.00748577, 0.05339006), tolerance = 1e-07)
if(requireNamespace("survival")){
expect_equal(unname(resi(mod.surv, nboot = 1, data = data.surv)$estimates),
c(0.194002944, 0.515785774, -0.080800275, 0.200549139, 0.105353954, -0.276392759,
0.046080344, 0.189247644, 0.081817903), tolerance = 1e-07)
expect_equal(unname(resi(mod.coxph, nboot = 1, data = data.surv)$estimates),
c(0.224354226, 0.136138740, -0.213293162, 0.008641209, 0.117732151,
0.202042262, 0.000000000))
}
if(requireNamespace("pscl")){
expect_equal(unname(resi(mod.hurdle, nboot = 1)$estimates[c(1, 2, 5, 8, 10)]),
c(0.28063439, 0.12770489, -0.06898884, 0.02686305, 0.06039791),
tolerance = 1e-07)
expect_equal(unname(resi(mod.zinf, nboot = 1)$estimates[c(1, 2, 5, 8, 10)]),
c(0.23725614, 0.11892936, -0.07014934, -0.03481673, -0.03461235),
tolerance = 1e-05)
}
if(requireNamespace("gee")){
expect_equal(unname(resi(mod.gee, nboot = 10, data = data.gee)$coefficients[,'L-RESI']),
c(-0.02585617, -0.48811210, 0.01414410, 0.56177861, -0.29398258),
tolerance = 1e-07)
expect_equal(unname(resi(mod.gee, nboot = 10, data = data.gee)$coefficients[,'CS-RESI']),
c(-0.007006305,-0.121182453,0.003769394,0.119402133,-0.069223594),
tolerance = 1e-07)
}
if(requireNamespace("geepack")){
expect_equal(unname(resi(mod.geeglm, nboot = 10, data = data.gee)$coefficients[,'L-RESI']),
c(-0.02585617, -0.48811210, 0.01414410, 0.56177861, -0.29398258),
tolerance = 1e-07)
expect_equal(unname(resi(mod.geeglm, nboot = 10, data = data.gee)$coefficients[,'CS-RESI']),
c(-0.007006305,-0.121182453,0.003769394,0.119402133,-0.069223594),
tolerance = 1e-07)
}
# if(requireNamespace("glmtoolbox")){
# for some reason test is saying mod.glmgee not found, can't resolve. manually tested
# # tests don't seem to like passing data argument. I think glmtoolbox deals with environments strangely
# mod.glmgee = glmtoolbox::glmgee(depression ~ diagnose + drug*time, data = data.gee,
# id = id, family = binomial, corstr = "independence")
# expect_equal(unname(resi(mod.glmgee, nboot = 10)$coefficients[,'L-RESI']),
# c(-0.02585617, -0.48811210, 0.01414410, 0.56177861, -0.29398258),
# tolerance = 1e-07)
# expect_equal(unname(resi(mod.glmgee, nboot = 10)$coefficients[,'CS-RESI']),
# c(-0.007006305,-0.121182453,0.003769394,0.119402133,-0.069223594),
# tolerance = 1e-07)
#
# # # SNV: temporary pasted in for now
# # library(glmtoolbox)
# # data(spruce)
# # m1 = glmgee(size ~ poly(days,4) + treat, id=tree, family=Gamma(log), data=spruces)
# # debug(resi.glmgee)
# # glmtoolboxRESI = resi(m1, nboot=10, data=spruces)
# # glmtoolboxRESI = resi_pe(m1, nboot=10, data=spruces)
# # ggmod = geeglm(size ~ poly(days,4) + treat, id=tree, family=Gamma(log), data=spruces)
# # geeglmRESI = resi_pe(ggmod)
# }
expect_equal(unname(suppressWarnings(resi(mod.lme, nboot = 10)$coefficients[,'RESI'])),
c(3.659090, 1.739166, 0.512371), tolerance = 1e-07)
if(requireNamespace("lme4")){
expect_equal(unname(suppressWarnings(resi(mod.lmerMod, nboot = 10)$coefficients[,'L-RESI'])),
c(8.434942, 1.533073), tolerance = 1e-07)}
if(requireNamespace("glmmTMB")){
expect_equal(unname(suppressWarnings(resi(model_glmmTMB, nboot = 10)$coefficients[,'RESI'])),
c( 6.22903, 3.06823), tolerance = 1e-05)
}
})
test_that("RESI estimates are in between the confidence limits", {
resi.obj = resi(mod, nboot = 500, store.boot = TRUE)
expect_true(all(resi.obj$coefficients$RESI >= resi.obj$coefficients$`2.5%`) & all(resi.obj$coefficients$RESI <= resi.obj$coefficients$`97.5%`))
expect_true(all(resi.obj$anova$RESI >= resi.obj$anova$`2.5%`) & all(resi.obj$anova$RESI <= resi.obj$anova$`97.5%`))
An.obj = car::Anova(resi.obj, alpha = 0.01)
expect_true(all(An.obj$RESI >= An.obj$`0.05%`) & all(An.obj$RESI <= An.obj$`99.5%`))
# ci.method='boot' needed: resi.lm now defaults to 'qf' which skips bootstrap
# and store.boot=TRUE would have no effect, breaking Anova.resi at a new alpha
resi.obj = resi(mod.lm, nboot = 500, store.boot = TRUE, ci.method = "boot")
expect_true(all(resi.obj$coefficients$RESI >= resi.obj$coefficients$`2.5%`) & all(resi.obj$coefficients$RESI <= resi.obj$coefficients$`97.5%`))
expect_true(all(resi.obj$anova[1:5, "RESI"] >= resi.obj$anova[1:5, "2.5%"]) & all(resi.obj$anova[1:5, "RESI"] <= resi.obj$anova[1:5, "RESI"]))
An.obj = car::Anova(resi.obj, alpha = 0.01)
expect_true(all(An.obj$RESI[1:5] >= An.obj$`0.05%`[1:5]) & all(An.obj$RESI[1:5] <= An.obj$`99.5%`[1:5]))
resi.obj = resi(mod.nls, data = data.nls)
expect_true(all(resi.obj$coefficients$RESI >= resi.obj$coefficients$`2.5%`) & all(resi.obj$coefficients$RESI <= resi.obj$coefficients$`97.5%`))
if(requireNamespace("survival")){
resi.obj = resi(mod.coxph, data = data.surv, nboot = 500)
expect_true(all(resi.obj$coefficients$RESI >= resi.obj$coefficients$`2.5%`) & all(resi.obj$coefficients$RESI <= resi.obj$coefficients$`97.5%`))
expect_true(all(resi.obj$anova$RESI >= resi.obj$anova$`2.5%`) & all(resi.obj$anova$RESI <= resi.obj$anova$`97.5%`))}
# commented for time
# resi.obj = resi(mod.hurdle, nboot = 500)
# expect_true(all(resi.obj$coefficients$RESI >= resi.obj$coefficients$`2.5%`) & all(resi.obj$coefficients$RESI <= resi.obj$coefficients$`97.5%`))
if(requireNamespace("geepack")){
resi.obj = resi(mod.geeglm, nboot = 500)
expect_true(all(resi.obj$coefficients$`L-RESI` >= resi.obj$coefficients$`L 2.5%`) & all(resi.obj$coefficients$RESI <= resi.obj$coefficients$`L 97.5%`))
expect_true(all(resi.obj$coefficients$`CS-RESI` >= resi.obj$coefficients$`CS 2.5%`) & all(resi.obj$coefficients$RESI <= resi.obj$coefficients$`CS 97.5%`))
expect_true(all(resi.obj$anova$`L-RESI` >= resi.obj$anova$`L 2.5%`) & all(resi.obj$anova$RESI <= resi.obj$anova$`L 97.5%`))
expect_true(all(resi.obj$anova$`CS-RESI` >= resi.obj$anova$`CS 2.5%`) & all(resi.obj$anova$RESI <= resi.obj$anova$`CS 97.5%`))}
})
# Regression test: bootstrap must succeed for models with transformed responses.
# Bug: resi.default previously fit the intercept-only reduced model using a
# backtick-quoted column name (e.g. `log10(charges)` ~ 1) from model.frame,
# which broke resi_stat's update(model.reduced, data = boot.data) because
# boot.data has 'charges' but not a pre-computed 'log10(charges)' column.
# Every replicate returned NA -> "Bootstrapping failed" error.
test_that("bootstrap succeeds for models with transformed responses", {
set.seed(1)
# lm with log-transformed response
r.lm <- resi(mod.lm.log10, data = data, nboot = 20, ci.method = "boot")
expect_s3_class(r.lm, "resi")
expect_false(anyNA(r.lm$coefficients$RESI))
expect_true(all(r.lm$coefficients$RESI >= r.lm$coefficients$`2.5%`))
expect_true(all(r.lm$coefficients$RESI <= r.lm$coefficients$`97.5%`))
# glm with I() binary response
r.glm <- suppressWarnings(
resi(mod.glm.I, data = data, nboot = 20, ci.method = "boot")
)
expect_s3_class(r.glm, "resi")
expect_false(anyNA(r.glm$coefficients$RESI))
expect_true(all(r.glm$coefficients$RESI >= r.glm$coefficients$`2.5%`))
expect_true(all(r.glm$coefficients$RESI <= r.glm$coefficients$`97.5%`))
})
if(requireNamespace("lme4")){
test_that("lmerMod: L-RESI and CS-RESI point estimates are correct", {
pe = suppressWarnings(resi_pe(mod.lmerMod))
expect_equal(unname(pe$coefficients[,'L-RESI']),
c(8.434942, 1.533073), tolerance = 1e-07)
expect_equal(unname(pe$coefficients[,'CS-RESI']),
c(3.6680844, 0.6185898), tolerance = 1e-06)
expect_equal(unname(pe$anova[,'L-RESI']), 1.533073, tolerance = 1e-06)
expect_equal(unname(pe$anova[,'CS-RESI']), 0.5719247, tolerance = 1e-06)
})
if(requireNamespace("geepack")){
test_that("geeglm (exchangeable, positive rho) L-RESI matches lmerMod (CR0) L-RESI", {
sdata = lme4::sleepstudy
sdata$subjid = as.integer(sdata$Subject)
mod_lmer_ri = lme4::lmer(Reaction ~ Days + (1 | Subject), sdata)
mod_gee_exch = geepack::geeglm(Reaction ~ Days, data = sdata, id = subjid,
family = gaussian, corstr = "exchangeable")
pe_lmer = suppressWarnings(resi_pe(mod_lmer_ri, vcov.args = list(type = "CR0")))
pe_gee = resi_pe(mod_gee_exch)
# for positive exchangeable correlation (alpha ~ 0.58), GEE and CR0 sandwich
# give the same L-RESI up to finite-sample differences
expect_equal(unname(pe_lmer$coefficients[,'L-RESI']),
unname(pe_gee$coefficients[,'L-RESI']), tolerance = 0.05)
# CS-RESI is identical: both compute the same weighted independence lm + HC0
expect_equal(unname(pe_lmer$coefficients[,'CS-RESI']),
unname(pe_gee$coefficients[,'CS-RESI']), tolerance = 1e-07)
})}}
if(requireNamespace("gee") & requireNamespace("geepack")){
test_that("vcov methods same for gee and geeglm",{
expect_equal(unname(mod.gee$robust.variance), mod.geeglm$geese$vbeta)
modg1 = glm(formula = formula(mod.gee), family = mod.gee$family, data = data.gee,
contrasts = mod.gee$contrasts)
modg2 = glm(formula = formula(mod.geeglm), family = mod.geeglm$family, data = data.gee,
contrasts = mod.geeglm$contrasts)
expect_equal(sandwich::vcovHC(modg1, type = "HC0"), sandwich::vcovHC(modg2, type = "HC0"))
})
test_that("boot.results same (approx) for gee and geeglm",{
set.seed(123)
resi.obj = resi(mod.geeglm, nboot = 5, store.boot = T, anova = F)
set.seed(123)
resi.obj2 = resi(mod.gee, data = data.gee, nboot = 5, store.boot = T)
expect_equal(resi.obj$boot.results$t[,-1], resi.obj2$boot.results$t, tolerance = 1e-09)
})}
if(requireNamespace("glmtoolbox") & requireNamespace("geepack")){
test_that("vcov and estimates are the same for glmgee and geeglm",{
expect_equal(unname(mod.glmgee$R), mod.geeglm$geese$vbeta)
})
test_that("boot.results same (approx) for gee and geeglm",{
set.seed(123)
resi.obj = resi(mod.geeglm, nboot = 5, store.boot = T, anova = F)
set.seed(123)
resi.obj2 = resi(mod.gee, data = data.gee, nboot = 5, store.boot = T)
expect_equal(resi.obj$boot.results$t[,-1], resi.obj2$boot.results$t, tolerance = 1e-09)
})}
test_that("unbiased = FALSE returns same abs. RESI as Chi-sq/F",{
resi.obj = resi(mod, nboot = 1, unbiased = FALSE)
expect_equal(abs(resi.obj$coefficients[6:7, 'RESI']), resi.obj$anova[3:4, 'RESI'], tolerance = 1e-03)
if(requireNamespace("survival")){
resi.obj = resi(mod.surv, data = data.surv, nboot = 1, unbiased = FALSE)
expect_equal(abs(resi.obj$coefficients[2:4, 'RESI']), resi.obj$anova[1:3, 'RESI'], tolerance = 1e-07)
resi.obj = resi(mod.coxph, data = data.surv, nboot = 1, unbiased = FALSE)
expect_equal(abs(resi.obj$coefficients[1:3, 'RESI']), resi.obj$anova[1:3, 'RESI'], tolerance = 1e-07)}
resi.obj = resi(mod.lm, nboot = 1, unbiased = FALSE)
expect_equal(abs(resi.obj$coefficients[6:7,'RESI']), resi.obj$anova[3:4,'RESI'], tolerance = 1e-07)
})
test_that("Using a matrix for data works", {
expect_true(class(resi(mod.nls, nboot = 10, data = as.matrix(data.nls))) == "resi")
})
test_that("Specifying a reduced model only changes overall output", {
resi.obj = resi(mod, nboot = 1)
resi.objr = resi(mod, model.reduced = mod.r, nboot = 1)
expect_equal(resi.obj$coefficients$RESI, resi.objr$coefficients$RESI)
expect_equal(resi.obj$anova$RESI, resi.objr$anova$RESI)
expect_false(resi.obj$overall$RESI[2] == resi.objr$overall$RESI[2])
resi.obj = resi(mod.lm.s, data = data, nboot = 1)
resi.objr = resi(mod.lm.s, model.reduced = mod.lm.s.r, data = data, nboot = 1)
expect_equal(resi.obj$coefficients$RESI, resi.objr$coefficients$RESI)
expect_equal(resi.obj$anova$RESI, resi.objr$anova$RESI)
expect_false(resi.obj$overall$RESI[2] == resi.objr$overall$RESI[2])
if(requireNamespace("survival")){
resi.obj = resi(mod.surv, data = data.surv, nboot = 1)
resi.objr = resi(mod.surv, model.reduced = mod.surv.red, data = data.surv[which(!(1:nrow(data.surv)%in% mod.surv$na.action)),], nboot = 10)
expect_false(resi.obj$overall$RESI[2] == resi.objr$overall$RESI[2])}
if(requireNamespace("pscl")){
resi.obj = resi(mod.hurdle, nboot = 1)
resi.objr = resi(mod.hurdle, model.reduced = mod.hurdle.r, nboot = 1)
expect_equal(resi.obj$coefficients$RESI, resi.objr$coefficients$RESI)
expect_equal(resi.obj$anova$RESI, resi.objr$anova$RESI)
expect_false(resi.obj$overall$RESI[2] == resi.objr$overall$RESI[2])}
})
test_that("specifying additional vcov args works",{
expect_false(any(resi(mod, nboot = 1, vcov.args = list(type = "HC0"))$coefficients[,'RESI'] == resi(mod, nboot = 1)$coefficients[,'RESI']))
})
test_that("specifying additional Anova args works",{
expect_false(all(resi(mod, nboot = 1, Anova.args = list(type = "3"))$anova[2:6,'RESI'] == resi(mod, nboot = 1)$anova[,'RESI']))
if(requireNamespace("survival")){
expect_false(length(resi(mod.surv, data = data.surv, nboot = 1, Anova.args = list(type = "3"))$anova$RESI) == length(resi(mod.surv, data = data.surv, nboot = 1)$anova$RESI))}
})
test_that("vcovfunc = vcov changes naive.var to TRUE",{
expect_true(resi(mod, vcovfunc = vcov, nboot = 10)$naive.var == TRUE)
})
if(requireNamespace("tibble")){
test_that("tibbles work", {
expect_equal(unname(resi(mod.glm.tib, nboot = 1)$estimates), c(0.35982590, -0.06733537, -0.02670248, -0.03341748, -0.00246893, 0.14980504, 0.15238719,
0.05839381, 0.01667149, 0.03610640, -0.01497171, 0.03655364, 0.29533571, 0.14991488,
0.05159896, 0.01617303), tolerance = 1e-07)
expect_equal(unname(resi(mod.lm.tib, nboot = 1)$estimates), c(0.3595407548, -0.0672973369,
-0.0266873974, -0.0333986021,
-0.0024675352, 0.1497204142,
0.1523011075, 0.0583608197,
0.0166620711, 0.0360860037,
-0.0149632506, 0.0364798507,
0.2951113263, 0.1497981921,
0.0515491712, 0.0160560332))
if(requireNamespace("survival")){
expect_equal(unname(resi(mod.surv.tib, nboot = 1, data = data.surv)$estimates), c(0.194002944, 0.515785774, -0.080800275, 0.200549139, 0.105353954, -0.276392759,
0.046080344, 0.189247644, 0.081817903), tolerance = 1e-07)
expect_equal(unname(resi(mod.coxph.tib, nboot = 1, data = data.surv)$estimates), c(0.224354226, 0.136138740, -0.213293162, 0.008641209, 0.117732151, 0.202042262, 0.000000000))
}
if(requireNamespace("pscl")){
expect_equal(unname(resi(mod.hurdle.tib, nboot = 1)$estimates[c(1, 2, 5, 8, 10)]), c(0.28063439, 0.12770489, -0.06898884, 0.02686305, 0.06039791), tolerance = 1e-07)
expect_equal(unname(resi(mod.zinf.tib, nboot = 1)$estimates[c(1, 2, 5, 8, 10)]), c(0.23725614, 0.11892936, -0.07014934, -0.03481673, -0.03461235), tolerance = 1e-05)
}
if(requireNamespace("gee")){
expect_equal(unname(resi(mod.gee.tib, nboot = 10, data = data.gee)$coefficients[,'L-RESI']), c(-0.02585617, -0.48811210, 0.01414410, 0.56177861, -0.29398258), tolerance = 1e-07)
}
if(requireNamespace("geepack")){
expect_equal(unname(resi(mod.geeglm.tib, nboot = 10, data = data.gee)$coefficients[,'L-RESI']), c(-0.02585617, -0.48811210, 0.01414410, 0.56177861, -0.29398258), tolerance = 1e-07)
}
# if(requireNamespace("glmtoolbox")){
# expect_equal(unname(resi(mod.glmgee.tib, nboot = 10)$coefficients[,'L-RESI']), c(-0.02585617, -0.48811210, 0.01414410, 0.56177861, -0.29398258), tolerance = 1e-07)
# }
if(requireNamespace("lme4")){
expect_equal(unname(suppressWarnings(resi(mod.lmerMod.tib, nboot = 10))$coefficients[,'L-RESI']),c(8.434942, 1.533073), tolerance = 1e-07)
}
})}
test_that("intercept-only full model works and doesn't have overall element", {
expect_true(is.null(resi(mod.glm.int, nboot = 1, Anova.args = list(type = 3))$overall))
expect_true(is.null(resi(mod.lm.int, nboot = 1, Anova.args = list(type = 3))$overall))
if(requireNamespace("survival")){
expect_true(is.null(resi(mod.surv.int, data = data.surv, nboot = 1, Anova.args = list(type = 3))$overall))
}
if(requireNamespace("pscl")){
expect_true(is.null(resi(mod.hurdle.int, nboot = 1, Anova.args = list(type = 3))$overall))
expect_true(is.null(resi(mod.zinf.int, nboot = 1, Anova.args = list(type = 3))$overall))
}
# to be implemented
# if(requireNamespace("geepack")){
# expect_true(is.null(resi(mod.geeglm.int, nboot = 1)$overall))
# }
})
test_that("same model.full and model.reduced produces error",{
expect_error(resi(model.full = mod, model.reduced = mod, nboot = 1))
})
# important because some resi_pe methods use waldtest and some use wald.test
test_that("wald.test and waldtest still consistent",{
expect_equal(unname(lmtest::waldtest(mod.lm, vcov = vcovHC, test = "Chisq")$Chisq[2]),
unname(aod::wald.test(vcovHC(mod.lm), coef(mod.lm),
Terms = 2:nrow(vcovHC(mod.lm)))$result$chi2["chi2"]))
})
# issue from github
test_that("overall = FALSE doesn't produce an error", {
expect_equal(unname(resi(mod.lm, nboot = 1, overall = F, coefficients = F)$estimates),
c(0.03647985, 0.29511133, 0.14979819, 0.05154917, 0.01605603), tolerance = 1e-05)
})
# robustbase: lmrob and glmrob
if (requireNamespace("robustbase", quietly = TRUE)) {
set.seed(7)
n_rob <- 150
df_rob <- data.frame(
x1 = rnorm(n_rob),
x2 = factor(sample(c("a", "b", "c"), n_rob, replace = TRUE)),
y = rnorm(n_rob)
)
df_rob$y <- df_rob$y + 2 * df_rob$x1 +
ifelse(df_rob$x2 == "b", 1.5, ifelse(df_rob$x2 == "c", -1, 0))
df_rob$y[1:5] <- df_rob$y[1:5] + 15 # outliers
mod.lmrob <- robustbase::lmrob(y ~ x1 + x2, data = df_rob)
mod.lmrob.r <- robustbase::lmrob(y ~ x1, data = df_rob)
mod.lmrob.int <- robustbase::lmrob(y ~ 1, data = df_rob)
set.seed(7)
n_grob <- 150
df_glmrob <- data.frame(
x1 = rnorm(n_grob),
x2 = factor(sample(c("a", "b", "c"), n_grob, replace = TRUE))
)
lp_grob <- 1.5 * df_glmrob$x1 +
ifelse(df_glmrob$x2 == "b", 1, -0.5)
df_glmrob$y <- rbinom(n_grob, 1, plogis(lp_grob))
mod.glmrob <- suppressWarnings(robustbase::glmrob(y ~ x1 + x2, family = binomial, data = df_glmrob))
mod.glmrob.r <- suppressWarnings(robustbase::glmrob(y ~ x1, family = binomial, data = df_glmrob))
mod.glmrob.int <- suppressWarnings(robustbase::glmrob(y ~ 1, family = binomial, data = df_glmrob))
test_that("lmrob: resi_pe produces correct estimates", {
pe <- suppressWarnings(resi_pe(mod.lmrob))
expect_equal(unname(pe$estimates),
c(2.1317677, 0.0321117, 1.7228418, 0.5869813, -0.2935287,
1.7298297, 0.9128559),
tolerance = 1e-06)
expect_equal(unname(pe$anova$RESI), c(1.7298297, 0.9128559), tolerance = 1e-06)
expect_equal(unname(pe$coefficients$RESI),
c(0.0321117, 1.7228418, 0.5869813, -0.2935287), tolerance = 1e-06)
})
test_that("lmrob: resi() returns RESI within CI bounds", {
set.seed(1)
ro <- suppressWarnings(resi(mod.lmrob, nboot = 100))
expect_true(all(ro$coefficients$RESI >= ro$coefficients$`2.5%`))
expect_true(all(ro$coefficients$RESI <= ro$coefficients$`97.5%`))
expect_true(all(ro$anova$RESI >= ro$anova$`2.5%`))
expect_true(all(ro$anova$RESI <= ro$anova$`97.5%`))
})
test_that("lmrob: specifying a reduced model only changes overall output", {
ro <- suppressWarnings(resi(mod.lmrob, nboot = 1))
ror <- suppressWarnings(resi(mod.lmrob, model.reduced = mod.lmrob.r, nboot = 1))
expect_equal(ro$coefficients$RESI, ror$coefficients$RESI)
expect_equal(ro$anova$RESI, ror$anova$RESI)
expect_false(ro$overall$RESI[2] == ror$overall$RESI[2])
})
test_that("lmrob: intercept-only model has no overall element", {
expect_true(is.null(suppressWarnings(resi(mod.lmrob.int, nboot = 1, anova = FALSE))$overall))
})
test_that("glmrob: resi_pe produces correct estimates", {
pe2 <- suppressWarnings(resi_pe(mod.glmrob))
expect_equal(unname(pe2$estimates),
c(0.4340252, -0.1257804, 0.4254354, 0.2549712, 0.0386988,
0.4175268, 0.2551174),
tolerance = 1e-06)
expect_equal(unname(pe2$anova$RESI), c(0.4175268, 0.2551174), tolerance = 1e-06)
})
test_that("glmrob: resi() returns RESI within CI bounds", {
set.seed(1)
ro2 <- suppressWarnings(resi(mod.glmrob, data = df_glmrob, nboot = 100))
expect_true(all(ro2$coefficients$RESI >= ro2$coefficients$`2.5%`))
expect_true(all(ro2$coefficients$RESI <= ro2$coefficients$`97.5%`))
expect_true(all(ro2$anova$RESI >= ro2$anova$`2.5%`))
expect_true(all(ro2$anova$RESI <= ro2$anova$`97.5%`))
})
test_that("glmrob: vcovHC redirects to stats::vcov with warning", {
expect_warning(
resi_pe(mod.glmrob, vcovfunc = sandwich::vcovHC),
"sandwich::vcovHC is not supported for glmrob"
)
})
test_that("glmrob: intercept-only model has no overall element", {
expect_true(is.null(suppressWarnings(resi(mod.glmrob.int, data = df_glmrob, nboot = 1, anova = FALSE))$overall))
})
}
Try the RESI package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.