tests/testthat/test_splines.R

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

library(mgcv)
library(tidyverse)

rm(list = ls())
devtools::load_all()
set.seed(1)

plot_design_matrix <- function(x, mat) {
    plot(0, xlim = range(x), ylim = range(mat), type = "n")
    for (i in seq_len(ncol(mat)))
        lines(x, mat[, i], col = i, lty = i)
    abline(h = 0, col = rgb(0, 0, 0, 0.2))
}

p <- 10
n_group <- 100
expected_obs <- 100
n_obs <- round(rexp(n_group, 1 / expected_obs))
group <- rep(1:n_group, n_obs)
knots <- seq(0, 1, length.out = p)
kd <- diff(knots)[1]
knots <- c(-2 * kd, -1 * kd, knots, 1 + 1 * kd, 1 + 2 * kd)

xx <- seq(0, 1, 0.01)

sc_xx <- smoothCon(s(xx, k = p + 1, bs = "ps", m = 1), data = data.frame(xx = xx),
                   absorb.cons = TRUE, knots = list(xx = knots))[[1]]
XX <- cbind(1, sc_xx$X)

plot_design_matrix(xx, sc_xx$X)
round(sc_xx$S[[1]], 4)
Sginv <- solve(sc_xx$S[[1]])

cov_weight <- Matrix::bdiag(1, Sginv)
covariance <- Matrix::bdiag(1, 2 * Sginv)
gamma <- chol(covariance)

slopes <- t(gamma) %*% matrix(rnorm(n_group * (p + 1)), ncol = n_group)

splines <- XX %*% slopes
plot(0, xlim = range(xx), ylim = range(splines), type = "n", xlab = "xx")
for (i in sample(n_group, 10))
    lines(xx, splines[, i], col = rgb(runif(1), runif(1), runif(1), 0.3))

x <- runif(length(group))
sc_x <- smoothCon(s(x, k = p + 1, bs = "ps", m = 1), data = data.frame(x = x),
                   absorb.cons = TRUE, knots = list(x = knots))[[1]]
X <- cbind(1, sc_x$X)
z <- list()
for (i in 1:n_group)
    z[[i]] <- X[group == i, , drop = FALSE] %*% slopes[, i, drop = FALSE]
z <- unlist(z)
noise <- rnorm(length(group))

y <- z + 2 * noise

i <- sample(n_group, 1)
plot(x[group == i], y[group == i])
lines(xx, XX %*% slopes[, i], col = 2, lwd = 2)

# test models -------------------------------------------------------------

eta_p <- 0.1
model1 <- gammmbest_fit(x = X, y = y, group = group, family = gaussian(),
                        fixef_index = 1:ncol(X), ranef_index = 1:ncol(X),
                        eta_precision = eta_p)

cov_supp <- diag(seq_len(p + 1))
cov_supp <- ifelse(cov_supp == 0, NA, cov_supp)
model2 <- gammmbest_fit(x = X, y = y, group = group, family = gaussian(),
                        fixef_index = 1:ncol(X), ranef_index = 1:ncol(X),
                        cov_supp = cov_supp, eta_precision = eta_p)

cov_supp <- as.matrix(Matrix::bdiag(1L, matrix(2L, p, p)))
cov_supp <- ifelse(cov_supp == 0, NA, cov_supp)

model3 <- gammmbest_fit(x = X, y = y, group = group, family = gaussian(),
                        fixef_index = 1:ncol(X), ranef_index = 1:ncol(X),
                        cov_supp = cov_supp, cov_weight = cov_weight,
                        eta_precision = eta_p)

# plot(slopes[1, ], model1$ranef[, 1][[1]]); abline(0, 1)
# plot(slopes[1, ], model2$ranef[, 1][[1]]); abline(0, 1)
# plot(slopes[1, ], model3$ranef[, 1][[1]]); abline(0, 1)
#
# cor(slopes[1, ], model2$ranef[, 1][[1]])
# cor(slopes[1, ], model3$ranef[, 1][[1]])
#
# plot(slopes[2, ], model1$ranef[, 2][[1]]); abline(0, 1)
# plot(slopes[2, ], model2$ranef[, 2][[1]]); abline(0, 1)
# plot(slopes[2, ], model3$ranef[, 2][[1]]); abline(0, 1)
#
# cor(slopes[2, ], model2$ranef[, 2][[1]])
# cor(slopes[2, ], model3$ranef[, 2][[1]])
#
# plot(slopes[3, ], model1$ranef[, 3][[1]]); abline(0, 1)
# plot(slopes[3, ], model2$ranef[, 3][[1]]); abline(0, 1)
# plot(slopes[3, ], model3$ranef[, 3][[1]]); abline(0, 1)
#
# cor(slopes[3, ], model2$ranef[, 3][[1]])
# cor(slopes[3, ], model3$ranef[, 3][[1]])
#
# plot(slopes[4, ], model1$ranef[, 4][[1]]); abline(0, 1)
# plot(slopes[4, ], model2$ranef[, 4][[1]]); abline(0, 1)
# plot(slopes[4, ], model3$ranef[, 4][[1]]); abline(0, 1)
#
# cor(slopes[4, ], model2$ranef[, 4][[1]])
# cor(slopes[4, ], model3$ranef[, 4][[1]])

model1$beta
model2$beta
model3$beta


Matrix::diag(model1$prior_cov) %>% round(2)
Matrix::diag(model2$prior_cov) %>% round(2)
Matrix::diag(model3$prior_cov) %>% round(2)
Matrix::diag(covariance) %>% round(2)