JM codes

Data generation

simulateJM <-   function (nsim, nsub, thetas, times, formulas, Data = NULL, censoring = NULL, max.FUtime = NULL) {

  if (is.null(max.FUtime))
    max.FUtime <-1e6

  # Function to compute the inverse survival function
  invS <- function (t, uu, i) {
    TD <- function (v) {
      # function to compute the time-dependent part for patient i at time v
      dd            <- Data[rep(i, length(v)), , drop = FALSE]
      dd[[timeVar]] <- pmax(v, 0)

      XX    <- model.matrix(formYx, data = dd)
      ZZ    <- model.matrix(formYz, data = dd)
      Y     <- as.vector(XX %*% betas +  rowSums(ZZ * b[rep(i, nrow(ZZ)), , drop = FALSE]))
      Y.aux <- 1/(1 + 1/exp(Y))
      out   <- (alpha ) * (nu + (1 - nu)*Y.aux)
      out
    }

    h <- function (s) {
      TD.i <- TD(s)
      exp(log(rho) + (rho - 1) * log(s) + eta.t[i] + TD.i)
    }
    integrate(h, lower = 0, upper = t)$value + log(uu)
  }

  # Coefficients
  betas  <- thetas$betas
  sigma  <- exp(thetas$sigma)
  D      <- exp(thetas$D)
  gamas  <- as.numeric(thetas$gamas)
  alpha  <- thetas$alpha
  Dalpha <- thetas$Dalpha
  rho    <- exp(thetas$rho)
  nu     <- exp(thetas$nu)/(1 + exp(thetas$nu))


  # Design matrices
  formYx  <- formulas$Yfixed
  formYz  <- formulas$Yrandom
  formT   <- formulas$Tfixed
  timeVar <- formulas$timeVar
  id      <- rep(1:nsub, each= length(times))
  times   <- rep(times, nsub)

  DD            <- Data[id, , drop = FALSE]
  DD[[timeVar]] <- times
  X   <- model.matrix(formYx, data = DD)
  Z   <- model.matrix(formYz, data = DD)
  W   <- model.matrix(formT, Data)
  ncz <- ncol(Z)


  # Simulate random effects
  b <- mvrnorm(nsub, rep(0, ncz), D)

  # Simulate event times
  eta.t <- if (!is.null(W)) as.vector(W %*% gamas) else rep(0, nsub)
  u     <- runif(nsub)
  trueTimes <- numeric(nsub)

  for (i in 1:nsub) {
    Root <- try(uniroot(invS, interval = c(1e-05, max.FUtime), uu = u[i], i = i)$root, TRUE)

    while(inherits(Root, "try-error")){
      b[i, ] <- c(mvrnorm(1, rep(0, ncz), D))
      Root   <- try(uniroot(invS, interval = c(1e-05, max.FUtime), uu =runif(1), i = i)$root, TRUE)
    }

    trueTimes[i] <- Root
  }

  # Simulate longitudinal responses
  eta.y    <- as.vector(X %*% betas + rowSums(Z * b[id, ]))
  mean.aux <- exp(eta.y)/(1 + exp(eta.y))
  y      <- rBEOI(nrow(DD), mean.aux, sigma, nu)
  Ctimes <- rep(censoring, length.out = nsub)
  Time   <- pmin(trueTimes, Ctimes)
  event  <- as.numeric(trueTimes <= Ctimes)

  ni            <- tapply(id, id, length)
  DD$EQ5D       <- y
  DD$tempo      <- rep(Time, ni)
  DD$delta      <- rep(event, ni)
  DD$X          <- id
  DD            <- DD[DD[[timeVar]] <= DD$tempo, ]
  row.names(DD) <- seq_len(nrow(DD))

  DD
}

Example

dadosSim = simulateJM(nsim       = 1,
                      nsub       = 500, 
                      thetas     = initi.theta.BEOI.sim, 
                      times      = c(0,15,90,180,360,540), 
                      formulas   = form, 
                      Data       = dadosLong,  
                      max.FUtime = NULL, 
                      censoring  = 570)

dadosSimB <- dadosSim

dadosSimB %<>%
    mutate(EQ5D = ifelse(EQ5D == 1, (nrow(dados) - 1 + 0.5)/nrow(dados), EQ5D), # substitui 1's por ((n-1) + 0.5)/n
           logito.EQ5D = log(EQ5D/(1-EQ5D))) %>%
    arrange(idpaciente, visita)


dadosSim.id <- dadosSim[!duplicated(dadosSim$X),]

head(dadosSim, n = 5)

Exploratory analysis:

Distribution of UI

Longitunial profiles

## `geom_smooth()` using formula = 'y ~ x'

Survival curve:

Model fitting

######################
# No scaling        #
######################
fit.betaJM.id  <- try(fit.JM(lmeObj.Beta, 
                              survObj, 
                              model = "beta", 
                              QH = 10, 
                              QL = 10, 
                              lag = 0,
                              timeVar = "visita", 
                              init.theta = initi.theta.beta, 
                              imp = FALSE), TRUE)

fit.BEOIJM.id  <- try(fit.JM(lmeObj, 
                             survObj,
                             model = "betaInf", 
                             QH = 10, 
                             QL = 10,
                             lag = 0, 
                             timeVar = "visita", 
                             init.theta = initi.theta.BEOI,
                             imp = FALSE), TRUE)

fit.lnJM.id    <- try(fit.JM(lmeObj.ln, 
                             survObj,
                             model = "normal",
                             QH = 10,
                             QL = 10, 
                             lag = 0, 
                             timeVar = "visita", 
                             init.theta = initi.theta.ln,
                             imp = FALSE), TRUE)

####################################
# Scaling factors: Typical values  #
#################################### 

fit.betaJM.vt  <- try(fit.JM(lmeObj.Beta,
                             survObj, 
                             model = "beta",
                             QH = 10, 
                             QL = 10, 
                             lag = 0, 
                             timeVar = "visita", 
                             init.theta = initi.theta.beta, 
                             imp = FALSE, 
                             precond = 'reescala'),TRUE)

fit.BEOIJM.vt  <- try(fit.JM(lmeObj, 
                             survObj,
                             model = "betaInf", 
                             QH = 10, 
                             QL = 10, 
                             lag = 0, 
                             timeVar = "visita", 
                             init.theta = initi.theta.BEOI, 
                             imp = FALSE,
                             precond = 'reescala'),TRUE)
  
fit.lnJM.vt    <- try(fit.JM(lmeObj.ln, 
                             survObj, 
                             model = "normal",
                             QH = 10, 
                             QL = 10, 
                             lag = 0, 
                             timeVar = "visita", 
                             init.theta = initi.theta.ln, 
                             imp = FALSE, 
                             precond = 'reescala'),TRUE)

############################
# Scaling factors: Jacobi  #
############################
  
fit.betaJM.jc  <- try(fit.JM(lmeObj.Beta, 
                             survObj, 
                             model = "beta",
                             QH = 10,
                             QL = 10, 
                             lag = 0, 
                             timeVar = "visita", 
                             init.theta = initi.theta.beta,
                             imp = FALSE,
                             precond = 'jacobi'),TRUE)

fit.BEOIJM.jc  <- try(fit.JM(lmeOb,
                             survObj,
                             model = "betaInf", 
                             QH = 10,
                             QL = 10, 
                             lag = 0, 
                             timeVar = "visita", 
                             init.theta = initi.theta.BEOI,
                             imp = FALSE,
                             precond = 'jacobi'),TRUE)

fit.lnJM.jc    <- try(fit.JM(lmeObj.ln, 
                             survObj, 
                             model = "normal",
                             QH = 10, 
                             QL = 10, 
                             lag = 0, 
                             timeVar = "visita", 
                             init.theta = initi.theta.ln,  
                             imp = FALSE, 
                             precond = 'jacobi'),TRUE)

Checking the output

fit.lnJM.jc 

## $fit
## $fit$par
##           betas.(Intercept)                betas.visita 
##                 0.586542727                 0.002556241 
##    betas.quimioSim-Curativa           gamas.(Intercept) 
##                 0.192065807                -2.360993353 
## gamas.tipo.adm.utiPlanejada                       alpha 
##                -0.914694044                -0.187746815 
##                         rho                           D 
##                -0.635201237                 1.771809639 
##                       sigma 
##                 1.427942553 
## attr(,"skeleton")
## $betas
##        (Intercept)             visita quimioSim-Curativa 
##        0.552344799        0.002523929        0.268761350 
## 
## $gamas
##           (Intercept) tipo.adm.utiPlanejada 
##            -2.5680063            -0.9497281 
## 
## $alpha
## [1] 1e-04
## 
## $rho
## [1] -0.5592707
## 
## $D
##                 [,1]
## (Intercept) 1.716548
## 
## $sigma
## [1] 1.431236
## 
## attr(,"class")
## [1] "relistable" "list"      
## 
## $fit$value
## [1] -7409.61
## 
## $fit$counts
## function gradient 
##       20       16 
## 
## $fit$convergence
## [1] 0
## 
## $fit$message
## NULL
## 
## $fit$hessian
##                             betas.(Intercept)  betas.visita
## betas.(Intercept)               -4.471656e+01    -4513.7898
## betas.visita                    -4.513790e+03 -3079043.9758
## betas.quimioSim-Curativa        -1.890355e+01    -1570.3777
## gamas.(Intercept)                2.679013e+01     2605.0495
## gamas.tipo.adm.utiPlanejada      1.089209e+01     1365.6734
## alpha                            5.584468e+01    25067.1316
## rho                              7.358830e+01    12417.5187
## D                               -7.826679e+00    -3566.3640
## sigma                            5.371893e-02     -173.3976
##                             betas.quimioSim-Curativa gamas.(Intercept)
## betas.(Intercept)                         -18.903548          26.79013
## betas.visita                            -1570.377712        2605.04953
## betas.quimioSim-Curativa                  -18.903548          13.24299
## gamas.(Intercept)                          13.242987        -325.99092
## gamas.tipo.adm.utiPlanejada                 2.263828        -147.24787
## alpha                                      31.457298        -225.79433
## rho                                        35.035417        -952.21364
## D                                          -3.357063         -11.19818
## sigma                                      -1.593202          11.54508
##                             gamas.tipo.adm.utiPlanejada       alpha         rho
## betas.(Intercept)                             10.892089    55.84468    73.58830
## betas.visita                                1365.673364 25067.13163 12417.51872
## betas.quimioSim-Curativa                       2.263828    31.45730    35.03542
## gamas.(Intercept)                           -147.247873  -225.79433  -952.21364
## gamas.tipo.adm.utiPlanejada                 -147.247874   -81.39342  -463.16754
## alpha                                        -81.393418 -1479.69051  -992.01336
## rho                                         -463.167544  -992.01336 -3111.50503
## D                                             -5.594832    90.73939   -13.96673
## sigma                                          4.186833  -148.09112    22.75321
##                                        D         sigma
## betas.(Intercept)              -7.826679  5.371893e-02
## betas.visita                -3566.364024 -1.733976e+02
## betas.quimioSim-Curativa       -3.357063 -1.593202e+00
## gamas.(Intercept)             -11.198182  1.154508e+01
## gamas.tipo.adm.utiPlanejada    -5.594832  4.186833e+00
## alpha                          90.739390 -1.480911e+02
## rho                           -13.966732  2.275321e+01
## D                             -62.630718 -1.247113e+02
## sigma                        -124.711309 -2.778801e+03
## 
## 
## $out
##                              Estimativas Erros padr?o      p-valor
## betas.(Intercept)            0.586542727 0.2126219618 5.804638e-03
## betas.visita                 0.002556241 0.0006633024 1.162945e-04
## betas.quimioSim-Curativa     0.192065807 0.3064218081 5.307890e-01
## gamas.(Intercept)           -2.360993353 0.1883708704 0.000000e+00
## gamas.tipo.adm.utiPlanejada -0.914694044 0.1165358677 4.218847e-15
## alpha                       -0.187746815 0.0366426830 2.995671e-07
## rho                         -0.635201237 0.0670636814 0.000000e+00
## D                            5.881487103 0.1467153368 0.000000e+00
## sigma                        4.170110580 0.0203847218 0.000000e+00
## 
## $cond
## [1] 112.8826