Nothing
context("Stepwise")
set.seed(1)
n <- 30
p <- 10
X <- matrix(rnorm(p*n), ncol = p)
colnames(X) <- letters[1:p]
y <- X[, 2] + 2*X[, 3] - X[, 6] + rnorm(n)
d <- prepare_data(y, X)
Xm <- cbind(1, X[, c(2, 3, 6)])
loglik <- loglik(y, Xm)
bic_v <- round(bic(loglik, k = 4, n = n), 2)
test_that("Typical data", {
res <- stepwise(d, crit = bic)
expect_equal(res$model, c("c", "f", "b"))
d$Xm <- cbind(1, X)
expect_equal(stepwise(d, crit = bic)$model, c("b", "c", "f"))
d$Xm <- cbind(1, X[, 1:5])
expect_equal(stepwise(d, crit = bic)$model, c("b", "c", "f"))
d$Xm <- cbind(1, X[, c(1, 4, 9)])
d$candidates <- (1:10)[-3]
expect_equal(stepwise(d, crit = bic)$model, c("b", "j"))
# is this model better than without "j"?
lm1 <- lm(y ~ X[, 2] + X[, 10])
lm2 <- lm(y ~ X[, 2])
expect_true(stats::BIC(lm1) < stats::BIC(lm2))
# is this model better than with "f" instead of "j?
lm1 <- lm(y ~ X[, 2] + X[, 10])
lm2 <- lm(y ~ X[, 2] + X[, 6])
expect_true(stats::BIC(lm1) < stats::BIC(lm2)) # strange, but it is TRUE
expect_equal(stepwise(d, crit = aic)$model, c("b", "j", "f")) # aic sees "f
d$candidates <- 10:1
d$Xm <- matrix(1, nrow = n)
expect_equal(stepwise(d, crit = bic)$model, c("c", "f", "b"))
})
test_that("Rare situations", {
d <- prepare_data(y, X)
d$candidates <- NULL
stepwise(d, crit = bic)
expect_message(stepwise(d, crit = bic), "Done.")
expect_equal(stepwise(d, crit = bic)$model, NULL)
})
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.