Splines Work
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 <- 3 n_group <- 100 expected_obs <- 100 n_obs <- round(rexp(n_group, 1 / expected_obs)) group <- rep(1:n_group, n_obs) 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 = seq(-1, 2, length.out = p + 4)))[[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 = seq(-1, 2, length.out = p + 4)))[[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 + 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 <- matrix(c(1L, NA, NA, NA, NA, 2L, NA, NA, NA, NA, 3L, NA, NA, NA, NA, 4L), ncol = 4) 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 <- matrix(c(1L, NA, NA, NA, NA, 2L, 2L, 2L, NA, 2L, 2L, 2L, NA, 2L, 2L, 2L), ncol = 4) 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 model1$prior_cov model2$prior_cov model3$prior_cov
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.