brea_gmrf.R
In brea: Bayesian Recurrent Events Analysis

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

## -----------------------------------------------------------------------------
time <- c(6,6,6,6,7,9,10,10,11,13,16,17,19,20,22,23,25,32,32,34,35,
          1,1,2,2,3,4,4,5,5,8,8,8,8,11,11,12,12,15,17,22,23)

event <- c(0,1,1,1,1,0,0,1,0,1,1,0,0,0,1,1,0,0,0,0,0,
           1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)

treat <- c(rep(1,21),rep(2,21))

## -----------------------------------------------------------------------------
N <- sum(time) # total number of person-time observations
x <- matrix(0,N,3)  # columns id, time, and treatment group
colnames(x) <- c("id","t","treat")
y <- matrix(0,N,1) # only one column since only one competing risk

## -----------------------------------------------------------------------------
next_row <- 1  # next available row in the person-time matrices
for (i in seq_along(time)) {
  rows <- next_row:(next_row + time[i] - 1)  # person-time rows to fill
  x[rows,1] <- rep(i,time[i]) # subject id number is constant for each person
  x[rows,2] <- seq_len(time[i])  # discrete time is integers 1,2,...,time[i]
  x[rows,3] <- rep(treat[i],time[i])  # treatment group is constant
  y[rows,1] <- c(rep(0,time[i] - 1),event[i])  # outcome is 0's until study time
  next_row <- next_row + time[i]  # increment the next available row pointer
}

## -----------------------------------------------------------------------------
x_brea <- matrix(0,N,2)
c_k <- seq(0,36,3) # step boundaries
x_brea[,1] <- cut(x[,2],c_k,labels=FALSE) # grouped time t
x_brea[,2] <- x[,3] # treatment group

## -----------------------------------------------------------------------------
priors <- list(list("gmrf",3,0.01),list("cat",4))

## ----eval=FALSE---------------------------------------------------------------
#  brea_mcmc(x, y, priors = NULL, S = 1000, B = 100, n = NULL, K = NULL,
#            store_re = FALSE,  block_mh = TRUE)

## -----------------------------------------------------------------------------
library(brea)
set.seed(1234)
fit <- brea_mcmc(x_brea,y,priors,10000,1000)

## -----------------------------------------------------------------------------
str(fit)

## -----------------------------------------------------------------------------
b_t_post_medians <- apply(fit$b_m_s[[1]][1,,],1,median)
b_t_post_medians

## -----------------------------------------------------------------------------
hr <- exp(b_t_post_medians)
hr

## ----fig.width=6.5,fig.height=4.9---------------------------------------------
par(cex=0.66,mgp=c(1.75,0.5,0),mai=c(0.4,0.4,0.5,0.1))

# setup the plotting window:
plot(NA,NA,xlim=range(c_k),ylim=range(hr),log="y",
     xlab="Time at Risk (weeks)",ylab="Relative Discrete Hazard",
     main="Time Effect on Hazard of Remission Cessation")

# add each of the K=12 steps:
for (k in seq_len(length(c_k)-1)) {
  lines(c(c_k[k],c_k[k+1]),rep(hr[k],2))
}


## -----------------------------------------------------------------------------
c_k6 <- seq(0,36,6) # 6-week apart step boundaries
x_brea[,1] <- cut(x[,2],c_k6,labels=FALSE) # grouped time t
set.seed(1234)
fit6 <- brea_mcmc(x_brea,y,priors,10000,1000)

c_k2 <- seq(0,36,2) # 2-week apart step boundaries
x_brea[,1] <- cut(x[,2],c_k2,labels=FALSE) # grouped time t
set.seed(1234)
fit2 <- brea_mcmc(x_brea,y,priors,10000,1000)

c_k1 <- seq(0,36,1) # 1-week apart step boundaries
x_brea[,1] <- cut(x[,2],c_k1,labels=FALSE) # grouped time t
set.seed(1234)
fit1 <- brea_mcmc(x_brea,y,priors,10000,1000)


## -----------------------------------------------------------------------------
hr6 <- exp(apply(fit6$b_m_s[[1]][1,,],1,median))
hr2 <- exp(apply(fit2$b_m_s[[1]][1,,],1,median))
hr1 <- exp(apply(fit1$b_m_s[[1]][1,,],1,median))

## ----fig.width=6.5,fig.height=4.9---------------------------------------------
par(cex=0.66,mgp=c(1.75,0.5,0),mai=c(0.4,0.4,0.5,0.1))

# setup the plotting window:
plot(NA,NA,xlim=range(c_k),ylim=range(c(hr,hr6,hr2,hr1)),log="y",
     xlab="Time at Risk (weeks)",ylab="Relative Discrete Hazard",
     main="Time Effect on Hazard of Remission Cessation")

for (k in seq_len(length(c_k6)-1)) {
  lines(c(c_k6[k],c_k6[k+1]),rep(hr6[k],2),col="red")
}

for (k in seq_len(length(c_k)-1)) {
  lines(c(c_k[k],c_k[k+1]),rep(hr[k],2))
}

for (k in seq_len(length(c_k2)-1)) {
  lines(c(c_k2[k],c_k2[k+1]),rep(hr2[k],2),col="green")
}

for (k in seq_len(length(c_k1)-1)) {
  lines(c(c_k1[k],c_k1[k+1]),rep(hr1[k],2),col="blue")
}

legend("topleft",
       legend=c("6-week Steps","3-week Steps","2-week Steps","1-week Steps"),
       col=c("red","black","green","blue"),lty=1)

## -----------------------------------------------------------------------------
c_k <- seq(0,36,3) # step boundaries
x_brea[,1] <- cut(x[,2],c_k,labels=FALSE) # grouped time t
set.seed(1234)
fit_rwm <- brea_mcmc(x_brea,y,priors,10000,1000,block_mh = FALSE)

## -----------------------------------------------------------------------------
# make the effectiveSize function from the coda package available:
library(coda)

# MCMC sample size:
S <- 10000

# MCMC efficiency for the default Metropolis-Hastings algorithm:
eff_mh <- apply(fit$b_m_s[[1]][1,,],1,effectiveSize)/S

# MCMC efficiency for the random walk Metropolis algorithm:
eff_rwm <- apply(fit_rwm$b_m_s[[1]][1,,],1,effectiveSize)/S

## -----------------------------------------------------------------------------
# efficiencies of the two algorithms:
eff_rwm
eff_mh
summary(eff_rwm)
summary(eff_mh)

# improvement factor for block Metropolis-Hastings over random walk Metropolis:
eff_mh/eff_rwm
summary(eff_mh/eff_rwm)

## ----fig.width=6.5,fig.height=4.9---------------------------------------------
par(cex=0.66,mgp=c(1.75,0.5,0),mai=c(0.4,0.4,0.5,0.1))

K <- 12
plot(NA,NA,xlim=c(1,K),ylim=c(0,0.7),xlab="Parameter Number",
     ylab="MCMC Efficiency",xaxt="n")
axis(1,1:K,1:K)
abline(h=0)
lines(1:K,eff_rwm,col="red",type="b",pch=20)
lines(1:K,eff_mh,col="blue",type="b",pch=20)
legend("topleft",
       legend=c("Block Metropolis-Hastings","Random Walk Metropolis"),
       col=c("blue","red"),lty=1,pch=20)