inst/doc/getting-started.R

In cevcmm: Communication-Efficient Varying Coefficient Mixed-Effects Models

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.width = 6,
  fig.height = 4,
  fig.align = "center"
)
set.seed(1)

## ----setup--------------------------------------------------------------------
library(cevcmm)

## ----simulate-----------------------------------------------------------------
set.seed(1)
N <- 500L
q <- 3L

t <- runif(N)
x <- runif(N)
Z <- matrix(rnorm(N * q), N, q)

beta_0     <- 2
beta_1_fun <- function(u) sin(2 * pi * u)
alpha_true <- rnorm(q, sd = 0.4)
sigma_eps  <- 0.5

y <- beta_0 + beta_1_fun(t) * x +
     as.vector(Z %*% alpha_true) +
     rnorm(N, sd = sigma_eps)

## ----fit----------------------------------------------------------------------
fit <- vcmm(y, X = x, Z = Z, t = t,
            method  = "auto",
            re_cov  = "diag",
            control = vcmm_control(sigma_eps       = 0.5,
                                   sigma_alpha     = 0.4,
                                   update_variance = TRUE))
fit

## ----inspect------------------------------------------------------------------
# Fixed-effects coefficient vector (intercept + spline basis coefs)
head(coef(fit))

# Same vector, reshaped: intercept + (basis x covariate) matrix
fx <- fixef(fit)
fx$intercept
dim(fx$varying)

# Random effects
ranef(fit)

# Sample size and residual SD
nobs(fit)
fit$sigma_eps

# Log-likelihood, AIC, BIC
logLik(fit)
AIC(fit); BIC(fit)

## ----summary------------------------------------------------------------------
summary(fit)

## ----curve, fig.cap = "Estimated and true varying coefficient."---------------
t_grid <- seq(0, 1, length.out = 100L)
vc     <- varying_coef(fit, t_new = t_grid, k = 1L, se.fit = TRUE)

plot(t_grid, vc$fit, type = "l", lwd = 2, col = "steelblue",
     xlab = "t", ylab = expression(hat(beta)[1](t)),
     ylim = range(vc$fit - 2 * vc$se.fit,
                  vc$fit + 2 * vc$se.fit,
                  beta_1_fun(t_grid)))
polygon(c(t_grid, rev(t_grid)),
        c(vc$fit + 2 * vc$se.fit, rev(vc$fit - 2 * vc$se.fit)),
        col = adjustcolor("steelblue", alpha.f = 0.25), border = NA)
lines(t_grid, beta_1_fun(t_grid), col = "red", lty = 2, lwd = 2)
legend("topright", c("estimate", "truth"),
       col = c("steelblue", "red"), lty = c(1, 2), lwd = 2, bty = "n")

## ----predict------------------------------------------------------------------
new_idx <- sample.int(N, 5L)
newdata <- list(t = t[new_idx],
                X = x[new_idx],
                Z = Z[new_idx, , drop = FALSE])

predict(fit, newdata = newdata)
y[new_idx]

## ----predict-se---------------------------------------------------------------
pred <- predict(fit, newdata = newdata, se.fit = TRUE)
cbind(fit = pred$fit, se = pred$se.fit)

## ----plot, fig.cap = "Built-in plot.vcmm_fit panel 1."------------------------
plot(fit, which = 1)