## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", warning = FALSE, message = FALSE ) # library(MEMWAS) ## ----install-tarball, eval=FALSE---------------------------------------------- # install.packages( # "MEMWAS_0.9.3.tar.gz", # repos = NULL, # type = "source" # ) # #library(MEMWAS) ## ----install-github, eval=FALSE----------------------------------------------- # remotes::install_github("EnochKang/MEMWAS") # library(MEMWAS) ## ----capabilities, eval=FALSE------------------------------------------------- # MEMWAS_capabilities() ## ----fit-interface, eval=FALSE------------------------------------------------ # fit <- fit_MEMWAS( # formula, # family = "gaussian", # data, # id, # time, # random = ~ 1, # autocor = "AR(1)", # serial = NULL, # predictor_autocor = NULL, # autocor_predictors = NULL, # L1_penalty = 0, # L2_penalty = 0, # control = list(), # method = "ML", # random_cov = "unstructured", # approximation = "laplace", # init_approximation = "variational_inference", # quadrature_points = 7L, # se_method = "hessian", # link = NULL, # offset = NULL, # weights = NULL, # theta = NULL, # shape = NULL, # screen_nonlinear = FALSE, # check_assumptions = FALSE, # bootstrap_inference = FALSE, # prediction_inference = FALSE, # verbose = TRUE # ) ## ----multi-serial, eval=FALSE------------------------------------------------- # serial_spec <- list( # outcome = serial_component( # structure = "OU", # name = "outcome_persistence" # ), # exposure = serial_component( # structure = "expOU", # predictor = "exposure", # name = "exposure_persistence" # ) # ) # # fit_multi <- fit_MEMWAS( # outcome ~ time_c + exposure + treatment, # family = "gaussian", # data = dat, # id = "subject_id", # time = "visit_day", # random = ~(1 + time_c | subject_id), # autocor = NULL, # serial = serial_spec # ) ## ----diagnostics, eval=FALSE-------------------------------------------------- # diagnose_approximation(fit) ## ----education-example, eval=FALSE-------------------------------------------- # set.seed(2026) # # m <- 80L # n_time <- 4L # rho <- 0.45 # sigma_u <- 3 # sigma_eps <- 2 # sigma_b <- 5 # # student_dat <- data.frame( # student_id = rep(seq_len(m), each = n_time), # semester = rep(seq_len(n_time), times = m) # ) # # student_dat$study_minutes <- rnorm(m * n_time, 90, 25) # student_dat$attendance_rate <- runif(m * n_time, 0.70, 1.00) # b0 <- rnorm(m, 0, sigma_b) # # simulate_ar1 <- function(n, rho, sigma) { # u <- numeric(n) # u[1L] <- rnorm(1L, sd = sigma) # innovation_sd <- sigma * sqrt(1 - rho^2) # if (n > 1L) { # for (tt in 2L:n) { # u[tt] <- rho * u[tt - 1L] + rnorm(1L, sd = innovation_sd) # } # } # u # } # # u <- unlist( # lapply(seq_len(m), function(i) simulate_ar1(n_time, rho, sigma_u)), # use.names = FALSE # ) # # eps <- rnorm(m * n_time, 0, sigma_eps) # # student_dat$reading_score <- # 55 + # 0.08 * student_dat$study_minutes + # 18 * student_dat$attendance_rate + # b0[student_dat$student_id] + # u + eps # # fit_student <- fit_MEMWAS( # reading_score ~ study_minutes + attendance_rate, # family = "gaussian", # data = student_dat, # id = "student_id", # time = "semester", # random = ~ 1, # autocor = "AR(1)", # method = "REML", # approximation = "laplace", # verbose = FALSE # ) # # summary(fit_student) # diagnose_approximation(fit_student) ## ----healthcare-example, eval=FALSE------------------------------------------- # set.seed(2026) # # m <- 120L # visit_times <- c(0, 2, 5, 9, 14, 20) # n_time <- length(visit_times) # # clinic_dat <- data.frame( # patient_id = rep(seq_len(m), each = n_time), # visit_week = rep(visit_times, times = m) # ) # # clinic_dat$visit_week_c <- # clinic_dat$visit_week - mean(visit_times) # clinic_dat$therapy_sessions <- rpois(m * n_time, lambda = 2) # clinic_dat$sleep_hours <- rnorm(m * n_time, mean = 6.5, sd = 1.1) # # # Illustrative response only; a realistic simulation would also generate # # the stated random-slope and OU latent processes. # clinic_dat$depression_score <- # 14 - # 0.12 * clinic_dat$visit_week_c - # 0.30 * clinic_dat$therapy_sessions - # 0.60 * clinic_dat$sleep_hours + # rnorm(m * n_time, sd = 4) # # fit_clinic <- fit_MEMWAS( # depression_score ~ visit_week_c + therapy_sessions + sleep_hours, # family = "gaussian", # data = clinic_dat, # id = "patient_id", # time = "visit_week", # random = ~(1 + visit_week_c | patient_id), # random_cov = "unstructured", # autocor = "OU", # control = list(serial_time_scale = 1), # method = "REML", # verbose = FALSE # ) # # summary(fit_clinic) ## ----nongaussian-example, eval=FALSE------------------------------------------ # set.seed(2026) # # m <- 100L # n_time <- 4L # count_dat <- data.frame( # patient_id = rep(seq_len(m), each = n_time), # visit = rep(seq_len(n_time), times = m) # ) # # count_dat$visit_c <- count_dat$visit - mean(seq_len(n_time)) # count_dat$therapy_sessions <- rpois(m * n_time, 2) # count_dat$sleep_hours <- rnorm(m * n_time, 6.5, 1.1) # b0 <- rnorm(m, 0, 0.45) # # eta <- # -0.4 - # 0.08 * count_dat$visit_c + # 0.10 * count_dat$therapy_sessions - # 0.06 * count_dat$sleep_hours + # b0[count_dat$patient_id] # # count_dat$emergency_visits <- rnbinom( # nrow(count_dat), # mu = exp(eta), # size = 2.5 # ) # # fit_count <- fit_MEMWAS( # emergency_visits ~ visit_c + therapy_sessions + sleep_hours, # family = "negative_binomial", # data = count_dat, # id = "patient_id", # time = "visit", # random = ~ 1, # autocor = "NONE", # approximation = "laplace", # init_approximation = "variational_inference", # se_method = "hessian", # verbose = FALSE # ) # # summary(fit_count) # diagnose_approximation(fit_count) # # approx_comparison <- compare_approximations( # approximations = c( # "laplace", # "adaptive_gaussian_quadrature", # "variational_inference", # "pql" # ), # formula = emergency_visits ~ visit_c + therapy_sessions + sleep_hours, # family = "negative_binomial", # data = count_dat, # id = "patient_id", # time = "visit", # random = ~ 1, # autocor = "NONE", # quadrature_points = 7L # ) # # print(approx_comparison) ## ----predictor-serial-example, eval=FALSE------------------------------------- # fit_exposure <- fit_MEMWAS( # outcome ~ time_c + exposure + treatment, # family = "gaussian", # data = dat, # id = "subject_id", # time = "visit_day", # random = ~(1 | subject_id), # autocor = NULL, # serial = list( # outcome_persistence = serial_component( # structure = "OU", # name = "outcome_persistence" # ), # exposure_effect = serial_component( # structure = "expOU", # predictor = "exposure", # name = "exposure_effect" # ) # ) # )