Nothing
if (isTRUE(as.logical(Sys.getenv("CI")))){
# If on CI
env_test <- "CI"
}else if (!identical(Sys.getenv("NOT_CRAN"), "true")){
# If on CRAN
env_test <- "CRAN"
set.seed(140) # CRAN SEED
}else{
# If on local machine
env_test <- 'local'
}
test_that("Test Robust for Linear (Unpenalized)", {
N <- 100
x <- rnorm(N)
y <- rnorm(N, sin(x)) * 100
est_linear_gam <- gam(y ~ x)
est_lm <- lm(y ~ x)
# Check fitted models are equivalent
expect_equivalent(sigma(est_lm), sigma(est_linear_gam))
expect_equivalent(fitted(est_lm), fitted(est_linear_gam))
expect_equivalent(vcov(est_lm), vcov(est_linear_gam))
hc_lm <- vcovHC(est_lm, type = 'HC0')
hc_linear_gam <- vcovHC(est_linear_gam, type = 'HC0')
ef_gam <- estfun(est_linear_gam, override_check = TRUE)
X_gam <- predict(est_linear_gam, newdata = data.frame(x = x), type = 'lpmatrix')
expect_equivalent(ef_gam, as.matrix(Diagonal(x = residuals(est_lm)) %*% X_gam))
expect_equivalent(ef_gam, estfun(est_lm))
tXXinv <- solve(crossprod(X_gam))
direct_estimation <- tXXinv %*%
t(X_gam) %*% Diagonal(x = residuals(est_linear_gam)^2) %*% X_gam %*%
tXXinv
expect_equivalent(as.matrix(direct_estimation), hc_lm)
expect_equivalent(hc_linear_gam, hc_lm)
expect_equivalent(vcovHC(est_lm, type = 'HC1'),
vcovHC(est_linear_gam, type = 'HC1'))
est_bam_linear <- bam(y ~ x)
expect_equivalent(vcov(est_bam_linear), vcov(est_lm))
expect_equivalent(vcovHC(est_bam_linear, type = 'HC1'),
vcovHC(est_lm, type = 'HC1'))
expect_equivalent(colSums(ef_gam), rep(0, ncol(ef_gam)))
})
test_that("Test Robust for Linear (Penalized)", {
N <- 1000
x <- rnorm(N)
z <- rnorm(N)
y <- rnorm(N, sin(x * z)) * 2
est_gam <- gam(y ~ te(x,z))
hc_gam <- vcovHC(est_gam, type = 'HC0')
ef_gam <- estfun(est_gam, override_check = TRUE)
X_gam <- predict(est_gam, newdata = data.frame(x = x, z = z), type = 'lpmatrix')
expect_equivalent(ef_gam, as.matrix(Diagonal(x = residuals(est_gam)) %*% X_gam))
expect_equivalent(ef_gam, estfun(est_gam, override_check = T))
S <- est_gam$smooth[[1]]$S[[2]] * est_gam$sp[2] + est_gam$smooth[[1]]$S[[1]] * est_gam$sp[1]
S <- bdiag(0, S)
m <- solve(crossprod(X_gam) + S, t(X_gam))
direct_estimation <- m %*% Diagonal(x = residuals(est_gam)^2) %*% t(m)
# Check direct formula for HC0
expect_equivalent(as.matrix(direct_estimation), hc_gam, tol = 1e-5)
# Check that HC1 works for "naive" correction
manual_hc1 <- vcovHC(est_gam, type = 'HC0') * N/(N-sum(est_gam$edf))
expect_equivalent(manual_hc1,
vcovHC(est_gam, type = 'HC1'))
expect_equivalent(as.vector(colSums(ef_gam) - S %*% coef(est_gam)),
rep(0, ncol(S)), tol = 1e-5)
})
test_that("Test Robust for GLM", {
N <- 100
x <- rnorm(N)
z <- rnorm(N)
# Check for poisson
y <- rpois(N, exp(sin(x)))
est_gam <- gam(y ~ s(x,z, bs = 'gKRLS'), family = poisson())
hc_gam <- vcovHC(est_gam, type = 'HC0')
ef_gam <- estfun(est_gam, override_check = TRUE)
X_gam <- predict(est_gam, newdata = data.frame(x = x, z = z), type = 'lpmatrix')
expect_equivalent(ef_gam, as.matrix(Diagonal(x = est_gam$y - fitted(est_gam)) %*% X_gam))
expect_equivalent(ef_gam, estfun(est_gam, override_check = T))
S <- est_gam$smooth[[1]]$S[[1]] * est_gam$sp[1]
S <- bdiag(0, S)
tXXinv <- solve(t(X_gam) %*% Diagonal(x = fitted(est_gam)) %*% X_gam + S)
direct_estimation <- tXXinv %*%
t(X_gam) %*% Diagonal(x = (y - fitted(est_gam))^2 ) %*% X_gam %*%
tXXinv
# Check direct formula for HC0
expect_equivalent(as.matrix(direct_estimation), hc_gam)
# Check that HC1 works for "naive" correction
manual_hc1 <- vcovHC(est_gam, type = 'HC0') * N/(N-sum(est_gam$edf))
expect_equivalent(manual_hc1,
vcovHC(est_gam, type = 'HC1'))
# Check for binomial with NONSTANDARD LINK
N <- 1000
x <- rnorm(N)
z <- rnorm(N)
y <- rbinom(N, 1, plogis( sqrt(abs(x)) + cos(z)))
est_gam <- gam(y ~ s(x) + s(z), family = binomial(link = 'cloglog'))
hc_gam <- vcovHC(est_gam, type = 'HC0')
ef_gam <- estfun(est_gam, override_check = TRUE)
X_gam <- predict(est_gam, newdata = data.frame(x = x, z = z), type = 'lpmatrix')
lp_gam <- X_gam %*% coef(est_gam)
score_weight <- y * exp(-exp(lp_gam) + lp_gam)/(1 - exp(-exp(lp_gam))) +
(1-y) * -exp(lp_gam)
hessian_weight <- exp(lp_gam) * (-1 - exp(exp(lp_gam)) * (-1 + exp(lp_gam)))/(-1 + exp(exp(lp_gam)))^2
e_y <- 1 - exp(-exp(lp_gam))
if (any(e_y %in% c(0,1))){
hessian_weight[which(e_y %in% c(0,1))] <- 0
}
hessian_weight <- as.vector(e_y * hessian_weight + (1-e_y) * -exp(lp_gam))
expect_equivalent(
score_weight, residuals(est_gam, 'pearson') * sqrt(weights(est_gam, 'working')),
tol = 1e-5
)
expect_equivalent(ef_gam, as.matrix(Diagonal(x = score_weight) %*% X_gam))
expect_equivalent(ef_gam, estfun(est_gam, override_check = T))
if (length(est_gam$smooth) == 0){
S <- Diagonal(x = rep(0, ncol(ef_gam)))
}else{
S <- bdiag(est_gam$smooth[[1]]$S[[1]] * est_gam$sp[1] , est_gam$smooth[[2]]$S[[1]] * est_gam$sp[2] )
S <- bdiag(0, S)
}
tXXinv <- solve(t(X_gam) %*% Diagonal(x = -hessian_weight) %*% X_gam + S)
expect_equivalent(as.matrix(tXXinv), vcov(est_gam))
direct_estimation <- tXXinv %*%
t(X_gam) %*% Diagonal(x = score_weight^2) %*% X_gam %*%
tXXinv
# Check direct formula for HC0
expect_equivalent(as.matrix(direct_estimation), hc_gam)
# Check that HC1 works for "naive" correction
manual_hc1 <- vcovHC(est_gam, type = 'HC0') * N/(N-sum(est_gam$edf))
expect_equivalent(manual_hc1, vcovHC(est_gam, type = 'HC1'))
expect_equivalent(as.vector(colSums(ef_gam) - S %*% coef(est_gam)), rep(0, ncol(S)),
tol = 1e-6)
})
test_that("Test Robust for bam", {
N <- 100
x <- rnorm(N)
z <- rnorm(N)
y <- rbinom(N, 1, plogis( sqrt(abs(x)) + cos(z)))
est_gam <- gam(y ~ x + z, family = binomial(link = 'cloglog'), method = 'REML')
est_bam <- bam(y ~ x + z, family = binomial(link = 'cloglog'), method = 'REML')
expect_equivalent(
vcovHC(est_gam, type = 'HC1'),
vcovHC(est_gam, type = 'HC0') * N/(N-sum(est_gam$edf))
)
est_gam <- gam(y ~ x + z, family = gaussian())
est_bam <- bam(y ~ x + z, family = gaussian())
c1 <- sample(1:15, size = N, replace = T)
expect_equivalent(
vcovCL(est_gam, cluster = c1),
vcovCL(est_bam, cluster = c1),
tol = 1e-5)
est_bam <- bam(y ~ te(x,z), family = gaussian())
mm <- predict(est_bam, newdata = data.frame(x, z), type = 'lpmatrix')
S <- est_bam$smooth[[1]]$S[[1]] * est_bam$sp[1] + est_bam$smooth[[1]]$S[[2]] * est_bam$sp[2]
meat <- t(mm) %*% Diagonal(x = residuals(est_bam)^2) %*% mm
inv <- solve(crossprod(mm) + bdiag(0, S))
# expect_equivalent(as.matrix(inv %*% meat %*% inv), vcovHC(est_bam, type = 'HC0'), scale = 1, tol = 1e-6)
expect_equivalent(as.matrix(inv %*% meat %*% inv), vcovHC(est_bam, type = 'HC0'), tol = 1e-5)
})
test_that("Test Robust for Complex Family", {
N <- 100
x <- rnorm(N)
z <- rnorm(N)
y <- rt(N,df = 2) + exp(x) + cos(z)
est_gam <- gam(y ~ x + z, family = scat, method = 'REML')
expect_error(vcovHC(est_gam), regexp = 'Robust SE from sandwich not set up')
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.