drizopoulos

set.seed(1234)
n <- 100 # number of subjects
K <- 8 # number of measurements per subject
t_max <- 15 # maximum follow-up time

# we constuct a data frame with the design: 
# everyone has a baseline measurment, and then measurements at random follow-up times
DF <- data.frame(id = rep(seq_len(n), each = K),
                 time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
                 sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))

drizopoulos

Ints <- list("log(serBilir)" = ~ drug, "log(serBilir)_slope" = ~ drug,
             "spiders_area" = ~ drug)

JMFit4 <- update(JMFit3, Interactions = Ints, priors = list(shrink_alphas = TRUE))
summary(JMFit4)

## Output
## [truncated]
##
## Survival Outcome:

drizopoulos

Forms <- list("log(serBilir)" = "value",
              "log(serBilir)" = list(fixed = ~ 1, random = ~ 1,
                                     indFixed = 2, indRandom = 2, name = "slope"),
              "spiders" = list(fixed = ~ 0 + year + I(year^2/2), random = ~ 0 + year,
                               indFixed = 1:2, indRandom = 1, name = "area"))

JMFit3 <- update(JMFit2, Formulas = Forms)
summary(JMFit3)

## Output

drizopoulos

MixedModelFit2 <- mvglmer(list(log(serBilir) ~ year + (year | id),
                               spiders ~ year + (1 | id)), data = pbc2,
                          families = list(gaussian, binomial))

JMFit2 <- mvJointModelBayes(MixedModelFit2, CoxFit1, timeVar = "year")
summary(JMFit2)

## Output
## Call:
## mvJointModelBayes(mvglmerObject = MixedModelFit2, coxphObject = CoxFit1, 

drizopoulos

JMFit1 <- mvJointModelBayes(MixedModelFit1, CoxFit1, timeVar = "year")
summary(JMFit1)

## Output
## Call:
## mvJointModelBayes(mvglmerObject = MixedModelFit1, coxphObject = CoxFit1, 
##     timeVar = "year")
## 
## Data Descriptives:
## Number of Groups: 312		Number of events: 169 (54.2%)

drizopoulos

pbc2.id$Time <- pbc2.id$years
pbc2.id$event <- as.numeric(pbc2.id$status != "alive")
CoxFit1 <- coxph(Surv(Time, event) ~ drug + age, data = pbc2.id, model = TRUE)

drizopoulos

MixedModelFit1 <- mvglmer(list(log(serBilir) ~ year + (year | id)), data = pbc2,
                          families = list(gaussian))

drizopoulos

plot(sfit, estimator = "mean", include.y = TRUE,
     conf.int = TRUE, fill.area = TRUE, col.area = "lightgrey")

drizopoulos

ND <- pbc2[pbc2$id == 2, ]
sfit <- survfitJM(jointFit, newdata = ND)
sfit

## Prediction of Conditional Probabilities for Event
## 	based on 200 Monte Carlo samples
##
## $`2`
##     times   Mean Median  Lower  Upper
## 0  8.8325 1.0000 1.0000 1.0000 1.0000

drizopoulos

runDynPred()