---
title: "Population-level marginal summaries from smooth two-slope mixed models with random changepoints"
author: "Divan A. Burger, Sean van der Merwe, Emmanuel Lesaffre"
format: 
  revealjs: 
    theme: [default, UFS1.scss]
    logo: UFSlogo2022SvdM.svg
    footer: "2026/08/14 - Two slope"
    slide-number: "c/t"
    chalkboard: true
    parallax-background-image: "pastel_curves_light.jpg"
    parallax-background-size: "2976px 1663px"
    parallax-background-vertical: 50
    margin: 0.05
execute:
  echo: false
bibliography: Bibliography.bib
---


```{r}
#| warning: false
#| message: false
#| results: hide
options(scipen = 12)
library(tidyverse)
theme_set(theme_bw())
library(runjags)
library(coda)
library(loo)
library(DHARMa)
```

```{r}
Actg315 <- ushr::actg315raw
# Data transformation
Actg315$CD4 <- Actg315$CD4/100
Actg315$RNA <- 10^(Actg315$log10.RNA)
Actg315$Threshold <- 1 - (Actg315$RNA <= 100)
Actg315$Log10RNAcens <- ifelse(Actg315$RNA <= 100, 2, Actg315$log10.RNA)
Actg315$Month <- Actg315$Day/30
Actg315$PatientID <- Actg315$Patid |> factor()
Actg315_long <- Actg315 |>
  pivot_longer(cols = c(log10.RNA., CD4), 
               names_to = "Measure", 
               values_to = "Value")
```

```{r data2}
Actg315_long_filtered <- Actg315_long[Actg315_long$Measure == "log10.RNA.", ]

RNAplot <- Actg315_long_filtered |> ggplot(aes(x = Day, y = Value, 
                                    group = PatientID, colour = PatientID)) +
  geom_line(linewidth = 0.5) +
  scale_colour_viridis_d(option = "B") +
  labs(x = "Day", y = expression(paste(log[10], " RNA Copies/mL"))) +
  theme(legend.position = "none") +
  geom_hline(yintercept = 2, colour = "purple", linetype = 2) + 
  scale_x_continuous(breaks = seq(0, 1000, by = 14))
```

```{r}
# Data list for JAGS
DataList <- list(
  N = length(unique(Actg315$Patid)),
  Ntotal = nrow(Actg315),
  Id = match(Actg315$Patid, sort(unique(Actg315$Patid))),
  Month = Actg315$Month,
  Y = Actg315$Log10RNAcens,
  Censor = Actg315$Threshold,
  Upper = rep(2, nrow(Actg315)),
  Ones = rep(1, nrow(Actg315)),
  C = 1e6,
  TauMax = 180/30,
  pi = pi
)


BilinearParams <- c(
  "Sigma_resid",
  "Beta_RE_bar", "L_chol",
  "Alpha0", "Beta1", "Beta2",
  "Tau", "Mu_tau", "Sigma_tau",
  "Lambda", "Nu"
)
```

```{r}
NumCores <- ceiling(parallel::detectCores(logical = FALSE)*0.8)

BilinearCacheFile <- "C:/temp/Simulations/Actg315_BilinearJags.rds"
Cached <- readRDS(BilinearCacheFile)
BilinearMcmc <- Cached$BilinearMcmc 
SumBilinearMcmc <- Cached$SumBilinearMcmc

SoftHingeCacheFile <- "C:/temp/Simulations/Actg315_SoftHingeJags.rds"
Cached <- readRDS(SoftHingeCacheFile)
SoftHingeMcmc <- Cached$SoftHingeMcmc
SumSoftHingeMcmc <- Cached$SumSoftHingeMcmc
rm(Cached)
```


```{r}
# E[(t - tau)+ | tau ~ TN(mu, sigma^2; L, R)]
BilinearITrunc <- function(t, mu, sigma, L, R) {
  CoreFun <- function(t, mu, sigma, L, R) {
    u <- (L - mu)/sigma
    v <- (R - mu)/sigma
    k <- pnorm(v) - pnorm(u)
    muT <- mu + sigma*(dnorm(u) - dnorm(v))/k
    if (t <= L) return(rep(0, length(mu)))
    if (t >= R) return(t - muT)
    z <- (t - mu)/sigma
    (t - mu)*(pnorm(z) - pnorm(u))/k +
      sigma*(dnorm(z) - dnorm(u))/k
  }
  if (length(t) == 1) {
    return(CoreFun(t, mu, sigma, L, R))
  }
  out <- numeric(length(t))
  idxL <- t <= L
  idxR <- t >= R
  idxM <- !(idxL | idxR)
  if (any(idxL)) out[idxL] <- 0
  if (any(idxR)) {
    u <- (L - mu[idxR])/sigma[idxR]
    v <- (R - mu[idxR])/sigma[idxR]
    k <- pnorm(v) - pnorm(u)
    muT <- mu[idxR] + sigma[idxR]*(dnorm(u) - dnorm(v))/k
    out[idxR] <- t[idxR] - muT
  }
  if (any(idxM)) {
    u <- (L - mu[idxM])/sigma[idxM]
    v <- (R - mu[idxM])/sigma[idxM]
    k <- pnorm(v) - pnorm(u)
    z <- (t[idxM] - mu[idxM])/sigma[idxM]
    out[idxM] <- (t[idxM] - mu[idxM])*(pnorm(z) - pnorm(u))/k +
      sigma[idxM]*(dnorm(z) - dnorm(u))/k
  }
  out
}

# G*(t) = d/dt E[(t - tau)+] under truncated normal
BilinearITruncPrime <- function(Time, Mu, Sigma, L, R) {
  U <- (L - Mu)/Sigma
  V <- (R - Mu)/Sigma
  K <- pnorm(V) - pnorm(U)
  Z <- (Time - Mu)/Sigma
  if (length(Time) == 1) {
    Out <- numeric(length(Mu))
    if (Time <= L) {
      Out[] <- 0
    } else if (Time >= R) {
      Out[] <- 1
    } else {
      Out <- (pnorm(Z) - pnorm(U))/K
    }
    return(Out)
  }
  Out <- numeric(length(Time))
  Out[Time <= L] <- 0
  Out[Time >= R] <- 1
  Mid <- Time > L & Time < R
  Out[Mid] <- (pnorm(Z[Mid]) - pnorm(U))/K
  Out
}

# Posterior mean of tau under TN(mu, sigma^2; L, R)
BilinearComputeTNMeanTau <- function(MuTau, SigmaTau, L, R) {
  MuTau + SigmaTau*
    (dnorm((L - MuTau)/SigmaTau) - dnorm((R - MuTau)/SigmaTau))/
    (pnorm((R - MuTau)/SigmaTau) - pnorm((L - MuTau)/SigmaTau))
}
```

```{r}
# Posterior summary bands
SummarizeBand <- function(Mat) {
  cbind(
    Mean = apply(Mat, 2, mean),
    Lower = apply(Mat, 2, quantile, 0.025),
    Upper = apply(Mat, 2, quantile, 0.975)
  )
}

# Extract posterior draws
BilinearExtract <- function(BilinearMcmc, DataList) {
  ChainsList <- BilinearMcmc
  Chains <- do.call(rbind, lapply(ChainsList, as.data.frame))
  list(
    Chains = Chains,
    Alpha0 = Chains[["Beta_RE_bar[1]"]],
    Beta1 = Chains[["Beta_RE_bar[2]"]],
    Beta2 = Chains[["Beta_RE_bar[3]"]],
    MuTau = Chains[["Mu_tau"]],
    SigmaTau = Chains[["Sigma_tau"]],
    Ltau = 0,
    Rtau = DataList$TauMax
  )
}

# Generic plotter for mean or slope
# type = "mean": plot mu_marg(t) vs mu_cond(t)
# type = "slope": plot d mu_marg/dt vs d mu_cond/dt
PlotBands <- function(Bands, type = c("mean", "slope")) {
  type <- match.arg(type)
  if (type == "mean") {
    PlotData <- data.frame(
      Time = Bands$Tgrid,
      MeanMarg = Bands$BandMargMean[, "Mean"],
      LowerMarg = Bands$BandMargMean[, "Lower"],
      UpperMarg = Bands$BandMargMean[, "Upper"],
      MeanPlug = Bands$BandPlugMean[, "Mean"],
      LowerPlug = Bands$BandPlugMean[, "Lower"],
      UpperPlug = Bands$BandPlugMean[, "Upper"]
    )
    ylabTxt <- expression("Mean "*log[10]*" RNA copies/mL")
    LegendBreaks <- c("Conditional", "Marginal")
    LegendLabels <- c(
      "Conditional" = expression(mu[cond*","*1](italic(t))),
      "Marginal" = expression(mu[marg*","*1](italic(t)))
    )
  } else {
    PlotData <- data.frame(
      Time = Bands$Tgrid,
      MeanMarg = Bands$BandMargSlope[, "Mean"],
      LowerMarg = Bands$BandMargSlope[, "Lower"],
      UpperMarg = Bands$BandMargSlope[, "Upper"],
      MeanPlug = Bands$BandPlugSlope[, "Mean"],
      LowerPlug = Bands$BandPlugSlope[, "Lower"],
      UpperPlug = Bands$BandPlugSlope[, "Upper"]
    )
    ylabTxt <- "Mean rate of change (log10 RNA copies/mL/mo)"
  }
  
  ggplotly(
    PlotData |> ggplot(aes(x = Time)) + 
      # Marginal band
      geom_ribbon(
        aes(
          ymin = LowerMarg,
          ymax = UpperMarg,
          fill = "Marginal",
          color = "Marginal"
        ),
        alpha = 0.35,
        linewidth = 0.7,
        linetype = "dashed"
      ) +
      geom_line(
        aes(y = MeanMarg, color = "Marginal"),
        linewidth = 1.1,
        linetype = "dashed"
      ) +
      # Conditional band
      geom_ribbon(
        aes(
          ymin = LowerPlug,
          ymax = UpperPlug,
          fill = "Plug-in",
          color = "Plug-in"
        ),
        alpha = 0.25,
        linewidth = 0.7,
        linetype = "solid"
      ) +
      geom_line(
        aes(y = MeanPlug, color = "Plug-in"),
        linewidth = 1.1,
        linetype = "solid"
      ) +
      xlim(0, 1) + 
      labs(
        x = "Months since baseline",
        y = ylabTxt
      ) + 
    guides(color = guide_legend(title = ""), 
           fill = guide_legend(title = "Line"))
    , width = 950, height = 400
  )
}

Bilinear <- list()
Bilinear$Extract <- BilinearExtract(
  BilinearMcmc,
  DataList
)
Bilinear$Bands <- readRDS("BilinearBandsShort.rds")
# Mean trajectories (marginal vs conditional)
# Bilinear$FigMean <- PlotBands(
#   Bands = Bilinear$Bands,
#   type = "mean"
# )
```

```{r}
# Core expectation under truncated normal tau, given kappa
# Returns:
#   Eh_vec(t) = E_tau[h_kappa(t - tau) | kappa]
#   EPhi_vec(t) = E_tau[Phi(kappa*(t - tau)) | kappa]
# where tau ~ TN(mu_tau, sigma_tau^2; L_tau, R_tau)
SoftHingeSlopeKernelGivenKappaTN <- function(t, mu_tau, sigma_tau, L_tau, R_tau, kappa) {
  t <- as.numeric(t)
  zL <- (L_tau - mu_tau)/sigma_tau
  zU <- (R_tau - mu_tau)/sigma_tau
  Ztr <- pnorm(zU) - pnorm(zL)
  m <- t - mu_tau
  a <- kappa*sigma_tau
  b <- kappa*m
  d <- sqrt(1 + a^2)
  I1 <- pbivnorm::pbivnorm(zU, b/d, rho = a/d) -
    pbivnorm::pbivnorm(zL, b/d, rho = a/d)
  I1/Ztr
}

CumTrapzZero <- function(X, Y) {
  X <- as.numeric(X)
  Y <- as.numeric(Y)
  Out <- numeric(length(X))
  if (length(X) <= 1L) {
    return(Out)
  }
  Out[-1] <- cumsum(0.5*diff(X)*(Y[-1] + Y[-length(Y)]))
  Out
}
```


```{r}
SoftHingeExtract <- function(SoftHingeMcmc, DataList) {
  ChainsList <- SoftHingeMcmc
  Chains <- do.call(rbind, lapply(ChainsList, as.data.frame))
  list(
    Chains = Chains,
    Alpha0 = Chains[["Beta_RE_bar[1]"]],
    Beta1 = Chains[["Beta_RE_bar[2]"]],
    Beta2 = Chains[["Beta_RE_bar[3]"]],
    MuTau = Chains[["Mu_tau"]],
    SigmaTau = Chains[["Sigma_tau"]],
    MuKappa = Chains[["Mu_kappa"]],
    SigmaKappa = Chains[["Sigma_kappa"]],
    Ltau = 0,
    Rtau = DataList$TauMax
  )
}

Soft <- list()
Soft$Extract <- SoftHingeExtract(
  SoftHingeMcmc,
  DataList
)
Soft$Bands <- readRDS("SoftHingeBands.rds")
```


```{r}
BuildIndex <- function(ParamNames, Nsubj, HasKappa = TRUE) {
  Out <- list(
    SigmaResid = "Sigma_resid" |> match(ParamNames),
    Lambda = "Lambda" |> match(ParamNames),
    Nu = "Nu" |> match(ParamNames),
    Alpha0 = paste0("Alpha0[", 1:Nsubj, "]") |> match(ParamNames),
    Beta1 = paste0("Beta1[", 1:Nsubj, "]") |> match(ParamNames),
    Beta2 = paste0("Beta2[", 1:Nsubj, "]") |> match(ParamNames),
    Tau = paste0("Tau[", 1:Nsubj, "]") |> match(ParamNames)
  )
  if (HasKappa) {
    Out$Kappa <- paste0("Kappa[", 1:Nsubj, "]") |> match(ParamNames)
  }
  Out
}

ExtractPosteriorFull <- function(McmcObj, Needed) {
  McmcList <- McmcObj
  ChainMats <- lapply(McmcList, function(Mc) as.matrix(Mc))
  CommonCols <- Reduce(intersect, lapply(ChainMats, colnames))
  KeepCols <- intersect(Needed, CommonCols)
  ChainMats <- lapply(
    ChainMats,
    function(M) as.matrix(M[, KeepCols, drop = FALSE])
  )
  Posterior <- do.call(rbind, ChainMats)
  NPerChain <- vapply(ChainMats, nrow, integer(1))
  ChainId <- rep(seq_along(ChainMats), times = NPerChain)
  list(
    Posterior = Posterior,
    ParamNames = colnames(Posterior),
    chain_id = ChainId
  )
}
```

```{r}
# Soft hinge helper function
SoftHingeCenteredHinge <- function(Time, Tau, Kappa, KappaTol = 1e-4) {
  DeltaT <- Time - Tau
  Delta0 <- -Tau
  Centered <- 0*Time
  UseApprox <- Kappa < KappaTol
  if (any(!UseApprox)) {
    DeltaTExact <- DeltaT[!UseApprox]
    Delta0Exact <- Delta0[!UseApprox]
    KappaExact <- Kappa[!UseApprox]
    Centered[!UseApprox] <- DeltaTExact*pnorm(KappaExact*DeltaTExact) +
      dnorm(KappaExact*DeltaTExact)/KappaExact -
      (Delta0Exact*pnorm(KappaExact*Delta0Exact) +
         dnorm(KappaExact*Delta0Exact)/KappaExact)
  }
  if (any(UseApprox)) {
    Phi0 <- dnorm(0)
    Centered[UseApprox] <- 0.5*Time[UseApprox] +
      0.5*Phi0*Kappa[UseApprox]*
      (DeltaT[UseApprox]^2 - Delta0[UseApprox]^2)
  }
  Centered
}

PointwiseLogLik <- function(PosteriorSamples, Ix, DrawChunk = 200L,
                            Model = "SoftHinge") {
  S <- nrow(PosteriorSamples)
  N <- length(DataList$Y)
  seq(1L, S, by = DrawChunk) |> lapply(\(SFrom) {
    STo <- min(S, SFrom + DrawChunk - 1L)
    NSamp <- STo - SFrom + 1L
    Theta <- PosteriorSamples[SFrom:STo, , drop = FALSE]
    SigmaResid <- Theta[, Ix$SigmaResid]
    Lambda <- Theta[, Ix$Lambda]
    Nu <- Theta[, Ix$Nu]
    Alpha0 <- Theta[, Ix$Alpha0, drop = FALSE]
    B1 <- Theta[, Ix$Beta1, drop = FALSE]
    B2 <- Theta[, Ix$Beta2, drop = FALSE]
    Tau <- Theta[, Ix$Tau, drop = FALSE]
    LogC <- lgamma((Nu + 1)/2) -
      0.5*log(pi*(Nu - 2)) -
      lgamma(Nu/2)
    Cc <- exp(LogC)
    Aconst <- 4*Lambda*Cc*(Nu - 2)/(Nu - 1)
    B2c <- 1 + 3*Lambda^2 - Aconst^2
    Bconst <- sqrt(B2c)
    Kfac <- sqrt(Nu/(Nu - 2))
    Knot <- -Aconst/Bconst
    Alpha0Block <- Alpha0[,DataList$Id]
    B1Block <- B1[,DataList$Id]
    B2Block <- B2[,DataList$Id]
    TauBlock <- Tau[,DataList$Id]
    TimeMat <- matrix(DataList$Month, nrow = NSamp, ncol = N, byrow = TRUE)
    if (Model == "SoftHinge") {
      Kappa <- Theta[, Ix$Kappa, drop = FALSE]
      H <- SoftHingeCenteredHinge(
        Time = TimeMat,
        Tau = TauBlock,
        Kappa = Kappa[,DataList$Id]
      )
    } else {
      H <- pmax(0, TimeMat - TauBlock)
    }
    Mu <- Alpha0Block + B1Block*TimeMat + B2Block*H
    Z <- (matrix(DataList$Y, NSamp, N, byrow = TRUE) - Mu)/SigmaResid
    Zup <- (matrix(DataList$Upper, NSamp, N, byrow = TRUE) - Mu)/SigmaResid
    Zt <- Z*Bconst + Aconst
    ZtUp <- Zup*Bconst + Aconst
    Left <- (Knot - Z) >= 0
    LeftUp <- (Knot - Zup) >= 0
    BaseL <- 1 + (Zt/(1 - Lambda))^2/(Nu - 2)
    BaseR <- 1 + (Zt/(1 + Lambda))^2/(Nu - 2)
    PartL <- Left*Bconst*Cc*BaseL^(-(Nu + 1)/2)
    PartR <- (!Left)*Bconst*Cc*BaseR^(-(Nu + 1)/2)
    Pdf <- (PartL + PartR)/SigmaResid
    tLUp <- (ZtUp/(1 - Lambda))*Kfac
    tRUp <- (-ZtUp/(1 + Lambda))*Kfac
    CdfLeft <- (1 - Lambda)*pt(tLUp, df = Nu)
    CdfRight <- 1 - (1 + Lambda)*pt(tRUp, df = Nu)
    Cdf <- ifelse(LeftUp, CdfLeft, CdfRight)
    CensorMat <- matrix(DataList$Censor, NSamp, N, byrow = TRUE)
    CensorMat*log(Pdf) + (1 - CensorMat)*log(Cdf)
  }) |> do.call(rbind, args = _)
}
```

```{r}
NSubj <- DataList$N

NeededBilin <- c(
  "Sigma_resid", "Lambda", "Nu",
  paste0("Alpha0[", 1:NSubj, "]"),
  paste0("Beta1[", 1:NSubj, "]"),
  paste0("Beta2[", 1:NSubj, "]"),
  paste0("Tau[", 1:NSubj, "]")
)
NeededSoft <- c(
  NeededBilin,
  paste0("Kappa[", 1:NSubj, "]")
)
ExtSoft <- ExtractPosteriorFull(SoftHingeMcmc, NeededSoft)
if (!file.exists("loo_calculations.rds")) {

  ExtBilin <- ExtractPosteriorFull(BilinearMcmc, NeededBilin)
  IxBilin <- BuildIndex(ExtBilin$ParamNames, NSubj, HasKappa = FALSE)
  LogLikMatBilin <- PointwiseLogLik(
    PosteriorSamples = ExtBilin$Posterior,
    Ix = IxBilin, 
    Model = "Bilinear"
  )
  REffBilin <- loo::relative_eff(exp(LogLikMatBilin), 
                                 chain_id = ExtBilin$chain_id)
  LooBilinConditional <- loo::loo(LogLikMatBilin, 
                                  r_eff = REffBilin, 
                                  moment_match = TRUE)

  IxSoft <- BuildIndex(ExtSoft$ParamNames, NSubj, HasKappa = TRUE)
  LogLikMatSoft <- PointwiseLogLik(
    PosteriorSamples = ExtSoft$Posterior, 
    Ix = IxSoft, 
    Model = "SoftHinge"
  )
  REffSoft <- loo::relative_eff(exp(LogLikMatSoft), 
                                chain_id = ExtSoft$chain_id)
  LooSoftConditional <- loo::loo(LogLikMatSoft, 
                                 r_eff = REffSoft, 
                                 moment_match = TRUE)
  loo_results <- loo::loo_compare(list(SoftHinge = LooSoftConditional, 
                 Bilinear = LooBilinConditional))
  loo_results |> saveRDS("loo_calculations.rds")
} else {
  loo_results <- readRDS("loo_calculations.rds")
}

```


```{r}
ParamToLatex <- function(ParamName) {
  Lookup <- c(
    "Sigma_resid" = " $\\sigma_{1}$",
    "Lambda" = " $\\Lambda_{1}$",
    "Nu" = " $\\nu_{1}$",
    "MeanTauTn" = "$\\bar{\\tau}_{1}$",
    "SdTauTn" = "$\\mathrm{SD}\\left(\\tilde{\\tau}_{1}\\right)$",
    "MeanKappaTn" = "$\\bar{\\kappa}_{1}$",
    "SdKappaTn" = "$\\mathrm{SD}\\left(\\tilde{\\kappa}_{1}\\right)$",
    "Beta_RE_bar\\[1\\]" = "$\\alpha_{0,1}$",
    "Beta_RE_bar[2]" = "  $\\beta_{1,1}$",
    "Beta_RE_bar[3]" = "  $\\beta_{2,1}$"
  )
  p_l <- Lookup[match(ParamName, names(Lookup))] |> unname()
  p_l[is.na(p_l)] <- ParamName[is.na(p_l)]
  sub(
    "^L_chol\\[([0-9]+),([0-9]+)\\]$",
    "$L_{1,\\1\\2}$",
    p_l
  )
}

ExtractSummariesActg <- function(McmcList, ModelName) {
  AllPars <- coda::varnames(McmcList)
  DropPattern <- "^(Alpha0\\[|Beta1\\[|Beta2\\[|Tau\\[|Kappa\\[)"
  KeepPars <- AllPars[!grepl(DropPattern, AllPars)]
  AllChains <- do.call(rbind, lapply(McmcList, \(Ch) as.matrix(Ch))
                       )[, KeepPars, drop = FALSE]
  GelOut <- lapply(McmcList, \(Ch) Ch[, KeepPars, drop = FALSE]) |>
    coda::gelman.diag(autoburnin = FALSE, multivariate = FALSE)
  data.frame(
    Model = ModelName,
    ParameterName = KeepPars |> ParamToLatex(), 
    Median = apply(AllChains, 2, median),
    Lower95 = apply(AllChains, 2, quantile, 0.025),
    Upper95 = apply(AllChains, 2, quantile, 0.975),
    PSRFPointEst = GelOut$psrf[, 1],
    PSRFUpperCI = GelOut$psrf[, 2],
    row.names = NULL
  ) |> filter(is.finite(PSRFPointEst))
}

AllSummaries <- rbind(
  ExtractSummariesActg(BilinearMcmc, "Abrupt"), 
  ExtractSummariesActg(SoftHingeMcmc, "Smooth")
  ) |> arrange(ParameterName)

```


# Introduction

In general, 

$$E[f(X)]\ne f(E[X])$$

and 

$$\hat{f}(\boldsymbol\theta)\ne f(\hat{\boldsymbol\theta})$$

yet statisticians keep assuming equality.

## Source and authors

This presentation is drawn from our paper that is under review at the Journal of Biopharmaceutical Statistics.

Divan A. Burger$^{1,2}$, Sean van der Merwe$^{2}$, Emmanuel Lesaffre$^{3,4}$

1.    Cytel Inc., Waltham, MA, USA
1.    Department of Mathematical Statistics and Actuarial Science, University of the Free State, Bloemfontein, South Africa
1.    I-BioStat, KU Leuven, Leuven, Belgium
1.    Department of Statistics and Actuarial Science, University of Stellenbosch, Stellenbosch, South Africa

```{r}
#| warning: false
#| message: false
#| results: hide
library(knitr)
opts_chunk$set(fig.ext = 'svg', dev = 'svg', fig.align = 'center')

UFScolours <- c(Blue="#0F204BFF", Red="#A71930FF", LightGray="#A7A8AAFF", DarkGray="#8D8D8EFF", NAS="#0039A7FF", EDU="#00675AFF", THEO="#C69317FF", HUM="#EA8400FF", SA="#A40084FF", SC="#9E83B7FF", QWA="#00B140FF", LAW="#BB133EFF", HLTH="#490E6FFF", Black="#000000FF", VeryLight="#D7D8DDFF")
cols <- UFScolours[c(8,11,10,2)]

wrap_strings <- function(vector_of_strings, width){ as.character( sapply( vector_of_strings, \(x) { paste( strwrap(x, width=width), collapse="\n")}))}
library(plotly)
plotly_table_UFS <- function(datafrm, height = "100%", width = "100%", colors=UFScolours,
    useRowNames = FALSE, rowNamesHeading = " ", highlights = NULL, 
    header_font_size = 14, cell_font_size = 12) {
  # For the colour palette: the first two and last two colours are used
  datafrm <- as.data.frame(datafrm)
  nms <- names(datafrm)
  if (useRowNames) {
    datafrm <- data.frame(row.names(datafrm), datafrm, row.names = NULL)
    nms <- c(rowNamesHeading, nms)
    names(datafrm) <- nms
  }
  k <- ncol(datafrm)
  algn <- rep('center', k)
  cell_values <- rbind(t(as.matrix(unname(datafrm))))
  col_list <- rep(c("#FFFFFFFF", colors[length(colors)]), 
                  nrow(datafrm) %/% 2 + 2)[1:nrow(datafrm)] |> rep(k)
  if (!is.null(highlights)) {
    col_list[highlights] <- "#FFBBBBFF"
  }
  col_list <- unname(col_list)
  col_list <- col_list |> matrix(nrow(cell_values), ncol(cell_values), byrow = TRUE)
  fig <- plot_ly(
  type = 'table',
  height = height,
  width = width, 
  header = list(
    values = nms,
    align = algn,
    line = list(width = 1, color = colors[2]),
    fill = list(color = colors[1]),
    font = list(family = "Arial", size = header_font_size, color = "#FFFFFFFF")
  ),
  cells = list(
    values = cell_values,
    align = algn,
    line = list(color = colors[2], width = 1),
    fill = list(color = col_list),
    font = list(family = "Arial", size = cell_font_size, color = colors[length(colors)-1])
  ))
fig
}
plotly_colours <- UFScolours |> unname()
```

## Outline

-   Motivation
-   Dataset to be used as example
-   Modelling aspects
-   More modelling aspects
-   Conclusions

## Motivation

::: {.incremental}
-   We want people to get better from illness
-   Understand the process of getting better
-   Understand the factors that affect the process
-   Understand the statistical uncertainty and variation
    -   Better modelling of the variation $\longrightarrow$ better modelling of the process
:::

## Spoiler

::: {.callout-caution}
# Consider the rate at which people get better from HIV when taking ARVs

The expected rate of healing for a random future person is not the same as for the *typical person* (doesn't exist)
:::

<!-- 
There's an old joke in statistics that says, "The average person has one testicle." 
--->

```{r}
#| warning: false
# Instantaneous slope (marginal vs conditional)
PlotBands(
  Bands = Bilinear$Bands,
  type = "slope"
)
```


# Data set

The ACTG 315 dataset, available in the *ushr* **R** package [@MORRIS2020], includes longitudinal measurements of HIV viral load (log$_{10}$ RNA copies/mL) over time. It features data on 46 patients, with the longest measurement recorded on Day 196 after baseline (Day 0). 

## Data set {.smaller}

```{css}
table.dataTable tbody td {
  font-size: 16px;
  padding: 6px 10px;
}
```

```{r}
actg315 <- ushr::actg315raw
display_data <- actg315 |> select(-Patid)
display_data$PatientID <- actg315$Patid |> factor()
library(DT)
datatable(display_data, rownames = FALSE, filter = "top", 
          options = list(pageLength = 10, dom = "tp"))
```

## RNA by time and subject

```{r}
RNAplot |> ggplotly(width = 950, height = 600)
```

## Censoring

We observe two kinds of censoring:

1.  Subjects stop taking part in the study at different lengths of time
1.  A proportion of observations censored below the lower detection limit

> We address these in different ways:

1.  By storing the data in long form and fitting a mixed effects model (each patient has their own curve)
1.  By putting censored observations into the likelihood via their CDF (not PDF), *i.e.* $P[X <= LDL]$ instead of $P[X = x]$; *e.g.* @lachos2011linear

# The modelling

We build up the model step by step.

## Mixed effects models

-   Mixed effects models include fixed effects and random effects (e.g. random intercepts or random slopes per subject)
-   Mixed effects models are used when we have sampling of nominal observation groups from a population of possible groups
-   Typically a random sample of people from a population of people, but with **multiple observations per person**
    -   Can also be a random sample of animals, random sample of fields, random sample of classes or schools

Random effects are important to help address the fact that our residuals are **not independent**.

## Non-linear regression {.smaller}

-   Looking at the log RNA plot we see a non-linear pattern. 
    -   First a fast curve down 
    -   then a slow curve down/flat/up
    -   with a smooth transition

::: {.callout-note}
# Each person has their own curve

Each person has an intercept, an initial slope, a turning point, an ultimate slope, and their own turning smoothness. Some turn gradually while others are more kinky (literally, not figuratively).
:::

Examples of papers implementing this idea in various fields include: @perelson1996hiv, @heerspink2021effect, @bacon1971estimating, @donald2008early, @burger2015bayesian, @hall2000change, @ghosh2007random, @hout2011smooth, *etc.*

## Non-linear curve illustration {.smaller}

While we did formally use the DHARMa approach in the paper to check the goodness-of-fit for the two models [@HARTIG2021B], I want to encourage all researchers to try to **plot the model and data together** if at all possible (sometimes some creativity and some conditioning is needed).

```{r}
Bands = Soft$Bands
PlotData <- data.frame(
  Time = Bands$Tgrid,
  MeanMarg = Bands$BandMargMean[, "Mean"],
  LowerMarg = Bands$BandMargMean[, "Lower"],
  UpperMarg = Bands$BandMargMean[, "Upper"]
)
ylabTxt <- expression("Mean log10 RNA copies/mL")
model_fit_plot <- Actg315_long_filtered |> ggplot() +
  geom_line(aes(x = Month, y = Value, group = PatientID, 
                colour = PatientID), 
            linewidth = 0.5, alpha = 0.4) +
  scale_colour_viridis_d(option = "B") +
  labs(x = "Months since baseline", y = ylabTxt) +
  theme(legend.position = "none") +
  geom_hline(yintercept = 2, colour = "purple", linetype = 2) + 
  geom_ribbon(aes(x = Time, ymin = LowerMarg, ymax = UpperMarg),
      data = PlotData, fill = "steelblue", colour = "blue",
    alpha = 0.4, linewidth = 1, linetype = "dashed"
  ) +
  geom_line(
    aes(x = Time, y = MeanMarg), data = PlotData, 
    linewidth = 1.2, linetype = "dashed", color = "darkblue"
  )
```

```{r}
model_fit_plot |> ggplotly(height = 500, width = 900)
```



## Model fit procedure

-   The tool we used for fitting this particular model is [JAGS](https://mcmc-jags.sourceforge.io/)
    -   JAGS stands for Just Another Gibbs Sampler
-   The fitting was done using the [R](https://cran.r-project.org/bin/windows/base/) package [runjags](https://cran.r-project.org/web/packages/runjags/index.html) [@DENWOOD2016]
-   JAGS uses clever Gibbs sampling to arrive at posterior simulations for all the parameters
-   Once you have posterior simulations you can calculate anything you want
    -   Fits and predictions for a particular subject
    -   Fits for the average subject
    -   Predictions for a random future subject, with full uncertainty

## Model code {.scrollable}

```{r}
#| echo: true
#| eval: false
model {
  for (r in 1:3) {
    L_chol[r, r] ~ dt(0, 0.5, 3)T(0, )
    for (c in 1:(r - 1)) {
      L_chol[r, c] ~ dnorm(0, 0.001)
    }
    for (c in (r + 1):3) {
      L_chol[r, c] <- 0
    }
  }

  for (r in 1:3) {
    for (c in 1:3) {
      Prec_RE[r, c] <- inprod(L_chol[r, 1:3], L_chol[c, 1:3])
    }
  }

  for (j in 1:3) {
    Beta_RE_bar[j] ~ dnorm(0, 0.001)
  }

  Sigma_resid ~ dt(0, 0.5, 3)T(0, )
  Tau_resid <- pow(Sigma_resid, -2)

  Lambda ~ dunif(-1, 1)
  Eps ~ dexp(1)
  NuMinus2 ~ dgamma(2, Eps)
  Nu <- NuMinus2 + 2

  LogC <- loggam((Nu + 1)/2) - 0.5*log(pi*(Nu - 2)) - loggam(Nu/2)
  Cc <- exp(LogC)
  Aconst <- 4*Lambda*Cc*(Nu - 2)/(Nu - 1)
  B2 <- 1 + 3*pow(Lambda, 2) - pow(Aconst, 2)
  Bconst <- sqrt(B2)
  Kfac <- sqrt(Nu/(Nu - 2))
  Knot <- -Aconst/Bconst

  Mu_tau ~ dnorm(0, 0.001)
  Sigma_tau ~ dt(0, 0.5, 3)T(0, )
  Tau_tau <- pow(Sigma_tau, -2)

  Mu_kappa ~ dnorm(0, 0.001)
  Sigma_kappa ~ dt(0, 0.5, 3)T(0, )
  Tau_kappa <- pow(Sigma_kappa, -2)

  Phi0 <- 1/sqrt(2*pi)
  KappaSwitch <- 1.0E-4

  for (i in 1:N) {
    v[i, 1:3] ~ dmnorm(Beta_RE_bar[1:3], Prec_RE[,])
    Alpha0[i] <- v[i, 1]
    Beta1[i] <- v[i, 2]
    Beta2[i] <- v[i, 3]

    Tau[i] ~ dnorm(Mu_tau, Tau_tau)T(0, TauMax)
    Kappa[i] ~ dnorm(Mu_kappa, Tau_kappa)T(0, )
  }

  for (n in 1:Ntotal) {
    DeltaT[n] <- Month[n] - Tau[Id[n]]
    Delta0[n] <- -Tau[Id[n]]
    
    Zk[n] <- Kappa[Id[n]]*DeltaT[n]
    Zk0[n] <- Kappa[Id[n]]*Delta0[n]

    PhiZ[n] <- pnorm(Zk[n], 0, 1)
    PhiZ0[n] <- pnorm(Zk0[n], 0, 1)
    phiZ[n] <- dnorm(Zk[n], 0, 1)
    phiZ0[n] <- dnorm(Zk0[n], 0, 1)

    HingeExact[n] <- DeltaT[n]*PhiZ[n] + phiZ[n]/Kappa[Id[n]]
    Hinge0Exact[n] <- Delta0[n]*PhiZ0[n] + phiZ0[n]/Kappa[Id[n]]
    CenteredHingeExact[n] <- HingeExact[n] - Hinge0Exact[n]

    CenteredHingeApprox[n] <- 0.5*Month[n] +
      0.5*Phi0*Kappa[Id[n]]*(pow(DeltaT[n], 2) - pow(Delta0[n], 2))

    UseApprox[n] <- step(KappaSwitch - Kappa[Id[n]])
    CenteredHinge[n] <- UseApprox[n]*CenteredHingeApprox[n] +
      (1 - UseApprox[n])*CenteredHingeExact[n]

    Mu[n] <- Alpha0[Id[n]] + Beta1[Id[n]]*Month[n] +
      Beta2[Id[n]]*CenteredHinge[n]

    Z[n] <- (Y[n] - Mu[n])/Sigma_resid
    Zup[n] <- (Upper[n] - Mu[n])/Sigma_resid
    Ztilde[n] <- Bconst*Z[n] + Aconst
    ZtildeUp[n] <- Bconst*Zup[n] + Aconst

    Left[n] <- step(Knot - Z[n])
    Right[n] <- 1 - Left[n]
    BaseL[n] <- 1 + pow(Ztilde[n]/(1 - Lambda), 2)/(Nu - 2)
    BaseR[n] <- 1 + pow(Ztilde[n]/(1 + Lambda), 2)/(Nu - 2)
    PartL[n] <- Left[n]*(Bconst*Cc)*pow(BaseL[n], -(Nu + 1)/2)
    PartR[n] <- Right[n]*(Bconst*Cc)*pow(BaseR[n], -(Nu + 1)/2)
    PDFstd[n] <- PartL[n] + PartR[n]
    PDF[n] <- PDFstd[n]/Sigma_resid

    LeftUp[n] <- step(Knot - Zup[n])
    RightUp[n] <- 1 - LeftUp[n]
    tL_up[n] <- (ZtildeUp[n]/(1 - Lambda))*Kfac
    tR_up[n] <- (-ZtildeUp[n]/(1 + Lambda))*Kfac
    CDF_left_up[n] <- (1 - Lambda)*pt(tL_up[n], 0, 1, Nu)
    CDF_right_up[n] <- 1 - (1 + Lambda)*pt(tR_up[n], 0, 1, Nu)
    CDFstd_up[n] <- LeftUp[n]*CDF_left_up[n] + RightUp[n]*CDF_right_up[n]
    CDF[n] <- CDFstd_up[n]

    Like[n] <- pow(PDF[n], Censor[n])*pow(CDF[n], 1 - Censor[n])
    p[n] <- Like[n]/C
    Ones[n] ~ dbern(p[n])
  }
}
```

# The two curve options considered

We consider two equations for $\mu$:

-   The classic bi-linear model where two lines meet at a point
    -   Note that this can also end up looking like a smooth transition if we allow for parameter uncertainty, subject uncertainty, or both (the main point of the paper)
-   A new model form which introduces a transition from one line to the other with flexible smoothness

## Two-line (bi-phasic) curve illustration {.scrollable}

-   A lot of real-world problems involve two straight-line trajectories meeting at a changepoint
    -   HIV, some kidney issues, and TB are cited in the paper
    -   I have seen it in farming (herbicide, insecticide, etc.)
-   The big question is, "How kinky is the transition from one curve to another?"

```{r}
x <- 0:100
f1 <- pmax(10 - x*0.2, 0)
f2 <- 1 + x*0.06
w_f2 <- seq(1, 3.3, l = 101) |> log() |> pmin(1)
y_ave <- f2*w_f2 + f1*(1-w_f2)
fig_data <- data.frame(x = x, f1 = f1, f2 = f2, y_ave, w_f2) |> 
  pivot_longer(-1, names_to = "Line", values_to = "y")
(fig_data |> ggplot(aes(x = x, y = y, colour = Line, linetype = Line)) + 
  geom_line(linewidth = 1) + geom_hline(yintercept = 1)) |> 
  ggplotly(width = 900, height = 400)
```


## Bilinear curve (standard model)

-   Early examples of the model include @bacon1971estimating, and especially @laird1982random - who include random effects, *i.e.* each subject has their own curve
-   The model is two lines that meet at an intersection point (for each subject)
-   This type of model is also called a bi-phasic regression
-   Given an unknown changepoint $\gamma$, the **simplest form of the model** can be formally defined as 

$$y_i\sim N(\beta_0+\beta_1(x_i-\gamma)I(x_i>\gamma)+\beta_2(x_i-\gamma)I(x_i<\gamma),\ \sigma^2)$$

## New Soft Hinge curve {.scrollable}

-   We model individual trajectories, optionally stratified by treatment, using a smooth transition around a subject-specific changepoint. 
-   Let $j = 1,\ldots,J$ index treatment groups, and let $i = 1,\ldots,N_j$ index subjects assigned to treatment $j$, where $N_j$ is the number of subjects in group $j$. 
-   For each subject $i$, let $k = 1,\ldots,n_{ij}$ index repeated measurements, where $n_{ij}$ is the number of observations available for that subject, and let the corresponding observation times be $t_{ijk}$. 
-   We write $N = \sum_{j=1}^{J} \sum_{i=1}^{N_j} n_{ij}$ for the total number of observations across all treatment groups and subjects. 
-   We take $t = 0$ as baseline for all subjects, so that changepoint times are measured relative to baseline.
-   We use the term 'hinge' for the standard broken-stick basis function $x^{+} = \max\left(x,0\right)$, which is zero for $x \le 0$ and linear for $x > 0$, so that including $\beta_{2ij} x^{+}$ in a linear predictor creates a kink at $x = 0$. 

To obtain a differentiable version, we work with the smooth hinge

$$
      h_{\kappa}\left(x\right) =
      x \Phi\left(\kappa x\right) +
      \frac{\phi\left(\kappa x\right)}{\kappa},
      \qquad \kappa > 0,
$$

which is continuous and differentiable, with derivative $h_{\kappa}'\left(x\right) = \Phi\left(\kappa x\right)$, where $\Phi$ and $\phi$ denote the standard normal cumulative distribution function (CDF) and density, respectively. For large $\kappa$, $h_{\kappa}\left(x\right)$ is close to the hard hinge $x^{+} = \max\left(x,0\right)$; for small $\kappa$, the bend is gradual and spans a wider time interval around $x = 0$.

The observed outcome is $Y_{ijk}$, with conditional mean for $k = 1,\ldots,n_{ij}$:

$$
      \mu_{ijk} = \alpha_{0ij} + \beta_{1ij} t_{ijk} +
      \beta_{2ij} \left\{h_{\kappa_{ij}}\left(t_{ijk} - \tau_{ij}\right) - h_{\kappa_{ij}}\left(-\tau_{ij}\right)\right\},
$$
where $\tau_{ij}$ is a subject-specific changepoint and $\kappa_{ij} > 0$ controls the smoothness of the transition.

The instantaneous slope with respect to $t$ is

$$
      \frac{d \mu_{ij}\left(t\right)}{d t} = \beta_{1ij} +
      \beta_{2ij} \Phi\left(\kappa_{ij}\left(t - \tau_{ij}\right)\right),
$$
which transitions smoothly from $\beta_{1ij}$ for $t \ll \tau_{ij}$ to $\beta_{1ij} + \beta_{2ij}$ for $t \gg \tau_{ij}$.

The parameter $\beta_{1ij}$ is the initial slope before the gate, and $\beta_{1ij} + \beta_{2ij}$ is the long-term slope after the gate.

The subject-specific smoothness parameter $\kappa_{ij}$ controls how quickly the slope changes; a larger $\kappa_{ij}$ means a sharper bend, and a smaller $\kappa_{ij}$ means a slower bend over time.

Conditional on the mean trajectory, the outcome follows an error distribution:

$$
      Y_{ijk} \left| \mu_{ijk}, \boldsymbol{\theta}_j \right.
      \sim f_{\varepsilon}\left(y; \mu_{ijk}, \boldsymbol{\theta}_j\right),
$$

where $\boldsymbol{\theta}_j$ collects dispersion, skewness, and tail-thickness parameters for treatment group $j$, and the parameterization of $f_{\varepsilon}$ ensures $E\left(\varepsilon_{ijk}\right) = 0$ for $\varepsilon_{ijk} = Y_{ijk} - \mu_{ijk}$. Examples include Gaussian, Student-$t$, and skew-$t$ families, parameterized so that $E\left(\varepsilon_{ijk}\right) = 0$; allowing heavy-tailed residuals in nonlinear mixed-effects models is standard practice, *e.g.* @davidian1995nonlinear @davidian2003nonlinear. In the Bayesian implementation for ACTG~315 we take $f_{\varepsilon}$ to be a skew-$t$ density.


## Curve comparison {.smaller}

-   For model comparison between the two curve options, we use the **leave-one-out cross-validation information criterion (LOOIC)** [@vehtari2024LOO]
    -   It approximates how well the model might have predicted each observation had it been fitted without ever having seen that observation
    -   Not true out-of-sample but a close approximation that is useful for smaller datasets and *relatively* fast

```{r}
loo_df <- loo_results |> as.data.frame() |> round(1)
loo_df |> mutate(
  Model = rownames(loo_results), .before = 1
) |> select(-c(6,7,9)) |> 
  plotly_table_UFS(height = 300, width = 900)
```


# Posterior predictive distribution

To enable accurate calculation of predictive quantities we must be able to simulate from the posterior predictive distribution. 

A simulation function for the residual distribution is created, which can then be applied for each combination of posterior vectors and new data matrices to arrive as predictions for any quantity of interest.

In this study the residual distribution used comes from @hansen1994autoregressive (<https://www.jstor.org/stable/2527081>).

```{r}
# Skew-t density simulation function
rskewt <- function (n = 1, nu = 4, l = 0, mu = 0, s = 1) {
  # Input checks
  l <- pmin(pmax(l, -0.999), 0.999)
  nu <- pmax(nu, 2 + 1e-6)
  # Left or right of mode
  u <- runif(n)
  Pleft <- (1 - l)/2
  pos <- u > Pleft
  result <- numeric(n)
  # Constants
  LogC <- lgamma((nu)/2) - 0.5*log(pi*(nu - 2)) - lgamma(nu/2)
  Aconst <- 4*l*exp(LogC)*(nu - 2)/(nu - 1)
  Bconst <- sqrt(1 + 3*l^2 - Aconst^2)
  Kfac <- sqrt(nu/(nu - 2))
  # Simulation of skew t, left then right
  result[!pos] <- (1-l)*qt(u[!pos]/(1-l), nu)/Kfac
  result[pos] <- -(1+l)*qt((1-u[pos])/(1+l), nu)/Kfac
  # Scale and shift
  (result - Aconst)/Bconst * s + mu
}
```

## Non-symmetric conditional distributions {.smaller}

-   For symmetric distributions such as the *normal*, *logistic*, and *t* distributions:
    -   $\mu\ \equiv$ mean $\equiv$ median $\equiv$ mode
-   For the rest we have to explicitly state what we are trying to explain

When building models, ask yourself whether you are trying to **model or predict**, then ask yourself whether you are trying to model or predict

-   The expected value (mean)
-   The middle value (median)
-   The most likely value (mode)

> Consider an ordinary normal regression on the log scale. What happens when you transform back? Are you still modelling the mean?

## Single subject exploration

-   Here we consider a single subject and illustrate their fit with uncertainty. 
-   This could easily be transformed into a web app on a company portal to allow a practitioner to interact with the model fit in future cases.
-   First code is created to illustrate a generic subject. This function calculates and plots the expected trajectory, the 95% credibility intervals of the expected trajectory, and 95% prediction intervals for a given subject.

```{r}
shortestinterval <- function(postsims, width=0.95) {
  sort(postsims) -> sorted.postsims
  round(length(postsims)*width) -> gap
  which.min(diff(sorted.postsims, gap)) -> pos
  sorted.postsims[c(pos, pos + gap)]
}

plot_subject <- function(sbjid) {
  # The subject data is isolated from the data set
  sbj_data <- Actg315 |> filter(Patid %in% sbjid) |> 
    select(-c(Obs.No, Patid, PatientID))
  # The relevant model simulations are extracted
  #   Note: using a lot of simulations can slow performance 
  #         so a subset of simulations is used here
  v_select <- c(1:3, which(ExtSoft$ParamNames |> 
                             endsWith(paste0("[", sbjid, "]"))))
  nsims <- 5000
  which_sims <- sample(seq_len(nrow(ExtSoft$Posterior)), nsims)
  sims_mat <- ExtSoft$Posterior[which_sims, v_select]
  vnms <- ExtSoft$ParamNames[v_select] |> str_remove("\\[[0-9]+\\]$")
  sims <- sims_mat |> as.data.frame() |> setNames(vnms)
  # The distributions over time for the subject are calculated
  Tgrid <- seq(0, 6, by = 0.02)
  npnts <- length(Tgrid)
  sbj_mu_dist <- seq_len(nsims) |> sapply(\(i) {
    d <- sims[i,]
    DeltaT <- Tgrid - d$Tau
    Delta0 <- -d$Tau
    Centered <- DeltaT*pnorm(d$Kappa*DeltaT) +
      dnorm(d$Kappa*DeltaT)/d$Kappa -
      (Delta0*pnorm(d$Kappa*Delta0) +
         dnorm(d$Kappa*Delta0)/d$Kappa)
    mu <- d$Alpha0 + d$Beta1*Tgrid + d$Beta2*Centered
  })
  sbj_pred_dist <- seq_len(nsims) |> sapply(\(i) {
    d <- sims[i,]
    rskewt(npnts, d$Nu, d$Lambda, sbj_mu_dist[,i], d$Sigma_resid)
  })
  # Calculate estimates and intervals from distributions
  get_stats <- function(sims) {
    int <- shortestinterval(sims)
    c(Pred = median(sims), Lower = int[1], Upper = int[2])
  }
  sbj_mu_stats <- sbj_mu_dist |> apply(1, get_stats) |> t()
  sbj_stats <- sbj_pred_dist |> apply(1, get_stats) |> t()
  sbj_df <- data.frame(Month = Tgrid, sbj_stats)
  sbj_mu_df <- data.frame(Month = Tgrid, sbj_mu_stats)
  # Draw plot
  sbj_mu_df |> ggplot(aes(x = Month)) + 
  geom_ribbon(aes(ymin = Lower |> smooth() |> c(), 
                  ymax = Upper |> smooth() |> c()), 
              data = sbj_df, fill = "lightblue", 
              colour = "lightblue", alpha = 0.1) + 
  geom_ribbon(aes(ymin = Lower |> smooth() |> c(), 
                  ymax = Upper |> smooth() |> c()), 
              fill = "cyan", 
              colour = "cyan", alpha = 0.1) + 
  geom_line(aes(y = Pred), linewidth = 1, colour = "blue") + 
  geom_point(aes(y = log10.RNA.), data = sbj_data, colour = "red") + 
  labs(x = "Months since baseline", y = ylabTxt)
}
```

For illustration we consider a few subjects. 

## Individual curve illustration  {.smaller}

:::: {.columns}

::: {.column width="50%"}

Subject 1 had few points and then dropped out so produced wide intervals.

```{r}
plot_subject(1) |> ggplotly(width = 440, height = 220)
```

Subject 2 did not follow the model pattern so produced even wider intervals which barely accommodate the observed pattern.

```{r}
plot_subject(2) |> ggplotly(width = 440, height = 220)
```

:::

::: {.column width="50%"}

Subject 8 fitted the expected data pattern nearly perfectly, showing neat alignment and smaller intervals.

```{r}
plot_subject(8) |> ggplotly(width = 440, height = 220)
```

Subject 9 is extra, pun intended 😉 

```{r}
plot_subject(9) |> ggplotly(width = 440, height = 220)
```

:::

::::

# Conclusion

The real data curves are smooth, not kinky, so we smooth out the kink in 3 ways ...

-   Bayesian parameter simulation
    -   Even with the simplest two-line model, you get a smooth transition if you 'integrate' over the uncertainty in the turning point
-   Allowing the turning point to vary between subjects
    -   Thus making the predictions valid for a random future subject, not just the middle subject from the observed data
-   Incorporating a smooth transition directly into the model curve
    -   Allowing each subject to have their own level of curviness 😜

## Technicalities

> This presentation was created using the **Reveal.js** format in [Quarto](https://quarto.org/), using the [RStudio IDE](https://posit.co/products/open-source/rstudio/). Font and line colours according to UFS branding, and background image using image editor [GIMP](https://www.gimp.org/) by compositing images from CoPilot.

-   Interactive plots were created by feeding *ggplot2* plots to the *plotly* package, via the *ggplotly* function.

## Parameter comparison {.scrollable}

```{r}
AllSummaries |> mutate(Rhat = PSRFPointEst) |> select(1:5, Rhat) |>
  kable(digits = 3)
```

## References {.scrollable}