tests/testthat/test_dev.R

In kschmaus/gammmbest: Method Of Moments Based Estimation for Generalized Additive Mixed Models

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)
})