test_that("DEV: test linear regression works", {
set.seed(1)
N <- 100L
x <- rnorm(N)
beta <- rnorm(2)
y <- cbind(1, x) %*% beta + rnorm(N)
fit_lm <- lm(y ~ x)
fit_lbfgs <- fit_lbfgs(cbind(1, x), y, matrix(0, 2, 2), gaussian())
coef_glm_round6 <- round(coef(fit_lm), 6)
coef_lbfgs_round6 <- round(fit_lbfgs$par, 6)
expect_equivalent(coef_glm_round6, coef_lbfgs_round6)
})
test_that("DEV: test linear regression with offset works", {
set.seed(1)
N <- 100L
x <- rnorm(N)
z <- rnorm(N)
beta <- rnorm(2)
y <- cbind(1, x) %*% beta + rnorm(N) + 2 * z
fit_lm <- lm(y ~ x + offset(2 * z))
fit_lbfgs <- fit_lbfgs(cbind(1, x), y, matrix(0, 2, 2), gaussian(),
offset = 2 * z)
coef_glm_round6 <- round(coef(fit_lm), 6)
coef_lbfgs_round6 <- round(fit_lbfgs$par, 6)
expect_equivalent(coef_glm_round6, coef_lbfgs_round6)
})
test_that("DEV: test logistic regression works", {
set.seed(2)
N <- 100L
x <- rnorm(N)
beta <- rnorm(2)
p <- plogis(cbind(1, x) %*% beta)
y <- as.integer(runif(N) < p)
fit_glm <- glm(y ~ x , family = binomial())
fit_lbfgs <- fit_lbfgs(cbind(1, x), y, matrix(0, 2, 2), binomial())
coef_glm_round6 <- round(coef(fit_glm), 6)
coef_lbfgs_round6 <- round(fit_lbfgs$par, 6)
expect_equivalent(coef_glm_round6, coef_lbfgs_round6)
})
test_that("DEV: test logistic regression with offset works", {
set.seed(2)
N <- 100L
x <- rnorm(N)
z <- rnorm(N)
beta <- rnorm(2)
p <- plogis(cbind(1, x) %*% beta + z)
y <- as.integer(runif(N) < p)
fit_glm <- glm(y ~ x + offset(z), family = binomial())
fit_lbfgs <- fit_lbfgs(cbind(1, x), y, matrix(0, 2, 2), binomial(),
offset = z)
coef_glm_round6 <- round(coef(fit_glm), 6)
coef_lbfgs_round6 <- round(fit_lbfgs$par, 6)
expect_equivalent(coef_glm_round6, coef_lbfgs_round6)
})
test_that("DEV: ranef works with lme4", {
set.seed(1)
n_each <- 100L
n_group <- 100L
group <- rep(seq_len(n_group), n_each)
x <- rnorm(n_each * n_group)
fixef_int <- 1
fixef_slope <- -0.6
rand_int <- rnorm(n_group, sd = 0.2)
rand_slope <- rnorm(n_group, sd = 0.1)
noise <- rnorm(n_each * n_group)
y <- fixef_int + rand_int[group] +
(rand_slope[group] + fixef_slope) * x + noise
model <- lme4::lmer(y ~ 1 + x + (1 + x | group))
model_ranef <- as.matrix(lme4::ranef(model)[[1]])
x_old <- cbind(1, x)
x <- split.data.frame(cbind(1, x), group)
fixef_index <- 1:2
ranef_index <- 1:2
xz_svd <- Map(compact_svd, x, list(fixef_index), list(ranef_index))
y_old <- y
y <- split(y, group)
beta <- lme4::fixef(model)
prior_cov <- lme4::VarCorr(model)[[1]]
family <- gaussian()
alpha <- 0
fit <- gammmbest_fit(x = x_old, y = y_old, group = group, family = family,
fixef_index = 1:2, ranef_index = 1:2,
cov_supp = matrix(c(1, 2, 2, 3), 2, 2))
plot(unlist(fit$ranef), unlist(lme4::ranef(model)[[1]]))
abline(0, 1)
})
test_that("DEV: WRK WRK WRK WRK WRK", {
x <- cbind(1, lme4::sleepstudy$Days)
y <- lme4::sleepstudy$Reaction
group <- lme4::sleepstudy$Subject
family <- gaussian()
fixef_index <- 1:2
ranef_index <- 1:2
cov_supp <- matrix(c(1, 2, 2, 3), 2, 2)
fit <- gammmbest_fit(x, y, group, family, fixef_index, ranef_index,
cov_supp)
mbest_model <- mbest::mhglm(Reaction ~ Days + (Days | Subject),
data = lme4::sleepstudy)
plot(mbest::fixef.mhglm(mbest_model), fit$beta)
abline(0, 1)
plot(mbest::VarCorr.mhglm(mbest_model)[[1]], fit$prior_cov)
abline(0, 1)
plot(unlist(mbest::ranef.mhglm(mbest_model)[[1]]), unlist(fit$ranef))
abline(0, 1)
fit1 <- gammmbest_fit(x, y, group, family, fixef_index, ranef_index,
cov_supp, alpha = 1)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.