HIV-1 viral load after unstructured treatment interruption

Here we revisit the data set analysed in [2], which concerns the study of 72 perinatally HIV-infected children and adolescents. Some of the subjects present unstructured treatment interruption, which is a common phenomenon among HIV-infected people, due mainly to treatment fatigue. The number of observations from baseline (month 0) to month 24 ranges from 13 to 71. Out of 362 observations, 26 observations (7%) were below the detection limits and were censored at these values. [2] proposed the model: \[\begin{eqnarray*} y_{ij} = \beta_j + u_i + \varepsilon_{ij}, \end{eqnarray*}\]

where \(y_{ij}\) is the \(\log_{10}\) HIV RNA (HIV Ribonucleic acid, which is used to monitor the status of HIV) for subject \(i\) at time \(t_j\), using \(t_1=0\), \(t_2=1\), \(t_3=3\), \(t_4=6\), \(t_5=9\), \(t_6=12\), \(t_7=18\) and \(t_8=24\) months.

We allow for Student-\(t\) tails in the random effects and the errors (with zero locations) [1], and thus use Model 4 in Subsection 4.2 from [1] with the same prior structure as in (15) from [1]: \[\begin{eqnarray}\label{PriorSimulationII} \pi({\bf \beta},\sigma_{\varepsilon},\delta_{\varepsilon},\sigma,\delta) \propto \dfrac{\pi(\delta_{\varepsilon})\pi(\sigma)\pi(\delta)}{\sigma_{\varepsilon}}, \end{eqnarray}\]

where \(\pi(\delta_{\varepsilon})\) and \(\pi(\delta)\) are the weakly informative priors for the degrees of freedom of a Student-\(t\) distribution proposed in [3] and implemented in [4], \(\pi(\sigma)\) is a half-Cauchy density with mode \(0\) and unit scale parameter. propriety of the posterior is then guaranteed by Theorem 2 from [1].

[1] showed that there is clear evidence for fat tails in both distributions, especially the residual errors. The predictive criterion LPML also favours the general model and particularly penalizes model with normal residual errors. Clearly, ignoring the heavy tails affects the inference on the fixed effects, especially for \(\beta_1\).

The format of the tables is better in Google Chrome and Mozilla Firefox in our experience

References

  1. Flexible linear mixed models with improper priors for longitudinal and survival data

  2. Fast implementation for normal mixed effects models with censored response

  3. Bayesian modelling of skewness and kurtosis with Two-Piece Scale and shape distributions

  4. A weakly informative prior for the degrees of freedom of the t distribution

  5. TPSAS R Package

rm(list=ls())

# Required packages
library(tlmec)
## Loading required package: mvtnorm
library(mvtnorm)
library(compiler)
library(spBayes)
## Loading required package: coda
## Loading required package: magic
## Loading required package: abind
## Loading required package: Formula
library(TPSAS)
library(knitr)
enableJIT(3)
## [1] 3
# Extracting and preparing the data
data(UTIdata)
o <- order(UTIdata$Patid, UTIdata$Fup)
UTIdata <- UTIdata[o,]

# Responses and design matrix
y <- log10(UTIdata$RNA)
x <- cbind((UTIdata$Fup==0)+0, (UTIdata$Fup==1)+0, (UTIdata$Fup==3)+0, (UTIdata$Fup==6)+0, (UTIdata$Fup==9)+0, (UTIdata$Fup==12)+0, (UTIdata$Fup==18)+0, (UTIdata$Fup==24)+0)

table(UTIdata$Patid)
## 
##   C1  C10  C11  C12  C13  C14  C15  C16  C17  C18  C19   C2   C3   C4   C5 
##    6    6    4    8    4    4    5    2    3    5    5    6    3    3    7 
##   C6   C7   C8   C9 LA10 LA11 LA12 LA13  LA2  LA3  LA4  LA5  LA6  LA7  LA8 
##    1    4    5    6    5    6    7    4    4    6    4    5    4    5    4 
##  LA9 SD10 SD11 SD12 SD13 SD14 SD15 SD16 SD17  SD2  SD3  SD4  SD5  SD6  SD7 
##    4    4    6    2    4    6    6    3    5    7    8    8    7    3    8 
##  SD8  SD9   T1  T10  T11  T12  T13  T14  T15  T16  T17  T18  T19   T2  T20 
##    7    6    6    4    5    4    5    5    7    1    2    4    3    6    4 
##  T21  T22  T23  T24  T25  T26   T3   T4   T5   T6   T7   T9 
##    8    4    3    8    8    3    6    8    6    4    8    5
ind <- rep(1:72,table(UTIdata$Patid))
ind1 <- unique(ind)

# Censoring indicators
cens <- UTIdata$RNAcens
indcl = which(cens==1)
indcu = which(cens==2)
indo = which(cens==0)

#-----------------------------------------------------------------------------------
#Prior on delta Student-T
#-----------------------------------------------------------------------------------
# Parameters obtained by fitting the sample from \pi(\nu)
MLE <- c(0.84886337, 1.44326981, 0.05384718, 0.89544807)

# Prior on nu
logkdeT <- Vectorize(function(x) dtpsas(log(x),MLE[1],MLE[2],MLE[3],MLE[4],param="eps")/x)
curve(logkdeT,0,20,n=10000)

# log posterior distribution
lp <- function(par){
  if(par[9]>0 & par[10]>0 & par[11]>0 & par[12]>0){
    varo = dt((y[indo]-x[indo,]%*%par[1:8] + par[12+ind][indo])/par[9],df=par[10],log=T) - log(par[9])
    varcl = pt((y[indcl]-x[indcl,]%*%par[1:8] + par[12+ind][indcl])/par[9],df=par[10],log=T)
    varcu = log(1-pt((y[indcu]-x[indcu,]%*%par[1:8] + par[12+ind][indcu])/par[9],df=par[10]))
    varre = dt(par[12+ind1]/par[11],df=par[12],log=T) - log(par[11])
    LOGLIK = sum(varo) + sum(varcl) + sum(varcu) + sum(varre)
    ifelse(is.na(LOGLIK),
           return(-Inf),
           return( LOGLIK - log(par[9]) + dcauchy(par[11],log=T) + log(logkdeT(par[10])) + log(logkdeT(par[12])) ))
  }
  else return(-Inf)
}

lpc <- cmpfun(lp)

# Initial points
inits = c(rep(4,8),1,5,1,5,rep(0,72))
lpc(inits)
## [1] -616.6615
# Sampling from the posterior
n.batch <- 2200
batch.length <- 25

# Here only fixed effects and shape parameters are saved
# To save the posterior samples of the random effects, just increase the index in p.theta.samples[,1:12]
set.seed(1234)
chain <- adaptMetropGibbs(ltd=lpc, starting=inits, accept.rate=0.44, batch=n.batch, batch.length=batch.length, report=100, verbose = FALSE)$p.theta.samples[,1:12]

# Thinning and burning
burn = 5000
thin = 25
NS = n.batch*batch.length
index = seq(burn,NS,thin)

beta1p = chain[index,1]
beta2p = chain[index,2]
beta3p = chain[index,3]
beta4p = chain[index,4]
beta5p = chain[index,5]
beta6p = chain[index,6]
beta7p = chain[index,7]
beta8p = chain[index,8]
sigmap = chain[index,9]
deltap = chain[index,10]
sigma1p = chain[index,11]
delta1p = chain[index,12]

# Summaries of the posterior samples
hist(beta1p)

hist(beta2p)

hist(beta3p)

hist(beta4p)

hist(beta5p)

hist(beta6p)

hist(beta7p)

hist(beta8p)

hist(sigmap)

hist(deltap)

hist(sigma1p)

hist(delta1p)

SUM <- apply(chain[index,],2,summary)
colnames(SUM) <- c("beta1", "beta2", "beta3", "beta4", "beta5", "beta6", "beta7", "beta7", "sigma", "delta", "sigma_e", "delta_e")
kable(SUM,digits = 3)
beta1 beta2 beta3 beta4 beta5 beta6 beta7 beta7 sigma delta sigma_e delta_e
Min. 3.795 4.082 4.127 4.245 4.343 4.318 4.438 4.479 0.121 0.869 0.213 0.799
1st Qu. 4.113 4.328 4.365 4.530 4.631 4.603 4.779 4.802 0.186 1.183 0.489 1.885
Median 4.187 4.391 4.433 4.601 4.698 4.674 4.860 4.891 0.200 1.302 0.561 2.455
Mean 4.188 4.393 4.433 4.601 4.701 4.673 4.860 4.893 0.201 1.312 0.561 2.779
3rd Qu. 4.259 4.459 4.500 4.671 4.769 4.743 4.937 4.984 0.216 1.423 0.631 3.302
Max. 4.520 4.693 4.769 4.927 5.029 5.042 5.203 5.372 0.287 2.036 0.974 13.710