## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----sim-sofr, eval=F---------------------------------------------------------
# set.seed(123)
# n <- 100
# M <- 50
# tgrid <- seq(0, 1, length.out = M)
# dt    <- tgrid[2] - tgrid[1]
# tmat  <- matrix(rep(tgrid, each = n), nrow = n)
# lmat  <- matrix(dt, nrow = n, ncol = M)
# 
# # Smooth latent trajectory + measurement noise on the functional predictor
# D_true <- t(apply(matrix(rnorm(n * M), n, M), 1, cumsum)) / sqrt(M)
# wmat   <- D_true + matrix(rnorm(n * M, sd = 0.3), nrow = n) ## noisy
# 
# beta_true <- sin(2 * pi * tgrid)
# X1        <- rnorm(n)
# eta       <- 0.5 * X1 + D_true %*% (beta_true * dt)   ## regression on D, not W
# prob      <- plogis(eta)
# y         <- rbinom(n, 1, prob)
# 
# data.SoFR <- data.frame(y = y, X1 = X1)
# data.SoFR$tmat <- tmat
# data.SoFR$lmat <- lmat
# data.SoFR$wmat <- wmat

## ----fit-sofr, eval=F---------------------------------------------------------
# library(refundBayes)
# 
# fit_sofr_joint <- sofr_bayes(
#   formula    = y ~ X1 + s(tmat, by = lmat * wmat, bs = "cc", k = 10),
#   data       = data.SoFR,
#   family     = binomial(),
#   joint_FPCA = c(TRUE),       ## << turn joint FPCA on for the (only) functional term
#   niter      = 1500,
#   nwarmup    = 500,
#   nchain     = 3,
#   ncores     = 3
# )

## ----fit-sofr-nojoint, eval=F-------------------------------------------------
# fit_sofr_plain <- sofr_bayes(
#   formula = y ~ X1 + s(tmat, by = lmat * wmat, bs = "cc", k = 10),
#   data    = data.SoFR,
#   family  = binomial(),
#   niter   = 1500, nwarmup = 500, nchain = 3, ncores = 3
# )

## ----compare-sofr, eval=F-----------------------------------------------------
# plot(fit_sofr_joint)
# plot(fit_sofr_plain)

## ----fit-fcox, eval=F---------------------------------------------------------
# library(refundBayes)
# 
# ## Use the example dataset shipped with the FCox vignette
# Func_Cox_Data <- readRDS("data/example_data_Cox.rds")
# Func_Cox_Data$wmat <- Func_Cox_Data$MIMS
# Func_Cox_Data$cens <- 1 - Func_Cox_Data$event
# 
# fit_cox_joint <- fcox_bayes(
#   formula    = survtime ~ X1 + s(tmat, by = lmat * wmat, bs = "cc", k = 10),
#   data       = Func_Cox_Data,
#   cens       = Func_Cox_Data$cens,
#   joint_FPCA = c(TRUE),
#   niter      = 5000,
#   nwarmup    = 2000,
#   nchain     = 1,
#   ncores     = 1
# )

## ----inspect-fcox, eval=F-----------------------------------------------------
# fcox_code <- fcox_bayes(
#   formula    = survtime ~ X1 + s(tmat, by = lmat * wmat, bs = "cc", k = 10),
#   data       = Func_Cox_Data,
#   cens       = Func_Cox_Data$cens,
#   joint_FPCA = c(TRUE),
#   runStan    = FALSE
# )
# cat(fcox_code$stancode)

## ----sim-fofr, eval=F---------------------------------------------------------
# library(refundBayes)
# set.seed(42)
# 
# n <- 200
# L <- 30
# M <- 30
# sindex <- seq(0, 1, length.out = L)
# tindex <- seq(0, 1, length.out = M)
# 
# # Smooth latent functional predictor + measurement noise
# D_true <- matrix(0, nrow = n, ncol = L)
# for (i in 1:n) {
#   D_true[i, ] <- rnorm(1) * sin(2 * pi * sindex) +
#                  rnorm(1) * cos(2 * pi * sindex) +
#                  rnorm(1) * sin(4 * pi * sindex)
# }
# X_func <- D_true + matrix(rnorm(n * L, sd = 0.3), nrow = n)   ## noisy
# 
# age <- rnorm(n)
# 
# # True coefficient functions
# beta_true  <- outer(sin(2 * pi * sindex), cos(2 * pi * tindex))
# alpha_true <- 0.5 * sin(pi * tindex)
# 
# # Generate response from the latent D_true (not from X_func!)
# signal_scalar <- outer(age, alpha_true)
# signal_func   <- (D_true %*% beta_true) / L
# epsilon       <- matrix(rnorm(n * M, sd = 0.3), nrow = n)
# Y_mat         <- signal_scalar + signal_func + epsilon
# 
# dat <- data.frame(age = age)
# dat$Y_mat  <- Y_mat
# dat$X_func <- X_func
# dat$sindex <- matrix(rep(sindex, n), nrow = n, byrow = TRUE)

## ----fit-fofr, eval=F---------------------------------------------------------
# fit_fofr_joint <- fofr_bayes(
#   formula     = Y_mat ~ age + s(sindex, by = X_func, bs = "cr", k = 10),
#   data        = dat,
#   joint_FPCA  = c(TRUE),
#   spline_type = "bs",
#   spline_df   = 10,
#   niter       = 2000,
#   nwarmup     = 1000,
#   nchain      = 3,
#   ncores      = 3
# )

## ----plot-fofr, eval=F--------------------------------------------------------
# beta_est <- apply(fit_fofr_joint$bivar_func_coef[[1]], c(2, 3), mean)
# 
# par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
# image(sindex, tindex, beta_true,
#       xlab = "s (predictor domain)", ylab = "t (response domain)",
#       main = expression("True " * beta(s,t)),
#       col = hcl.colors(64, "Blue-Red 3"))
# image(sindex, tindex, beta_est,
#       xlab = "s (predictor domain)", ylab = "t (response domain)",
#       main = expression("Joint-FPCA " * hat(beta)(s,t)),
#       col = hcl.colors(64, "Blue-Red 3"))