## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(lmmprobe)
library(MASS) # For mvrnorm if needed

## ----simulate-----------------------------------------------------------------
set.seed(123)

n <- 50 # Number of subjects
m <- 4 # Observations per subject
p <- 20 # Number of high-dimensional predictors (Z)
sigma_e <- 1.5 # Residual standard deviation
sigma_b <- 1.0 # Random intercept standard deviation

# Subject IDs
ID <- rep(1:n, each = m)
N <- n * m

# High-dimensional predictors (Z) - simple random normal
Z <- matrix(rnorm(N * p), ncol = p)
colnames(Z) <- paste0("z", 1:p)

# Sparse coefficients for Z
beta_true <- rep(0, p)
beta_true[1:3] <- c(0.5, -0.5, 0.3) # First 3 are active

# Random Intercepts
b <- rnorm(n, 0, sigma_b)
b_vec <- rep(b, each = m)

# Response variable Y
# Y = Z*beta + b + e
Y <- Z %*% beta_true + b_vec + rnorm(N, 0, sigma_e)

# Fixed effect covariates (V)
# For scenario 1 (random intercept), V is a 1-column matrix of IDs.
V <- matrix(ID, nrow = N, ncol = 1)

## ----lmmprobe_run-------------------------------------------------------------
# Preparing inputs
# Y is vector/matrix
# Z is matrix of high-dim predictors
# ID_data is vector of IDs
# V is matrix of IDs (for random intercept)

fit <- lmmprobe(
    Y = matrix(Y, ncol = 1),
    Z = Z,
    V = V,
    ID_data = ID,
    Y_test = matrix(Y, ncol = 1),
    Z_test = Z,
    V_test = V,
    ID_test = ID,
    alpha = 0.05,
    maxit = 10, # Short run for vignette
    ep = 0.5,
    cpus = 1, # Serial for vignette
    B = 3,
    adj = 1,
    LR = 0,
    C = 1,
    sigma_init = NULL
)


# Inspect results
print(fit$c_coefs)
print(head(fit$beta))