# joineRmeta

#### Maria Sudell

#### 2018-03-02

# Introduction

The joineRmeta package implements methods for analysing joint longitudinal and time-to-event data from multiple studies. The responses for each subject should consist of a single continuous longitudinal outcome and a single possibly censored time-to-event outcome. This package provides methods to perform the second stage of a two stage meta-analysis of joint model fits, where model parameters from study specific fits are pooled using standard meta-analytic techniques, and methods to perform a one stage meta-analysis of the data, where data from all studies is analysed simultaneously. Methods are also included to easily plot joint data from multiple studies, and to simulate joint data from multiple studies.

To load the library, use code

    library("joineRmeta")

    ## Loading required package: lme4

    ## Loading required package: Matrix

    ## Loading required package: survival

    ## Loading required package: JM

    ## Loading required package: MASS

    ## Loading required package: nlme

    ##
    ## Attaching package: 'nlme'

    ## The following object is masked from 'package:lme4':
    ##
    ##     lmList

    ## Loading required package: splines

# Data

This package contains several functions beneficial when first dealing with multi-study joint data in a project. These include a function to change multi-study data supplied in a range of formats to the `jointdata` format introduced in the `joineR` package, and a function to remove longitudinal information recorded after the event of interest if the event of interest is not terminal.

## Data Simulation

The `joineRmeta` package contains the `simjointmeta` function that allows simulation of three level joint longitudinal and survival data (longitudinal measurements (level 1) nested within individuals (level 2) nested within studies (level 3)). To see the help file for the `simjointmeta` function use the following command:

    help("simjointmeta", package = "joineRmeta")

This function is an extension of one written by the `joineR` package team to the multi-study case. This function simulates data for multiple studies, for one continuous longitudinal and one survival outcome, under a joint model with a random effects only association structure. The association structure can either be proportional (one common association parameter for all shared random effects at a given level) or separate (separate association parameters for each shared random effects) by setting `sepassoc` to be `FALSE` or `TRUE` respectively. Random effects can be specified at just the individual level, or at both the individual and the study levels. At the individual level, a random intercept can be specified, or both a random intercept and a random slope. If study level random effects are included, a study level random intercept, a study level random treatment effect, or both can be specified. The covariance matrices for the individual level and the study level random effects (if specified) are supplied to arguments `sigb_ind` and `sigb_stud`. The variation of the longitudinal measurement error term is set using argument `vare`.

The longitudinal outcome is generated under a linear mixed effects model, with a fixed intercept, time (slope) and binary indicator variable (which is labelled `treat` to represent treatment assignment to one of two treatment groups). The coefficients for these fixed effects are supplied in the function call to argument `beta1`, in order intercept, slope, treatment effects. The only fixed effect influencing the survival times is again the binary treatment indicator variable, the coefficient for which is supplied in the function call to argument `beta2`.

If separate association parameters have been set for each included random effects, the association parameters for each of the individual level random effects should be supplied as a vector to `gamma_ind`, and the association parameters for each of the study level random effects should be supplied as a vector to `gamma_stud`. However if a proportional association has been selected, then a single association parameter for the individual level random effects should be supplied, and a single association parameter for the study level random effects.

The survival data is generated under an exponential distribution if the individual level random effects do not contain a random slope (time) term, and a Gompertz distribution otherwise (see Austin (2012), and Bender, Augustin, and Blettner (2005)). These distributions can be controlled using the `theta0` and `theta1` parameters. If the individual level random effects contain only a random intercept, then the mean of the exponentially distributed survival times is equal to one over the log of `theta0`. However if the individual level random effects contain both a random intercept and a random slope (time) term, then the survival times are generated under a Gompertz distribution, and so the values of `theta0` and `theta1`, can be determined using the extreme value distribution (a description of how to do this is given in Bender, Augustin, and Blettner (2005)). The survival times can also be censored by setting `censoring` equal to `TRUE`. If censored, the censoring times are expoentially distributed, with lambda parameter equal to `censlam`. The survival times could also be set to be truncated by setting `truncation` to `TRUE`, the default is to truncate the survival times at the largest value of `longmeasuretimes` i.e. the last longitudinal time-point. Note that longitudinal observations simulated at times after the simulated event time for each individual are discarded by the function and not returned in the final datasets (i.e. the event being simulated is assumed terminal)

The number of studies to simulate data for is supplied to argument `k`. The number of individuals to simulate in each study is supplied to argument `n`, which should be a numeric of length equal to `k`. The number of longitudinal measurement times is supplied to argument `ntms`. The time points for the longitudinal measurements can either be supplied to argument `longmeasuretimes`, or if they are not, the function sets them to integer values from 0 to `ntms`.

Various ways exist to introduce between study heterogeneity into the data. Firstly study level random effects can be specified in the function call. Secondly it is possible to allow different distributions of survival times between simulated studies by suppling a vector of `theta0` and `theta1` values of length equal to the number of studies `k` rather than single values common across studies for each argument. Thirdly it is possible to allow different censoring time distributions between simulated studies by suppling a vector of `censlam` values of length equal to the number of studies `k` rather than a single value common across studies for each argument. Finally whether separate or proportional association has been selected, different association parameters can be supplied for each study in a list of length equal to the number of studies in the simulated data. If separate association, this is a list of vectors for each of `gamma_ind` and `gamma_stud`, where each vector is of length equal to the number of random effects at each level. If proportional association this is a list of single values for each of `gamma_ind` and `gamma_stud`.

We will give examples of simulation of several different datasets to demonstrate the capabilities of this data simulation function. The first example contains no study level random effects, with proportional assocation, with common parameters across studies.

    exampledat1 <- simjointmeta(k = 5, n = rep(500, 5), sepassoc = FALSE,
                                ntms = 5,longmeasuretimes = c(0, 1, 2, 3, 4),
                                beta1 = c(1, 2, 3), beta2 = 1,
                                rand_ind = "intslope", rand_stud = NULL,
                                gamma_ind = 1,
                                sigb_ind = matrix(c(1, 0.5, 0.5, 1.5),nrow = 2),
                                vare = 0.01, theta0 = -3, theta1 = 1,
                                censoring = TRUE, censlam = exp(-3),
                                truncation = FALSE,
                                trunctime = max(longmeasuretimes))

    ## 73.8 % experienced event
    ## 72.8 % experienced event
    ## 72.4 % experienced event
    ## 73.8 % experienced event
    ## 74.6 % experienced event

When run, the function prints the percentage of events observed in each study. This data is also returned by the function and can be accessed using the following code, which returns a list of the percentage events in each study.

    exampledat1$percentevent

The longitudinal and survival data are returned seperately by this function. They are returned as lists of datasets, one for each study specified to be simulated. These can be extracted, and using the `tojointdata` function, later can easily be transformed into a `jointdata` object, which is necessary to be able to use many of the functions available in this package. The function `tojointdata` is further discussed below.

    str(exampledat1$longitudinal)

    ## List of 5
    ##  $ :'data.frame':    1325 obs. of  7 variables:
    ##   ..$ id       : num [1:1325] 1 2 2 2 2 3 3 4 4 5 ...
    ##   ..$ Y        : num [1:1325] 2.62 1.29 4.21 7.25 10.26 ...
    ##   ..$ time     : num [1:1325] 0 0 1 2 3 0 1 0 1 0 ...
    ##   ..$ study    : int [1:1325] 1 1 1 1 1 1 1 1 1 1 ...
    ##   ..$ intercept: num [1:1325] 1 1 1 1 1 1 1 1 1 1 ...
    ##   ..$ treat    : num [1:1325] 1 0 0 0 0 1 1 0 0 1 ...
    ##   ..$ ltime    : num [1:1325] 0 0 1 2 3 0 1 0 1 0 ...
    ##  $ :'data.frame':    1404 obs. of  7 variables:
    ##   ..$ id       : num [1:1404] 501 501 501 501 502 503 503 503 503 504 ...
    ##   ..$ Y        : num [1:1404] 3.49 5.74 7.87 10.32 3.61 ...
    ##   ..$ time     : num [1:1404] 0 1 2 3 0 0 1 2 3 0 ...
    ##   ..$ study    : int [1:1404] 2 2 2 2 2 2 2 2 2 2 ...
    ##   ..$ intercept: num [1:1404] 1 1 1 1 1 1 1 1 1 1 ...
    ##   ..$ treat    : num [1:1404] 1 1 1 1 1 1 1 1 1 1 ...
    ##   ..$ ltime    : num [1:1404] 0 1 2 3 0 0 1 2 3 0 ...
    ##  $ :'data.frame':    1373 obs. of  7 variables:
    ##   ..$ id       : num [1:1373] 1001 1001 1002 1003 1003 ...
    ##   ..$ Y        : num [1:1373] 3.44 7.99 4 4.12 8.15 ...
    ##   ..$ time     : num [1:1373] 0 1 0 0 1 0 0 0 1 2 ...
    ##   ..$ study    : int [1:1373] 3 3 3 3 3 3 3 3 3 3 ...
    ##   ..$ intercept: num [1:1373] 1 1 1 1 1 1 1 1 1 1 ...
    ##   ..$ treat    : num [1:1373] 1 1 1 1 1 1 1 0 0 0 ...
    ##   ..$ ltime    : num [1:1373] 0 1 0 0 1 0 0 0 1 2 ...
    ##  $ :'data.frame':    1398 obs. of  7 variables:
    ##   ..$ id       : num [1:1398] 1501 1501 1501 1501 1501 ...
    ##   ..$ Y        : num [1:1398] 2.32 3.83 5.43 6.75 8.35 ...
    ##   ..$ time     : num [1:1398] 0 1 2 3 4 0 1 2 0 0 ...
    ##   ..$ study    : int [1:1398] 4 4 4 4 4 4 4 4 4 4 ...
    ##   ..$ intercept: num [1:1398] 1 1 1 1 1 1 1 1 1 1 ...
    ##   ..$ treat    : num [1:1398] 1 1 1 1 1 0 0 0 0 1 ...
    ##   ..$ ltime    : num [1:1398] 0 1 2 3 4 0 1 2 0 0 ...
    ##  $ :'data.frame':    1345 obs. of  7 variables:
    ##   ..$ id       : num [1:1345] 2001 2002 2002 2002 2003 ...
    ##   ..$ Y        : num [1:1345] 1.734 2.768 5.945 8.878 0.953 ...
    ##   ..$ time     : num [1:1345] 0 0 1 2 0 1 2 3 0 0 ...
    ##   ..$ study    : int [1:1345] 5 5 5 5 5 5 5 5 5 5 ...
    ##   ..$ intercept: num [1:1345] 1 1 1 1 1 1 1 1 1 1 ...
    ##   ..$ treat    : num [1:1345] 0 1 1 1 1 1 1 1 1 1 ...
    ##   ..$ ltime    : num [1:1345] 0 0 1 2 0 1 2 3 0 0 ...

    str(exampledat1$survival)

    ## List of 5
    ##  $ :'data.frame':    500 obs. of  5 variables:
    ##   ..$ id      : num [1:500] 1 2 3 4 5 6 7 8 9 10 ...
    ##   ..$ survtime: num [1:500] 1.81e-05 3.25 1.73 1.67 8.94e-01 ...
    ##   ..$ cens    : num [1:500] 1 1 0 1 1 0 1 1 1 0 ...
    ##   ..$ study   : int [1:500] 1 1 1 1 1 1 1 1 1 1 ...
    ##   ..$ treat   : int [1:500] 1 0 1 0 1 0 1 1 0 1 ...
    ##  $ :'data.frame':    500 obs. of  5 variables:
    ##   ..$ id      : num [1:500] 501 502 503 504 505 506 507 508 509 510 ...
    ##   ..$ survtime: num [1:500] 3.9 0.663 3.612 0.282 5.561 ...
    ##   ..$ cens    : num [1:500] 1 1 1 1 1 1 1 1 1 0 ...
    ##   ..$ study   : int [1:500] 2 2 2 2 2 2 2 2 2 2 ...
    ##   ..$ treat   : int [1:500] 1 1 1 1 0 1 0 1 1 0 ...
    ##  $ :'data.frame':    500 obs. of  5 variables:
    ##   ..$ id      : num [1:500] 1001 1002 1003 1004 1005 ...
    ##   ..$ survtime: num [1:500] 1.146 0.1972 1.2637 0.4497 0.0959 ...
    ##   ..$ cens    : num [1:500] 1 1 1 1 1 1 1 0 1 1 ...
    ##   ..$ study   : int [1:500] 3 3 3 3 3 3 3 3 3 3 ...
    ##   ..$ treat   : int [1:500] 1 1 1 1 1 0 0 0 0 1 ...
    ##  $ :'data.frame':    500 obs. of  5 variables:
    ##   ..$ id      : num [1:500] 1501 1502 1503 1504 1505 ...
    ##   ..$ survtime: num [1:500] 6.972 2.44 0.335 1.335 23.39 ...
    ##   ..$ cens    : num [1:500] 0 1 0 1 0 1 1 0 1 1 ...
    ##   ..$ study   : int [1:500] 4 4 4 4 4 4 4 4 4 4 ...
    ##   ..$ treat   : int [1:500] 1 0 0 1 1 1 1 0 0 1 ...
    ##  $ :'data.frame':    500 obs. of  5 variables:
    ##   ..$ id      : num [1:500] 2001 2002 2003 2004 2005 ...
    ##   ..$ survtime: num [1:500] 0.417 2.252 3.714 0.692 0.322 ...
    ##   ..$ cens    : num [1:500] 1 1 1 1 1 1 1 1 1 0 ...
    ##   ..$ study   : int [1:500] 5 5 5 5 5 5 5 5 5 5 ...
    ##   ..$ treat   : int [1:500] 0 1 1 1 1 1 0 0 1 1 ...

The second example shows both the inclusion of study level random effects, and also the ability to assign separate association parameters for each random effect at each level in the data.

    exampledat2 <- simjointmeta(k = 5, n = rep(500, 5), sepassoc = TRUE, ntms = 5,
                              longmeasuretimes = c(0, 1, 2, 3, 4), beta1 = c(1, 2, 3),
                              beta2 = 1, rand_ind = "intslope", rand_stud = "inttreat",
                              gamma_ind = c(0.5, 1), gamma_stud = c(0.5, 1),
                              sigb_ind = matrix(c(1, 0.5, 0.5, 1.5), nrow = 2),
                              sigb_stud = matrix(c(1, 0.5, 0.5, 1.5), nrow = 2),
                              vare = 0.01, theta0 = -3, theta1 = 1, censoring = TRUE,
                              censlam = exp(-3), truncation = FALSE,
                              trunctime = max(longmeasuretimes))

    ## 73.2 % experienced event
    ## 85.6 % experienced event
    ## 73.8 % experienced event
    ## 76.6 % experienced event
    ## 69.4 % experienced event

The final example extends the simulated data assigned to `exampledat2` to the case where separate parameters for the association parameters, the distribution of the survival times and the distribution of the censoring times are supplied for each study.

    gamma_ind_set <- list(c(0.5, 1), c(0.4, 0.9), c(0.6, 1.1), c(0.5, 0.9), c(0.4, 1.1))
    gamma_stud_set <- list(c(0.6, 1.1), c(0.5, 1), c(0.5, 0.9), c(0.4, 1.1), c(0.4, 0.9))
    censlamset <- c(exp(-3), exp(-2.9), exp(-3.1), exp(-3), exp(-3.05))
    theta0set <- c(-3, -2.9, -3, -2.9, -3.1)
    theta1set <- c(1, 0.9, 1.1, 1, 0.9)

    exampledat2 <- simjointmeta(k = 5, n = rep(500, 5), sepassoc = TRUE, ntms = 5,
                                longmeasuretimes = c(0, 1, 2, 3, 4), beta1 = c(1, 2, 3),
                                beta2 = 1, rand_ind = "intslope", rand_stud = "inttreat",
                                gamma_ind = gamma_ind_set, gamma_stud = gamma_stud_set,
                                sigb_ind = matrix(c(1, 0.5, 0.5, 1.5), nrow = 2),
                                sigb_stud = matrix(c(1, 0.5, 0.5, 1.5), nrow = 2),
                                vare = 0.01, theta0 = theta0set, theta1 = theta1set,
                                censoring = TRUE, censlam = censlamset, truncation = FALSE,
                                trunctime = max(longmeasuretimes))

    ## 77.2 % experienced event
    ## 69.4 % experienced event
    ## 82 % experienced event
    ## 75.8 % experienced event
    ## 65.6 % experienced event

## Example Data

Throughout this vignette we will refer to simulated example datasets `simdat2`, and `simdat3`, which are available in the package. Details of these datasets can be found by running the commands:

    help("simdat2", package = "joineRmeta")
    help("simdat3", package = "joineRmeta")

Briefly, `simdat2` contains simulated joint data for 3 studies each containing 100 individuals. This is a subset of simulated dataset `simdat`, available in the `joineRmeta` package. The data in `simdat3` is available solely to demonstrate the `removeafter` function below, whilst the data in `simdat` and `simdat2` exist in the package as simulated example datasets for cases with a small or larger number of identified studies.

We will proceed with the `simdat2` dataset unless otherwise stated. The `simdat2` dataset is supplied as a list, the elements of which are a list of long format longitudinal datasets, a list of time-to-event datasets, and a list of the percentage of individuals experiencing an event in each study. Both of `simdat`, and `simdat3` were generated using the function `simjointmeta` available in the `joineRmeta` package, and have been supplied in the format output by this function. Note that `simdat3` was produced by manually stepping through the `simjointmeta` function, to avoid longitudinal observations occuring after the event time from being discarded. Also note that `simdat2` is a subset of `simdat` rather than a completely independently simulated dataset. The structure of the data can be viewed using:

    data("simdat2")
    str(simdat2)

## Reformatting to jointdata format

We will reformat this data into the `jointdata` format required for the functions in this dataset. The `jointdata` function is available as part of the `joineR` package. A function `tojointdata` is available in `joineRmeta` that aids multistudy data provided in a variety of forms to be reformatted into a `jointdata` object. In this case we have `simdat2` whose first element is a list of long format longitudinal datasets, and whose second element is a list of time-to-event datasets. We change this to a `jointdata` object using the following code:

    jointdat<-tojointdata(longitudinal = simdat2$longitudinal,
      survival = simdat2$survival, id = "id", longoutcome = "Y",
      timevarying = c("time","ltime"), survtime = "survtime", cens = "cens",
      time = "time")

Because we had separate longitudinal and survival datasets we used the parameters `longitudinal` and `survival` in the `tojointdata` function call. If we had had one dataset in wide or long format containing both the longitudinal and survival information we would have supplied this to an argument `dataset`. If we had had a dataset or datasets of baseline data in addition to the longitudinal and survival datasets, these would have been supplied to argument `baseline`. Data can be supplied as lists of identically formatted (long or wide) datasts where the list is of length equal to the number of studies in the dataset, or as metadatasets containing all data, or all longitudinal or survival or baseline data from each study. The function assumes that variables that are identically named between datasets concern the same variables, therefore datasets should be checked before reformatting to ensure consistency in naming of variables, and of format (long or wide) if separate datasets for each study are supplied. For more information concerning the `tojointdata` function run

    help("tojointdata", package = "joineRmeta")

If we now examine the structure of this joint data, we can see that two variables (`study` and `treat`, respectively the study membership and treatment assignment variables) that should be factors are currently numeric.

    str(jointdat)

    ## List of 6
    ##  $ subject     : num [1:300] 1 2 3 4 6 12 13 17 20 22 ...
    ##  $ longitudinal:'data.frame':    810 obs. of  4 variables:
    ##   ..$ id   : num [1:810] 1 2 3 4 6 6 6 6 6 12 ...
    ##   ..$ Y    : num [1:810] 0.847 3.04 -0.252 1.976 1.434 ...
    ##   ..$ time : num [1:810] 0 0 0 0 0 1 2 3 4 0 ...
    ##   ..$ ltime: num [1:810] 0 0 0 0 0 1 2 3 4 0 ...
    ##  $ survival    :'data.frame':    300 obs. of  3 variables:
    ##   ..$ id      : num [1:300] 1 2 3 4 6 12 13 17 20 22 ...
    ##   ..$ survtime: num [1:300] 0.346 0.595 0.489 0.323 6 ...
    ##   ..$ cens    : num [1:300] 1 0 0 1 0 1 1 1 0 0 ...
    ##  $ baseline    :'data.frame':    300 obs. of  4 variables:
    ##   ..$ id       : num [1:300] 1 2 3 4 6 12 13 17 20 22 ...
    ##   ..$ study    : int [1:300] 1 1 1 1 1 1 1 1 1 1 ...
    ##   ..$ treat    : num [1:300] 0 1 0 1 0 1 1 1 0 0 ...
    ##   ..$ intercept: num [1:300] 1 1 1 1 1 1 1 1 1 1 ...
    ##  $ time.col    : chr "time"
    ##  $ subj.col    : chr "id"
    ##  - attr(*, "class")= chr "jointdata"

We can solve this by running the following code, which changes these variables to factors.

    jointdat$baseline$study <- as.factor(jointdat$baseline$study)
    jointdat$baseline$treat <- as.factor(jointdat$baseline$treat)

Re-checking the structure of the data we can see that these variables are now factors. Care should be taken to ensure that all variables in the `jointdata` object are of the correct format, otherwise they may not be treated correctly during later model fits.

The main dataset we are examining is now in `jointdata` format in object `jointdat`. We will now examine how to remove any longitudinal measurements recorded after the event of interest from the dataset, however the main portion of this document will use `jointdat` for examples, not the objects from the following section.

## Removing any longitudinal information recorded after the event

If a dataset concerns a non-terminal event, it is possible that longitudinal information could be recorded after an individual’s event. For example if the event is time until first seizure, longitudinal information could be available after this first seizure time. However it should not contribute to the joint model as the joint model is investigating the effect of the longitudinal outcome on the risk of an event. As such the `joineRmeta` package contains the function `removeafter` to remove any longitudinal information recorded after each individual’s event time from the `jointdata` object. To see the help file for this function run

    help("removeafter", package = "joineRmeta")

We will demonstrate this function on the simulated dataset `simdat3` with contains longitudinal data not capped at each individual’s event time. If we firstly examine the structure `simdat3` data we notice that it is not yet in `jointdata` format.

    str(simdat3)

To change this dataset into `jointdata` format again we run the `tojointdata` function

    jointdat3<-tojointdata(longitudinal = simdat3$longitudinal,
      survival = simdat3$survival, id = "id", longoutcome = "Y",
      timevarying = c("time","ltime"), survtime = "survtime", cens = "cens",
      time = "time")

Now that we have `simdat3` in `jointdata` format we can run `removeafter` on it to remove any longitudinal information recorded after an individual’s event.

    jointdat3.1<-removeafter(data = jointdat3, longitudinal = "Y",
                             survival = "survtime", id = "id", time = "time")

If we now examine the structures of `jointdat3` and `jointdat3.1` we can see that some of the longitudinal information, namely any that was recorded after an individual’s event time, is present in the `longitudinal` data frame in `jointdat3`, but has been removed from the `longitudinal` data frame in `jointdat3.1`. Note that again there are variables in these datasets (namely `study` and `treat`) that would need to be converted to factors if these datasets were being further used.

    str(jointdat3)

    ## List of 6
    ##  $ subject     : num [1:2500] 1 2 3 4 5 6 7 8 9 10 ...
    ##  $ longitudinal:'data.frame':    12500 obs. of  4 variables:
    ##   ..$ id   : num [1:12500] 1 1 1 1 1 2 2 2 2 2 ...
    ##   ..$ Y    : num [1:12500] 1.5 4.18 6.79 9.39 12.02 ...
    ##   ..$ time : num [1:12500] 0 1 2 3 4 0 1 2 3 4 ...
    ##   ..$ ltime: num [1:12500] 0 1 2 3 4 0 1 2 3 4 ...
    ##  $ survival    :'data.frame':    2500 obs. of  3 variables:
    ##   ..$ id      : num [1:2500] 1 2 3 4 5 6 7 8 9 10 ...
    ##   ..$ survtime: num [1:2500] 5.31 1.17 1.59 5.14 4.6 ...
    ##   ..$ cens    : num [1:2500] 1 0 0 1 1 1 1 1 0 1 ...
    ##  $ baseline    :'data.frame':    2500 obs. of  4 variables:
    ##   ..$ id       : num [1:2500] 1 2 3 4 5 6 7 8 9 10 ...
    ##   ..$ study    : int [1:2500] 1 1 1 1 1 1 1 1 1 1 ...
    ##   ..$ treat    : num [1:2500] 0 1 1 1 0 0 1 1 0 1 ...
    ##   ..$ intercept: num [1:2500] 1 1 1 1 1 1 1 1 1 1 ...
    ##  $ time.col    : chr "time"
    ##  $ subj.col    : chr "id"
    ##  - attr(*, "class")= chr "jointdata"

    str(jointdat3.1)

    ## List of 6
    ##  $ subject     : num [1:2500] 1 2 3 4 5 6 7 8 9 10 ...
    ##  $ longitudinal:'data.frame':    5846 obs. of  4 variables:
    ##   ..$ id   : num [1:5846] 1 1 1 1 1 2 2 3 3 4 ...
    ##   ..$ Y    : num [1:5846] 1.5 4.18 6.79 9.39 12.02 ...
    ##   ..$ time : num [1:5846] 0 1 2 3 4 0 1 0 1 0 ...
    ##   ..$ ltime: num [1:5846] 0 1 2 3 4 0 1 0 1 0 ...
    ##  $ survival    :'data.frame':    2500 obs. of  3 variables:
    ##   ..$ id      : num [1:2500] 1 2 3 4 5 6 7 8 9 10 ...
    ##   ..$ survtime: num [1:2500] 5.31 1.17 1.59 5.14 4.6 ...
    ##   ..$ cens    : num [1:2500] 1 0 0 1 1 1 1 1 0 1 ...
    ##  $ baseline    :'data.frame':    2500 obs. of  4 variables:
    ##   ..$ id       : num [1:2500] 1 2 3 4 5 6 7 8 9 10 ...
    ##   ..$ study    : int [1:2500] 1 1 1 1 1 1 1 1 1 1 ...
    ##   ..$ treat    : num [1:2500] 0 1 1 1 0 0 1 1 0 1 ...
    ##   ..$ intercept: num [1:2500] 1 1 1 1 1 1 1 1 1 1 ...
    ##  $ time.col    : chr "time"
    ##  $ subj.col    : chr "id"
    ##  - attr(*, "class")= chr [1:2] "jointdata" "list"

# Plotting the data

Plotting the longitudinal trajectories, and graphs of survival probabilities, are a useful way to examine both whether there appears to be heterogeneity between studies, and what types of random effects structures might be most appropriate for the longitudinal sub-model.

The `joineRmeta` package contains a plotting function that allows the longitudinal and / or the survival data to be plotted for each study in the dataset.

To plot both the survival and the longitudinal data we use the following code

    sepplots <- jointmetaplot(dataset = jointdat, study = "study", longoutcome = "Y",
                              longtime = "time", survtime = "survtime", cens = "cens",
                              id = "id", smoother = TRUE,
                              studynames = c("Study 1", "Study 2", "Study 3"), type = "Both")

The above code asked the plotting function to provide plots for both the longitudinal data and the survival data, by setting argument `type = "Both"`. If just one type of graph were required then `type` is set to one of `Longitudinal` or `Event`. A smoother was requested for the longitudinal information, but this can be turned off with `smoother = FALSE`. The `studynames` parameter is optional, but allows the graphs for each study to be clearly labelled. The names supplied in this parameter are used in order of the unique study names in the supplied data. The object `sepplots` is then a list of two elements in this case `longplots` and `eventplots`. Individual plots for each study can then be called by name or by their place in the list for example

    sepplots$longplots$`studyplot.Study 3`

    ## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, : pseudoinverse used
    ## at -0.02

    ## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, : neighborhood
    ## radius 2.02

    ## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, : reciprocal
    ## condition number 2.8426e-016

    ## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, : There are other
    ## near singularities as well. 1

    ## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, : pseudoinverse used
    ## at -0.02

    ## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, : neighborhood
    ## radius 1.02

    ## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, : reciprocal
    ## condition number 0

    ## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric, : There are other
    ## near singularities as well. 1

    sepplots$eventplots[[1]]

In the above, warning messages have been produced by including the loess smoother on the longitudinal trajectory plots. We refer issues with the loess smoother to the `ggplot2` helpfiles, and note that the smoother here is used as a guide to the mean trajectory in an exploratory analysis rather than the end result of an analysis.

Other options also exist for the event plot. The event plot survival probability can be plotted for different groups, by providing a grouping variable to `eventby` in the function call. Also confidence intervals for the survival plot can be included by setting `eventconfint = TRUE`.

    sepplots2 <- jointmetaplot(dataset = jointdat, study = "study", longoutcome = "Y",
                               longtime = "time", survtime = "survtime", cens = "cens",
                               id = "id", smoother = TRUE,
                               studynames = c("Study 1", "Study 2", "Study 3"),
                               type = "Both", eventby = "treat")

    sepplots3 <- jointmetaplot(dataset = jointdat, study = "study", longoutcome = "Y",
                               longtime = "time", survtime = "survtime", cens = "cens",
                               id = "id", smoother = TRUE,
                               studynames = c("Study 1", "Study 2", "Study 3"),
                               type = "Event", eventconfint =  TRUE)

    sepplots2$eventplots$`studyplot.Study 3`

    sepplots3$eventplots[[2]]

For more information concerning this plotting function see

    help("jointmetaplot", package = "joineRmeta")

It is also useful to see the plots side by side instead of viewing each separately. The function `jointmetaplotall` provides a way to easily do this. For ease of output, we have wrapped this function call in `suppressWarnings` to prevent printout of warnings relating to the loess smoother inclusion in the longitudinal trajectory plot.

    allplot2 <- suppressWarnings(jointmetaplotall(plotlist = sepplots2, ncol = 2,
                                                  top = "All studies",
                                                  type = "Both"))

The above code has generated arranged grids of all the plots held in the `sepplots2` object, for both the survival and the longitudinal information. To examine these plots use the code

    allplot2$longall
    allplot2$eventsall

    allplot2$longall

    allplot2$eventsall

This prints out the arranged plot for the longitudinal trajectories and for the survival data respectively. For more help on the `jointmetaplotall` function see

    help("jointmetaplotall", package = "joineRmeta")

# Two stage methods

## Pre-analysis - data

When conducting a two stage meta-analysis of joint longitudinal and survival data using the methods available in this package we firstly assume that the data is in joint data format (explanation for putting in this format is above). The data should be examined and plotted, for example using the `jointmetaplot` function. From these plots, and with consideration for the type of data and the research questions of the investigation, a decision should be made as to the most appropriate joint modelling structure for each study.

The association structure used in the joint models fitted to the data from each study should be chosen based on the aims of the investigation. For example, if the investigation wishes to establish whether the current value (fixed and random) effects of the longitudinal trajectory significantly effects the risk of an event, then a current value sharing structure as available in the `JM` package must be employed (see Dimitris Rizopoulos (2012), D. Rizopoulos (2010)). Alternatively if the effect of the rate of change of the longitudinal trajectory on the risk of an event is of interest, then the slope or first derivative of the longitudinal trajectory with respect to time must inserted into the survival sub-model to link the longitudinal and survival components. Again this option is fitted with the `JM` package, which also allows both the current value and the slope of the longitudinal trajectory to be inserted into the survival sub-model each with their own association parameters. Alternatively if the effect of the deviation of an individual from the population mean longitudinal trajectory on the risk of an event is of interest, then an association structure that shares only the random effects between the sub-model should be used (see Gould et al. (2015), Wulfsohn and Tsiatis (1997) and Henderson, Diggle, and Dobson (2000)). Such models can be fitted using the `joineR` package (see Philipson et al. (2012)).

Different studies may require different fixed effects and/or random effects to be included in their joint model. For example the data from one study might display individual variation in the slope of the longitudinal trajectory, whereas the data from another study may not. Additionally different fixed effects might be available in different studies, or might have a significant effect on the response in one study but not another. Dependent on the type of association structure selected for the investigation, the interpretations of parameters can differ between joint model fits depending on their fixed effect and random effect specifications. If the random effects specification differs between the study specific joint model fits, and either a fixed and random effect association structure or random effects only association structure has been selected, then the interpretation of the association parameters will differ between studies. Additionally for fixed and random effect association structures, if the fixed effects included in the longitudinal trajectory differ between the study specific fits, then the association structure will differ between the joint model fits, so again the association parameters between model fits would have different interpretations. Finally if one study fit involves one type of association structure (e.g. the current value of the longitudinal trajectory) and another study fit used a different association structure (e.g. random effects only) then the interpretation of the association parameters again will be different between studies.

To summarise, only parameters with the same interpretation should be pooled in a meta-analysis. Care needs to be taken to ensure that the association parameters have the same interpretation. Differing interpretations could stem from differing random effects specifications, or differing longitudinal fixed effects specifications if a fixed and random effect association structure is used, or from different association structures being employed in different model fits. Special care should be taken for the case where the first derivative of the longitudinal trajectory is inserted into the survival sub-model, as the association structure will only involve the derivative of fixed or random effects that are some function of the longitudinal time variable.

If different joint modelling structures are appropriate for different studies identified in a two stage meta-analysis we recommend grouping studies by identical joint modelling structures and only quantitatively pooling data from those studies within the same group.

## First stage - study specific joint model fits

In the first stage of a two stage meta-analysis of joint longitudinal and survival data identical joint models are fitted to the data from each study. Depending on the type of joint models to be fitted, and the number of studies included in the meta-analysis, these model fits can take some time to fit. For example the models fitted using the `joineR` package require a bootstrapping procedure to estimate the standard errors, which can take several hours (see Hsieh, Tseng, and Wang (2006) for reasons for bootstrapping). For ease of use of this vignette, we provide a range of joint model fits, and where appropriate the bootstrapped standard error estimates, from the `joineR` and the `JM` packages to illustrate the second stage of the meta-analysis. We will briefly introduce these fits and list the code needed to fit them, but we highlight that the fits are available in this package so that the user does not need to take time to bootstrap model fits whilst completing this tutorial. As noted earlier these study specific model fits are fitted using the data held in the `simdat2` dataset.

### Fits using the joineR package

Two sets of example fits from the `joineR` package are available in the `joineRmeta` package. The first is available in the `joineRfits` object, the help file for which can be examined using:

    help("joineRfits", package = "joineRmeta")

This object contains a list of `joineR` fits and the results of the bootstrapping procedure from the package to estimate the standard errors of the model fits. We can extract these fits using the following code:

    joineRmodels <- joineRfits[c("joineRfit1", "joineRfit2", "joineRfit3")]
    joineRmodelsSE <- joineRfits[c("joineRfit1SE", "joineRfit2SE", "joineRfit3SE")]

The result of this code is that the model fits for each study in the `simdat2` data are held in the object `joineRmodels`, and the bootstraps for these model fits are held in the object `joineRmodelsSE`. Each of these model fits has the same structure which can be viewed using the following code:

    summary(joineRmodels[[1]])

    ##
    ## Call:
    ## joint(data = jointdat.1, long.formula = Y ~ 1 + time + treat,
    ##     surv.formula = Surv(survtime, cens) ~ treat, model = "intslope")
    ##
    ## Random effects joint model
    ##  Data: jointdat.1
    ##  Log-likelihood: -484.796
    ##
    ## Longitudinal sub-model fixed effects: Y ~ 1 + time + treat
    ## (Intercept) 0.294367
    ## time        2.707325
    ## treat1      2.055769
    ##
    ## Survival sub-model fixed effects: Surv(survtime, cens) ~ treat
    ## treat1 2.064662
    ##
    ## Latent association:
    ## gamma_0 0.9333956
    ##
    ## Variance components:
    ##        U_0        U_1   Residual
    ## 0.90780689 1.19371405 0.01057928
    ##
    ## Convergence at iteration: 27
    ##
    ## Number of observations: 282
    ## Number of groups: 100

We can see from the above code that a model with a random intercept and random slope formulation has been fitted to the data from each study, and that the longitudinal sub-model contains a fixed intercept, time (slope) and treatment term, and finally that the survival sub-model contains a fixed treatment term.

The second set of example fits from the `joineR` package is held in the object `joineRfits2`. The help file for these fits can be viewed using the following code:

    help("joineRfits2", package = "joineRmeta")

This object contains a list of `joineR` fits and the results of the bootstrapping procedure from the package to estimate the standard errors of the model fits. We can extract these fits using the following code:

    joineRmodels2 <- joineRfits2[c("joineRfit6", "joineRfit7", "joineRfit8",
                                   "joineRfit9", "joineRfit10")]
    joineRmodels2SE <- joineRfits2[c("joineRfit6SE", "joineRfit7SE", "joineRfit8SE",
                                     "joineRfit9SE", "joineRfit10SE")]

The result of this code is that the model fits for each study in the `simdat2` data are held in the object `joineRmodels2`, and the bootstraps for these model fits are held in the object `joineRmodels2SE`. Each of these model fits has the same structure which can be viewed using the following code:

    summary(joineRmodels2[[1]])

    ##
    ## Call:
    ## joint(data = jointdat.1, long.formula = Y ~ 1 + time + treat,
    ##     surv.formula = Surv(survtime, cens) ~ treat, model = "int")
    ##
    ## Random effects joint model
    ##  Data: jointdat.1
    ##  Log-likelihood: -812.2205
    ##
    ## Longitudinal sub-model fixed effects: Y ~ 1 + time + treat
    ## (Intercept) 0.642934
    ## time        2.005585
    ## treat1      2.004105
    ##
    ## Survival sub-model fixed effects: Surv(survtime, cens) ~ treat
    ## treat1 1.806079
    ##
    ## Latent association:
    ## gamma_0 1.10049
    ##
    ## Variance components:
    ##      U_0 Residual
    ## 1.898174 1.331146
    ##
    ## Convergence at iteration: 40
    ##
    ## Number of observations: 282
    ## Number of groups: 100

We can see from the above summary of the joint model fit that this model contained only a random intercept, but the same fixed effects in the longitudinal and the survival sub-models as the fits held in the `joineRmodels` object. These two sets of fits (`joineRfits` and `joineRfits2`) have been supplied so that we can demonstrate what happens in `joineR` fits with differing random effects structures are supplied to the two stage pooling function in the `joineRmeta` package.

### Fits using the JM package

Two example sets of fits from the `JM` package are available in the `joineRmeta` package. The first is available in the `JMfits` object, the help file for which can be examined using:

    help("JMfits", package = "joineRmeta")

This object contains a list of `JM` fits. We do not need to extract the fits from this object, as all the objects in `JMfits` are of type `jointModel`, a joint model fit. The `JM` package does not require bootstrapping to estimate the standard errors. Each of the elements in the `JMfits` object is the result of fitting a model from `JM` to the data from each study in turn from the example dataset `simdat2`. Each of these model fits has the same structure which can be viewed using the following code:

    summary(JMfits[[1]])

    ##
    ## Call:
    ## jointModel(lmeObject = lmefit1, survObject = survfit1, timeVar = "time",
    ##     parameterization = "value", method = "spline-PH-aGH")
    ##
    ## Data Descriptives:
    ## Longitudinal Process     Event Process
    ## Number of Observations: 282  Number of Events: 65 (65%)
    ## Number of Groups: 100
    ##
    ## Joint Model Summary:
    ## Longitudinal Process: Linear mixed-effects model
    ## Event Process: Relative risk model with spline-approximated
    ##      baseline risk function
    ## Parameterization: Time-dependent
    ##
    ##    log.Lik      AIC      BIC
    ##  -301.5223 637.0446 681.3325
    ##
    ## Variance Components:
    ##              StdDev    Corr
    ## (Intercept)  0.9534  (Intr)
    ## time         1.1027  0.3163
    ## Residual     0.1025
    ##
    ## Coefficients:
    ## Longitudinal Process
    ##              Value Std.Err z-value p-value
    ## (Intercept) 0.2735  0.0853  3.2080  0.0013
    ## time        2.7473  0.0369 74.3892 <0.0001
    ## treat1      2.0508  0.1309 15.6620 <0.0001
    ##
    ## Event Process
    ##           Value Std.Err z-value p-value
    ## treat1   0.1941  0.3004  0.6462  0.5182
    ## Assoct   0.9368  0.1028  9.1149 <0.0001
    ## bs1     -6.6541  2.0017 -3.3242  0.0009
    ## bs2     -3.2374  1.0247 -3.1592  0.0016
    ## bs3     -5.2178  0.9645 -5.4099 <0.0001
    ## bs4     -6.4022  0.9332 -6.8604 <0.0001
    ## bs5     -7.7865  1.0686 -7.2864 <0.0001
    ## bs6     -9.7900  1.6444 -5.9535 <0.0001
    ## bs7    -12.9364  2.3964 -5.3982 <0.0001
    ## bs8    -14.4693  2.2515 -6.4265 <0.0001
    ##
    ## Integration:
    ## method: (pseudo) adaptive Gauss-Hermite
    ## quadrature points: 5
    ##
    ## Optimization:
    ## Convergence: 0

From the above summary we can see a range of details about the model fit. The longitudinal sub-model contains a fixed intercept, time (slope) and treatment assignment term. The survival sub-model contains only a fixed treatment assignment term. The baseline hazard was modelled using splines, whose parameters are listed as `bs1` to `bs8` in the results for the Event process. A random intercept and random slope were included in the model. The current value (fixed and random effects) of the longitudinal trajectory was inserted into the survival sub-model. The association parameter for this association structure is denoted by `Assoct`.

The second set of joint model fits from the `JM` package is available in the `JMfits2` object, the help file for which can be examined using:

    help("JMfits2", package = "joineRmeta")

This object contains a list of `JM` fits. Again, we do not need to extract the fits from this object, as all the objects in `JMfits` are of type `jointModel`, a joint model fit. The `JM` package does not require bootstrapping to estimate the standard errors. Each of the elements in the `JMfits2` object is the result of fitting a model from `JM` to the data from each study in turn from the example dataset `simdat2`. Each of these model fits has the same structure which can be viewed using the following code:

    summary(JMfits2[[1]])

    ##
    ## Call:
    ## jointModel(lmeObject = lmefit6, survObject = survfit6, timeVar = "time",
    ##     parameterization = "both", method = "spline-PH-aGH", derivForm = dform)
    ##
    ## Data Descriptives:
    ## Longitudinal Process     Event Process
    ## Number of Observations: 282  Number of Events: 65 (65%)
    ## Number of Groups: 100
    ##
    ## Joint Model Summary:
    ## Longitudinal Process: Linear mixed-effects model
    ## Event Process: Relative risk model with spline-approximated
    ##      baseline risk function
    ## Parameterization: Time-dependent + time-dependent slope
    ##
    ##    log.Lik      AIC      BIC
    ##  -300.9606 639.9212 689.4194
    ##
    ## Variance Components:
    ##              StdDev    Corr
    ## (Intercept)  0.9531  (Intr)
    ## time         1.0939  0.3193
    ## Residual     0.1026
    ##
    ## Coefficients:
    ## Longitudinal Process
    ##              Value Std.Err z-value p-value
    ## (Intercept) 0.2713  0.1111  2.4424  0.0146
    ## time        2.6658  0.0671 39.7521 <0.0001
    ## treat1      2.0687  0.1492 13.8670 <0.0001
    ## time:treat1 0.1842  0.0834  2.2088  0.0272
    ##
    ## Event Process
    ##             Value Std.Err z-value p-value
    ## treat1     0.0836  0.3606  0.2318  0.8167
    ## Assoct     0.9699  0.1483  6.5393 <0.0001
    ## Assoct.s  -0.1466  0.3276 -0.4474  0.6546
    ## bs1       -6.1035  2.1203 -2.8787  0.0040
    ## bs2       -2.7906  1.3288 -2.1001  0.0357
    ## bs3       -4.7828  1.2226 -3.9120  0.0001
    ## bs4       -6.0285  1.0695 -5.6367 <0.0001
    ## bs5       -7.4763  1.0986 -6.8050 <0.0001
    ## bs6       -9.8464  1.6902 -5.8256 <0.0001
    ## bs7      -12.8981  2.4032 -5.3671 <0.0001
    ## bs8      -14.5395  2.3272 -6.2475 <0.0001
    ##
    ## Integration:
    ## method: (pseudo) adaptive Gauss-Hermite
    ## quadrature points: 5
    ##
    ## Optimization:
    ## Convergence: 1

From the above summary we can extract a range of details about the model fit. The longitudinal sub-model contains fixed effects for the intercept, the time (slope) term, the treatment assignment, and an interaction term between the time and the treatment. As before, a random intercept and random time (slope) has been included in the model. The survival sub-model contains a fixed effect for the treatment assignment. The baseline hazard is modelled using splines, the parameters for which are represented by the terms `bs1` to `bs8` in the Event process results. The joint model contains a current value of the longitudinal trajectory and current slope of the longitudinal trajectory association structure, the association parameters for which are denoted by `Assoct` and `Assoct.s` respectively. We have deliberately included a fixed interaction term with the time and the treatment variables to demonstrate the pooling of parameters from a range of joint models of varying complexity.

We should highlight that we have fitted a range of models to the example data to demonstrate the methods in this package. We are not examining or presenting the model which is necessarily the most appropriate for the data.

## Second stage - pooling results

### Pooling results from joineR fits

In the second stage of a two stage meta-analysis of joint longitudinal and survival data, we assume that the study specific joint model fits are available to the user. We then use the `jointmeta2` function from the `joineRmeta` package to pool the model parameters. To access the help file for this function use:

    help("jointmeta2", package = "joineRmeta")

Firstly consider a case of pooling results from joint model fits estimated using the `joineR` package that were supplied in the `joineRfits` object (described above). Code to perform the second stage of this two stage meta-analyis is:

    MAjoineRfits <- jointmeta2(fits = joineRmodels, SE = joineRmodelsSE,
                               longpar = c("time", "treat1"),
                               survpar = "treat1", assoc = TRUE,
                               studynames = c("Study 1","Study 2", "Study 3"))

In the above code, because we are supplying fits from the `joineR` package, we have supplied model fits to the argument `fits`, but have also supplied the results of the bootstrapping procedure to `SE`. We then specify the names of the longitudinal fixed parameters that we want to perform meta-analyses for to argument `longpar`, similarly for the survival sub-model to `survpar`. If we want to perform a meta-analysis for the association parameter(s), we set `assoc = TRUE`, otherwise this is set to `FALSE`. An additional optional parameter is `studynames`. This is an opportunity to supply character strings denoting the names of the studies in the dataset, in the order of the factor levels for the study membership variable. These character strings will be used to label the meta-analysis results, and will appear on any forest plots produced using these results. A list of lists of meta-analyses is returned by this function. The meta-analyses of longitudinal fixed parameters are held in the list labelled `longMA`, the meta-analyses of survival fixed parameters are held in the list labelled `survMA.direct`, and the meta-analyses of association parameters, if requested, is held in `assocMA`.

These lists can each be of length greater than one if multiple meta-analyses of that type were requested. For example the names of the `longMA` returned are:

    names(MAjoineRfits$longMA)

    ## [1] "time"   "treat1"

So the results are labelled by the parameters they concern. If required, forest plots can easily be generated from these results by applying the `forest` function from the `meta` package to each separate meta-analysis:

    library(meta)
    forest(MAjoineRfits$longMA$treat1)

We can see from the above plot that both fixed and random meta-analyses are performed for each requested meta-analysis. This plot can be saved from the R console to file.

The importance of not pooling parameters with different interpretations has been addressed earlier in this vignette. As such, the `jointmeta2` function will throw errors if joint models of different structures are supplied to be pooled. In this example, the random effects specification differs, and so an error is printed. An error would also have been printed if some of the supplied joint model fits had separate association parameters for each random effect, and some supplied joint model fits had a single association parameter common to all the shared random effects:

    MAjoineRfits2 <- jointmeta2(fits = c(joineRmodels[1:3], joineRmodels2[1:2]),
                                SE = c(joineRmodelsSE[1:3],joineRmodels2SE[1:2]),
                                longpar = c("time", "treat1"), survpar = "treat1",
                                assoc = TRUE,
                                studynames = c("Study 1","Study 2", "Study 3"))

    ## Error in jointmeta2(fits = c(joineRmodels[1:3], joineRmodels2[1:2]), SE = c(joineRmodelsSE[1:3], : Some of the joint model fits have differing random
    ##            effects structures

### Pooling results from JM fits

As the models fitted using the `joineR` package share only the random effects between the sub-models, when pooling model parameters, values come from one source only. However for joint models that involve an association structure containing some function of both fixed and random effects, such as those provided in the `JM` package, survival sub-model fixed effects could be decomposed into direct components and indirect components. Direct components stem from parameters being included in the survival sub-model as fixed effects. Indirect components stem from parameters with fixed effects in the survival sub-model also being present in the association structure. The overall estimate of the parameter then comes from the sum of the direct component plus the product of the relevant association parameter and the indirect component (see Ibrahim, Chu, and Chen (2010)). The `jointmeta2` function takes this into consideration, and provides estimates for the direct estimates for fixed effects of interest from the survival sub-model, as well as overall estimates (that contain the combination of the direct and indirect information).

The first set of fits from the `JM` package held in the object `JMfits` could be pooled in the following way:

    MAJMfits <- jointmeta2(fits = JMfits, longpar = c("time", "treat1"),
                           survpar = "treat1", assoc = TRUE,
                           studynames = c("Study 1","Study 2", "Study 3"))

The above function returns a list of lists. The names of the lists included in the list are `longMA`, `survMA.direct`, `survMA.overall`, and `assocMA`, which contain respectively the meta-analyses of the specified longitudinal parameters of interest, the meta-analyses of the direct values of the survival parameters of interst, the meta-analyses of the overall values of the survival parameters of interest, and the meta-analyses of the association parameters. Note that `survMA.overall` does not appear in meta-analyses of `joineR` fits because fixed effects do not appear in the association structure. We can produce forest plots of the meta-analyses in the same way as before.

The second set of fits from the `JM` package held in the object `JMfits2` has an association structure that shares both the current value of the longitudinal trajectory and the first derivative of the longitudinal trajectory with respect to time between the sub-models. The `jointmeta2` function can handle indirect components of a survival parameter of interest stemming from both the current value and the first derivative parts of the association structure, and the code is identical to that above:

    MAJMfits2 <- jointmeta2(fits = JMfits2, longpar = c("time", "treat1"),
                            survpar = "treat1", assoc = TRUE,
                            studynames = c("Study 1","Study 2", "Study 3"))

Again, the above function returns a list of lists namely `longMA`, `survMA.direct`, `survMA.overall`, and `assocMA`, which contain respectively the meta-analyses of the specified longitudinal parameters of interest, the meta-analyses of the direct values of the survival parameters of interst, the meta-analyses of the overall values of the survival parameters of interest, and the meta-analyses of the association parameters.

As with the `joineR` model fits, the fits from the `JM` package will only be pooled by the `jointmeta2` if they have the same specification. For example the following will throw an error because the joint models contain different longitudinal fixed effects and have different association structures.

    MAtest <- jointmeta2(fits = c(JMfits2[1:3], JMfits[1:2]),
                         longpar = c("time", "treat1"),
                         survpar = "treat1", assoc = TRUE,
                         studynames = c("Study 1","Study 2", "Study 3"))

    ## Error in jointmeta2(fits = c(JMfits2[1:3], JMfits[1:2]), longpar = c("time", : Some of the joint model fits have differing association structures

The `jointmeta2` function will also throw errors if fits from both `joineR` and `JM` are supplied to be pooled:

    MAtest <- jointmeta2(fits = c(JMfits2[1:3], joineRfits[1:2]),
                         longpar = c("time", "treat1"),
                         survpar = "treat1", assoc = TRUE,
                         studynames = c("Study 1","Study 2", "Study 3"))

    ## Error in jointmeta2(fits = c(JMfits2[1:3], joineRfits[1:2]), longpar = c("time", : Some of the joint modelling fits are different classes -
    ##          consider subgrouping

# One stage methods

The package also allows one stage joint longitudinal and survival models to be fitted for a single continuous longitudinal outcome and a single possibly censored time-to-event outcome. These models take in data from all the studies in the meta-analysis, and allow a variety of ways for differences between the studies to be investigated. Between study heterogeneity can be accounted for by including study level random effects, by including fixed study membership indicator variables and their interactions between variables that could vary between studies, or in the survival sub-model by stratifying the baseline hazard by study.

The longitudinal sub-model is a linear mixed effects model. The survival sub-model is a Cox proportional hazards model with an unspecified baseline hazard. The random effects in the model are updated iteratively using pseudo-adaptive quadrature (see “Rizopoulos, D.” (2012)). The model is fitted using an EM algorithm. Standard errors are estimated via a bootstrapping method, as Hsieh, Tseng, and Wang (2006) point out that otherwise standard errors could be underestimated as the baseline hazard of the model is unspecified.

The function `jointmeta1` is used to fit these one stage models. To view the help file for this function, run code:

    help("jointmeta1", package = "joineRmeta")

A `jointdata` object should be supplied to the `data` argument of the function (see `jointdata` in package `joineR` for more information). A `formula` object for the longitudinal sub-model needs to be supplied to argument `long.formula`. A `formula` object for the survival sub-model needs to be supplied to argument `surv.formula`.

A range of other arguments need to be specified to the `jointmeta` function call. The argument `gpt` sets the number of quadrature points used in the Gauss Hermite quadrature procedure to estimate the random effects (the same number of quadrature points is used for each level of the random effects). This defaults to 5. The argument `lgpt` is the number of quadrature points used in the estimation of the log-likelihood, and defaults to 7. The `max.it` specifies the maximum number of iterations that the EM algorithm will perform during model fitting, this defaults to 350 although more iterations could be needed for complex datasets. The argument `tol` sets the tolerance level used for checking the convergence in the EM algorithm, and defaults to `0.001`. The name of the variable holding the study membership for individuals in the dataset must be supplied to argument `study.name`. The argument `strat` should be set to `TRUE` or `FALSE` depending on whether the baseline hazard should be stratified by study or not. The arguments `longsep` and `survsep` should be set to `TRUE` if the separate results for the initial longitudinal and survival fits used to generate the starting values for the EM algorithm respectively should be returned or not (this includes the model fit object that would be returned by fitting the longitudinal or time-to-event sub-models as separate longitudinal or time-to-event models using the `lmer` or `coxph` functions from packages `lme4` and `survival` respectively). The argument `bootrun` should be set to `FALSE` if the log-likelihood calculation is required, `TRUE` otherwise - this option exists to speed up the boostrapping of the joint model, during which it is not necessary to calculate the log-likelihood for each bootstrap. If details of the current parameter values are required to be printed at each iteration of the EM algorithm then argument `print.detail` should be set to `TRUE`.

Random effects are included in the model by supplying vectors of characters strings representing the variables to assign random effects to, to arguments `long.rand.ind` and `long.rand.stud` for the individual level and the study level random effects respectively. The one stage multilevel joint model must contain at least one individual level random effect. Due to the complex nature of the model structure, a maximum of three random effects are permitted at each level. Note that the number of random effects could be altered e.g. by assigning random effects to a factor with more than two levels. Note that assigning a random effect to study membership in `long.rand.stud` is including a random intercept at the study level of the model. Currently the function does not support assigning random effects to interaction terms unless they have been pre-calculated in the dataset. These random effects are zero mean, and so a random effect should only be assigned to a variable in the model if a corresponding fixed or population effect is also present in `long.formula`.

If a random intercept is to be included at the individual level, then `long.rand.ind` must contain the string `"int"`. If no random intercept is to be included at the individual level then `long.rand.ind` must contain the string `"noint"`. Other variables to include as random effects at the individual level should be supplied as character strings of the variable names to this argument.

If a random intercept is to be included at the study level, then the name of the variable holding the study membership should be supplied to the argument `long.rand.stud`. If no study level random intercept is to be included in the model then the name of the study membership variable should not be supplied to this argument. If no study level random effects are to be included in the model, then this argument should be set to `NULL` or not specified at all in the function call.

Currently, the only sharing structure available in this model is `"randprop"` namely random effects only proportional association (see Wulfsohn and Tsiatis (1997)). In this association structure, the zero mean random effects are shared between the longitudinal and the survival sub-models, and the random effects at each level share a common association parameter. It is planned to expand the code to allow for additional association structures or sharing structures in the future, namely random effects only association structure with separate association parameters for each random effect (`"randsep"`), current value of the longitudinal trajectory (`"value"`), current first derivative or slope of the longitudinal trajectory (`"slope"`), and inserting both the current value and current slope of the longitudinal trajectory (`"valandslope"`). However currently argument `sharingstrct` must be set to `"randprop"`, choosing any of the other options for this argument currently results in an error message.

We will now examine some examples of fitting these one stage joint models. The package contains example fitted models and bootstrap results `onestage0`, `onestage1`, `onestage2`, `onestage3`, and `onestage4`. These are models applied to the `simdat2` simulated dataset encountered earlier in this document, available in the package. The data was transformed to a `jointdata` object using the `tojointdata` function, as described earlier, to give object `jointdat`, which was then supplied to the function calls. These fits and bootstrap results are provided so that users can quickly examine what format the results of using the functions in this package will look like without having to wait for bootstraps to be completed, for cases with interaction terms, with stratified baseline hazard or with study level random effects. These examples are respectively:

- onestage0 - naive model with individual level random effects (random intercept and slope), but with no accounting for between study heterogeneity (unstratified baseline hazard, no study level random effects, no inclusion of fixed study membership variable or its interactions)
- onestage1 - model with individual level random effects (random intercept and slope), unstratified baseline hazard, between study heterogeneity accounted for using fixed interaction terms with study membership
- onestage2 - model with individual (random intercept and slope) and study level random effects (random intercept and treatment effect), unstratified baseline hazard, between study heterogeneity accounted for using study level random effects
- onestage3 - model with individual level random effects (random intercept and slope), baseline hazard stratified by study, fixed interaction terms with study membership in the longitudinal sub-model, between study heterogeneity accounted for using stratified baseline hazard in the time-to-event sub-model and fixed interaction terms in the longitudinal sub-model
- onestage4 - model with individual (random intercept and slope) and study level random effects (treatment effect), baseline hazard stratified by study, study membership as a fixed effect in the longitudinal sub-model, between study heterogeneity accounted for using the stratified baseline hazard, study level random effects and fixed study membership variable.

Further information about these model fits can be found by calling the relevant help file, for example:

    help("onestage1", package = "joineRmeta")

It is important when fitting a one stage meta-analytic model to joint longitudinal and time-to-event data not to account for between study heterogeneity in a certain aspect of the data multiple times. For example, including study membership as a fixed effect as well as including a study level random intercept in the longitudinal sub-model accounts would account for between study heterogeneity in the model intercept twice. Similarly, stratifying the baseline hazard by study, but also including a study level random intercept in the longitudinal sub-model would account for between study heterogeneity in the time-to-event sub-model twice, due to the presence of the random effects in the time-to-event sub-model through the association structure. As such, care needs to be taken when specifying a model to fit, to ensure that it is reasonable and sensible.

It is possible, but not advisable, to fit a naive model using `jointmeta` that does not take into account any potential differences between studies. The code for fitting such as model is:

    onestagefit0 <- jointmeta1(data = jointdat,
                               long.formula = Y ~ 1 + time + treat,
                               long.rand.ind = c("int", "time"),
                               sharingstrct = "randprop",
                               surv.formula = Surv(survtime, cens) ~ treat,
                               study.name = "study", strat = F)

The results of this model fit and the results of its bootstrap are given in the data `onestage0` in this package.

One method of taking study level differences into account is to include interaction terms between the study membership variable and other variables of interest in the model. In model `onestagefit1` (fitted using the code below) we include interaction terms in both sub-models between the study membership variable and the treatment assignment variable. In this model the baseline hazard is not stratified by study, and we have not specified any random effects at the study level, only a random intercept and slope at the individual level through argument `long.rand.ind`.

    onestagefit1 <- jointmeta1(data = jointdat,
                               long.formula = Y ~ 1 + time + treat*study,
                               long.rand.ind = c("int", "time"),
                               sharingstrct = "randprop",
                               surv.formula = Surv(survtime, cens) ~ treat*study,
                               study.name = "study", strat = F)

The `jointmeta1` function returns an object of class `jointmeta1`, information for which can be found by running:

    help("jointmeta1.object", package = "joineRmeta")

Once the model has been fitted by running function `jointmeta1`, we would run the boostrapping function `jointmetaSE` to estimate the standard errors for the model parameters. The `jointmetaSE` function allows the standard errors for combinations of estimated model parameters to be calculated. For example, for the model fitted in `onestagefit1` we have a treatment effect estimate, and study specific adjustments for this treatment effect estimate. If it would be of interest to know the estimated treatment effect in each study in the dataset, we can calculate the overall value (main estimate plus study specific treatment effect estimate) for each study and its standard error using the bootstrapping function. To bootstrap this model, and estimate the standard errors for the study specific treatment effect estimates where `study1` is the baseline study, we run the following code:

    onestagefit1SE <- jointmetaSE(fitted = onestagefit1, n.boot = 200,
                                  overalleffects = list(long = list(c("treat1", "treat1:study2"),
                                                                    c("treat1", "treat1:study3")),
                                                        surv = list(c("treat1", "treat1:study2"),
                                                                    c("treat1", "treat1:study3"))))

In the above code we can see that 200 bootstraps are requested, and that overall effects are listed for both the longitudinal and survival sub-model as we have interaction terms of interest for both sub-models. If the standard errors of the overall estimates are not of interest, then the `overalleffects` argument can simply be left `NULL`. Additionally, it is not necessary to supply both the `long` and `surv` lists in the `overalleffects` argument, for example if overall effects are only required from the longitudinal sub-model, then the list `long` would be present in the argument `overalleffects` but the list `surv` would not be present or would be set to `NULL`. The `jointmetaSE` function returns an object of class `jointmeta1SE1`, the help file for which can be examined by running:

    help("jointmeta1SE.object", package = "joineRmeta")

In a normal research project we would then have two objects, `onestagefit1` and `onestagefit1SE` in the R workspace. We have saved these two objects and made them available in the package so that this tutorial can be completed without waiting for bootstraps to complete. We will extract these now, but these following lines are not necessary when fitting a one stage model to your own data.

    #extract the saved model fit and bootstrap results
    onestagefit1 <- onestage1$onestagefit1
    onestagefit1SE <- onestage1$onestagefit1SE

Once the one stage model has been fitted, we can start to examine the model fit and the bootstrapped SE results to extract information. The first function of interest is the `summary` function which is applied to objects of class `jointmeta1` - the one stage model fits.

    summary(onestagefit1)

    ##
    ## Call:
    ## jointmeta1(data = jointdat, long.formula = Y ~ 1 + time + treat *
    ##     study, long.rand.ind = c("int", "time"), sharingstrct = "randprop",
    ##     surv.formula = Surv(survtime, cens) ~ treat * study, study.name = "study",
    ##     strat = F)
    ##
    ## Random effects joint meta model
    ##  Data: jointdat
    ##  Log-likelihood: -1802.121
    ##  AIC: 3640.242
    ##
    ## Longitudinal sub-model fixed effects: Y ~ 1 + time + treat * study
    ## (Intercept)    0.26972369
    ## time           2.81195261
    ## treat1         2.02660233
    ## study2         0.60142733
    ## study3         0.03523701
    ## treat1:study2  0.34771000
    ## treat1:study3 -1.69817672
    ##
    ## Time-to-event sub-model fixed effects: Surv(survtime, cens) ~ treat * study
    ## Strat: FALSE
    ##        treat1        study2        study3 treat1:study2 treat1:study3
    ##     2.2928924     0.4181184    -0.1056273     1.1901143     0.5724516
    ##
    ##
    ## Latent association:
    ## gamma_ind_0 1.152764
    ##
    ## Variance components:
    ##               Type        Name     Value
    ## 1 Individual level (Intercept)  1.218564
    ## 2                         time 1.0122324
    ## 3         Residual             0.0199034
    ##
    ## Convergence at iteration: 16
    ##
    ## Number of studies: 3
    ##
    ## Number of individuals per study:
    ##   1   2   3
    ## 100 100 100
    ##
    ## Number of longitudinal observations:
    ##   1   2   3
    ## 282 279 249

This function tells us various information including the log-likelihood of the model, the AIC, model parameter estimates for the fixed effects from each sub-model and the association parameters, the estimates for the variance components (i.e. the random effects and the residual or the measurement error term). The number of iterations taken to converge, the number of studies in the dataset and the number of individuals per study, and the number of longitudinal observations per study are also reported. A print function also exists for `jointmeta1` objects, which prints similar information:

    print(onestagefit1)

    ##
    ## Call:
    ## jointmeta1(data = jointdat, long.formula = Y ~ 1 + time + treat *
    ##     study, long.rand.ind = c("int", "time"), sharingstrct = "randprop",
    ##     surv.formula = Surv(survtime, cens) ~ treat * study, study.name = "study",
    ##     strat = F)
    ##
    ## Random effects joint meta model
    ##  Data: jointdat
    ##
    ## Longitudinal sub-model fixed effects: Y ~ 1 + time + treat * study
    ## (Intercept)    0.26972369
    ## time           2.81195261
    ## treat1         2.02660233
    ## study2         0.60142733
    ## study3         0.03523701
    ## treat1:study2  0.34771000
    ## treat1:study3 -1.69817672
    ##
    ## Time-to-event sub-model fixed effects: Surv(survtime, cens) ~ treat * study
    ## Strat: FALSE
    ##        treat1        study2        study3 treat1:study2 treat1:study3
    ##     2.2928924     0.4181184    -0.1056273     1.1901143     0.5724516
    ##
    ##
    ## Latent association:
    ## gamma_ind_0 1.152764
    ##
    ## Variance components:
    ##               Type        Name     Value
    ## 1 Individual level (Intercept)  1.218564
    ## 2                         time 1.0122324
    ## 3         Residual             0.0199034
    ##
    ## Number of studies: 3
    ##
    ## Number of individuals per study:
    ##   1   2   3
    ## 100 100 100
    ##
    ## Number of longitudinal observations:
    ##   1   2   3
    ## 282 279 249

The fixed effects can be extracted from the model fit using the `fixef` function. The source of the fixed effects should be specified. For example to extract the fixed effects from the longitudinal sub-model set `type = "Longitudinal"`, to extract the fixed effects from the survival sub-model set `type = "Survival"`, and finally to extract the association parameters set `type = "Latent"`:

    fixef(onestagefit1, type = "Longitudinal")

    ##   (Intercept)          time        treat1        study2        study3 treat1:study2 treat1:study3
    ##    0.26972369    2.81195261    2.02660233    0.60142733    0.03523701    0.34771000   -1.69817672

    fixef(onestagefit1, type = "Survival")

    ##        treat1        study2        study3 treat1:study2 treat1:study3
    ##     2.2928924     0.4181184    -0.1056273     1.1901143     0.5724516

    fixef(onestagefit1, type = "Latent")

    ## gamma_ind_0
    ##    1.152764

It is also possible to extract the estimates of the modes of any random effects specified in the model using the function `ranef`. To extract the estimates of the individual level random effects for individuals within each study run the following code:

    ranef(onestagefit1, type = "individual")

This will return a list of length equal to the number of studies in the dataset, each element of the list will be a matrix with number of rows equal to the number of individuals in the dataset, and number of columns equal to the number of random effects at that level. To extract the study level random effects, the argument `type` would be set to `"study"`, however for this particular model fit this would result in an error message being printed, as no study level random effects were included in the model.

Another important piece of information about the random effects in a joint model is their estimated covariance matrix. This can easily be extracted from the model fit using the function `rancov`. If the estimated covariance matrix for the individual level random effect is required then the following code is run:

    rancov(onestagefit1, type = "individual")

    ##             (Intercept)      time
    ## (Intercept)   1.2185640 0.4570225
    ## time          0.4570225 1.0122324

Again, if the estimated covariance matrix for the study level random effects is required then this code is used, but `type` is set to `"study"`. However this would result in an error message in this case because study level random effects were not included in `onestagefit1`.

It also may be of interest to extract formulas for the model fit. This is done using function `formula`. This takes a `jointmeta1` object, and a character string specifying the type of formula required (`"Longitudinal"` to return the formula for the longitudinal sub-model, `"Survival"` to return the formula for the survival sub-model, `"Rand_ind"` to return the formula for the individual level random effects, and `"Rand_stud"` to return the formula for the study level random effects if included in the model). For example:

    formula(onestagefit1, type = "Longitudinal")

    ## Y ~ 1 + time + treat * study

    formula(onestagefit1, type = "Survival")

    ## Surv(survtime, cens) ~ treat * study

    formula(onestagefit1, type = "Rand_ind")

    ## ~1 + time

Several functions also exist to use with `jointmeta1SE` objects. A `jointmeta1SE` object consists of three elements; `results` (a data frame containing the estimates, standard errors and confidence intervals for model parameters), `covmat` (the covariance matrix for the model parameters), and `bootstraps` (a data frame holding the results of each bootstrap conducted). To make it easier to extract some of this information, using the `print` function on a `jointmeta1SE` object just prints out the `results` element of the object:

    print(onestagefit1SE)

    ##               Component                Parameter Estimate     SE 95%Lower 95%Upper
    ## 1          Longitudinal              (Intercept)   0.2697  0.129   0.0236    0.566
    ## 2                                           time    2.812 0.0818   2.6714   2.9745
    ## 3                                         treat1   2.0266 0.1838   1.6452   2.3608
    ## 4                                         study2   0.6014 0.1913   0.1871   0.9397
    ## 5                                         study3   0.0352 0.1825  -0.3539   0.3906
    ## 6                                  treat1:study2   0.3477 0.2659  -0.1841    0.898
    ## 7                                  treat1:study3  -1.6982 0.2687  -2.2538  -1.1355
    ## 8  Longitudinal Overall treat1 and treat1:study2   2.3743 0.1917   2.0053    2.789
    ## 9                       treat1 and treat1:study3   0.3284 0.1786   -0.069   0.6434
    ## 10             Survival                   treat1   2.2929 0.3541   1.5947   2.9743
    ## 11                                        study2   0.4181 0.3759  -0.4712    1.054
    ## 12                                        study3  -0.1056 0.3932  -0.9556   0.5458
    ## 13                                 treat1:study2   1.1901 0.5639   0.1221   2.3117
    ## 14                                 treat1:study3   0.5725 0.5142  -0.3908   1.7028
    ## 15     Survival Overall treat1 and treat1:study2    3.483 0.5219   2.6751   4.5564
    ## 16                      treat1 and treat1:study3   2.8653 0.4383   2.0492   3.8231
    ## 17          Association              gamma_ind_0   1.1528 0.0914   0.9915   1.3675
    ## 18             Variance                     b2_0   1.2186 0.0925   1.0308    1.386
    ## 19                                          b2_1   1.0122  0.118   0.7899   1.2427
    ## 20                                      Residual   0.0199 0.0037   0.0137   0.0276

From this we can see we have the model estimates, standard errors and confidence intervals for the model parameters, and overall effects, from both sub-models, the association parameter and the variance parameters (random effects and residual term). Because the print function exists, simplying typing `onestagefit1SE` into the console would also have had the same result, as would typing `onestagefit1SE$results`. We can also just extract the estimate and the confidence intervals in a data frame from the bootstrapped results using `confint`:

    confint(onestagefit1SE)

    ##               Component                Parameter Estimate 95%Lower 95%Upper
    ## 1          Longitudinal              (Intercept)   0.2697   0.0236    0.566
    ## 2                                           time    2.812   2.6714   2.9745
    ## 3                                         treat1   2.0266   1.6452   2.3608
    ## 4                                         study2   0.6014   0.1871   0.9397
    ## 5                                         study3   0.0352  -0.3539   0.3906
    ## 6                                  treat1:study2   0.3477  -0.1841    0.898
    ## 7                                  treat1:study3  -1.6982  -2.2538  -1.1355
    ## 8  Longitudinal Overall treat1 and treat1:study2   2.3743   2.0053    2.789
    ## 9                       treat1 and treat1:study3   0.3284   -0.069   0.6434
    ## 10             Survival                   treat1   2.2929   1.5947   2.9743
    ## 11                                        study2   0.4181  -0.4712    1.054
    ## 12                                        study3  -0.1056  -0.9556   0.5458
    ## 13                                 treat1:study2   1.1901   0.1221   2.3117
    ## 14                                 treat1:study3   0.5725  -0.3908   1.7028
    ## 15     Survival Overall treat1 and treat1:study2    3.483   2.6751   4.5564
    ## 16                      treat1 and treat1:study3   2.8653   2.0492   3.8231
    ## 17          Association              gamma_ind_0   1.1528   0.9915   1.3675
    ## 18             Variance                     b2_0   1.2186   1.0308    1.386
    ## 19                                          b2_1   1.0122   0.7899   1.2427
    ## 20                                      Residual   0.0199   0.0137   0.0276

Finally, the function `vcov` allows the covariance matrix for the bootstrapped estimates to be extracted:

    vcov(onestagefit1SE)

Alternatively, typing `onestagefit1SE$covmat` into the console would have had the same result.

We have demonstrated the methods for the one stage portion of the `joineRmeta` package using a model without study level random effects. An example of a fit with random effects at both the individual and the study level is held in data `onestage2`, for the help file:

    help("onestage2", package = "joineRmeta")

Again this data is the result of fitting a one stage joint model and then bootstrapping the model fit using the following code to fit the model:

    onestagefit2 <- jointmeta1(data = jointdat,
                               long.formula = Y ~ 1 + time + treat,
                               long.rand.ind = c("int", "time"),
                               long.rand.stud = c("study", "treat"),
                               sharingstrct = "randprop",
                               surv.formula = Surv(survtime, cens) ~ treat,
                               study.name = "study", strat = F)

And to bootstrap the standard errors:

    onestagefit2SE<-jointmetaSE(fitted = onestagefit2, n.boot = 200)

From the above function calls, we can see that this model has included two study level random effects, namely a study level random intercept included in the model by supplying the name of the study membership variable `study` to the argument for the study level randome effects `long.rand.stud` and also a study level random treatment effect `treat` which would quantify variation between studies in true treatment effect. As earlier, we still have a random intercept and time at the individual level through `long.rand.ind = c("int","time")`, and the remainder of the arguments remain the same.

In the bootstrap function, unlike for `onestagefit1`, we no longer have a fixed study membership variable or interactions with a fixed study membership variable in either sub-model, and so we do not specify the `overalleffects` argument in the call to `jointmetaSE`.

Again, the results of this model fit and bootstrapping procedure are available in the package. We can extract the results of this code using the following:

    #extract the saved model fit and bootstrap results
    onestagefit2<-onestage2$onestagefit2
    onestagefit2SE<-onestage2$onestagefit2SE

Again note that these fits and bootstrap results are provided in this way so that this tutorial can be completed without waiting for bootstrap code. In reality, the model would be fitted and the results bootstrapped using the code available in the help file for `onestage2`, which could take several hours.

As the model fit `onestagefit2` contains both study and individual level random effects, we can now extract estimates for modes of the study level random effects conditional on the data and the estimated model parameters:

    ranef(onestagefit2, type = "study")

We can extract the study level random effects covariance matrix:

    rancov(onestagefit2, type = "study")

    ##             (Intercept)    treat1
    ## (Intercept)  0.07677683 0.1740077
    ## treat1       0.17400773 0.7615021

And we can extract the formula for the study level random effects:

    formula(onestagefit2, type = "Rand_stud")

    ## ~1 + treat1

Information for the individual level random effects can be extracted from this model using the same methods as earlier shown in this document. The study level random effects also appear on the summary of the model fits, and the bootstrapping results, the code for extracting which is identical to the examples presented for the case without study level random effects.

The model results held in `onestage3` introduce a baseline hazard stratified by study. The model fit used the following function call:

    onestagefit3 <- jointmeta1(data = jointdat,
                               long.formula = Y ~ 1 + time + treat*study,
                               long.rand.ind = c("int", "time"),
                               sharingstrct = "randprop",
                               surv.formula = Surv(survtime, cens) ~ treat,
                               study.name = "study", strat = TRUE)

Here we can see that again we have a random intercept and slope at the individual level, but as `long.rand.stud` is not included in the function call, we do not include any study level random effects. As argument `strat` is set to `T` or `TRUE`, we have included a baseline hazard stratified by study. The fixed study membership variable, and its interaction with treatment, is again included in the longitudinal sub-model. Study membership is not included as a fixed effect in the time-to-event sub-model, because we have a stratified baseline hazard, and so it’s inclusion would doubly account for between study heterogeneity.

The bootstrapping for this model used the following function call:

    onestagefit3SE <- jointmetaSE(fitted = onestagefit3, n.boot = 200,
                                  overalleffects = list(long = list(c("treat1", "treat1:study2"),
                                                                    c("treat1", "treat1:study3")))))

Here we can see that overall effects were requested for the longitudinal sub-model (where a fixed effect was included for study membership, as well as its interaction with the treatment assignment variable). However no overall effects were requested for the time-to-event sub-model (`surv` was not specified in argument `overalleffects`, it could have also been set to `NULL` in the function call).

As before, the results of the model fit and the bootstrapping function are available in the package, and can be extracted using the following code:

    #extract the saved model fit and bootstrap results
    onestagefit3<-onestage3$onestagefit3
    onestagefit3SE<-onestage3$onestagefit3SE

We can view a summary of the model fit using the following code:

    summary(onestagefit3)

    ##
    ## Call:
    ## jointmeta1(data = jointdat, long.formula = Y ~ 1 + time + treat *
    ##     study, long.rand.ind = c("int", "time"), sharingstrct = "randprop",
    ##     surv.formula = Surv(survtime, cens) ~ treat, study.name = "study",
    ##     strat = T)
    ##
    ## Random effects joint meta model
    ##  Data: jointdat
    ##  Log-likelihood: -1591.281
    ##  AIC: 3210.562
    ##
    ## Longitudinal sub-model fixed effects: Y ~ 1 + time + treat * study
    ## (Intercept)    0.26828090
    ## time           2.81143387
    ## treat1         2.02858806
    ## study2         0.60729680
    ## study3         0.03720349
    ## treat1:study2  0.35006878
    ## treat1:study3 -1.70683174
    ##
    ## Time-to-event sub-model fixed effects: Surv(survtime, cens) ~ treat
    ## Strat: TRUE
    ##   treat1
    ## 2.903544
    ##
    ##
    ## Latent association:
    ## gamma_ind_0 1.145833
    ##
    ## Variance components:
    ##               Type        Name     Value
    ## 1 Individual level (Intercept)  1.216498
    ## 2                         time 1.0122119
    ## 3         Residual             0.0192927
    ##
    ## Convergence at iteration: 19
    ##
    ## Number of studies: 3
    ##
    ## Number of individuals per study:
    ##   1   2   3
    ## 100 100 100
    ##
    ## Number of longitudinal observations:
    ##   1   2   3
    ## 282 279 249

The summary is similar to previous model fit print outs. The fact that baseline hazard is stratified by study can be seen both from the function call at the top of the summary, but also through the information about the time-to-event sub-model, where the message `Strat: TRUE` is printed.

Extraction of fixed effects, random effects, covariance matrices, boostrapping results, confidence intervals etc. is performed in the same way for models with stratified baseline hazard as the models without stratified baseline hazard, using the functions and methods already shown in this document. However, as with some earlier examples, as no study level random effects were included in the model again, attempting to extract study level random effects results from the model fit would result in error messages:

    rancov(fitted = onestagefit3, type = "individual")

    ##             (Intercept)      time
    ## (Intercept)   1.2164980 0.4569782
    ## time          0.4569782 1.0122119

    rancov(fitted = onestagefit3, type = "study")

    ## Error in rancov(fitted = onestagefit3, type = "study"): No study level random effects specified in this model

The final set of example models availabe in the package are saved in file `onestage4`. These models include a baseline hazard stratified by study, study level random effect for the treatment effect, and in the longitudinal sub-model a fixed effect for the study membership variable. The call to fit this model is:

    onestagefit4 <- jointmeta1(data = jointdat, long.formula = Y ~ 1 + time + treat + study,
                               long.rand.ind = c("int", "time"), long.rand.stud = c("treat"),
                               sharingstrct = "randprop", surv.formula = Surv(survtime, cens) ~ treat,
                               study.name = "study", strat = TRUE)

Whilst the call to run the bootstraps to obtain the standard errors for the model fit is:

    onestagefit4SE <- jointmetaSE(fitted = onestagefit4, n.boot = 200)

In this model, as we had no interactions between fixed effects, there was no need to use the `overalleffects` argument in the bootstrapping function. As earlier the results for these model fit and bootstrapping functions are available in the package, although the same results can be obtained by running the above function calls. To extract the precalculated model fit and bootstraps run:

    #extract the saved model fit and bootstrap results
    onestagefit4 <- onestage4$onestagefit4
    onestagefit4SE <- onestage4$onestagefit4SE

Again, as with the results from the data file `onestage3` we can see that this model contains a baseline stratified by study by examining the results of the summary function:

    summary(onestagefit4)

    ##
    ## Call:
    ## jointmeta1(data = jointdat, long.formula = Y ~ 1 + time + treat +
    ##     study, long.rand.ind = c("int", "time"), long.rand.stud = c("treat"),
    ##     sharingstrct = "randprop", surv.formula = Surv(survtime,
    ##         cens) ~ treat, study.name = "study", strat = T)
    ##
    ## Random effects joint meta model
    ##  Data: jointdat
    ##  Log-likelihood: -1465.433
    ##  AIC: 2958.867
    ##
    ## Longitudinal sub-model fixed effects: Y ~ 1 + time + treat + study
    ## (Intercept) 0.281210687
    ## time        2.811866494
    ## treat1      1.579158588
    ## study2      0.599279631
    ## study3      0.004874285
    ##
    ## Time-to-event sub-model fixed effects: Surv(survtime, cens) ~ treat
    ## Strat: TRUE
    ##   treat1
    ## 2.928285
    ##
    ##
    ## Latent association:
    ## gamma_ind_0   1.15085692
    ## gamma_stud_0 -0.01016081
    ##
    ## Variance components:
    ##               Type        Name     Value
    ## 1 Individual level (Intercept) 1.2190563
    ## 2                         time 1.0151515
    ## 3      Study level      treat1 0.7694369
    ## 4         Residual             0.0196275
    ##
    ## Convergence at iteration: 20
    ##
    ## Number of studies: 3
    ##
    ## Number of individuals per study:
    ##   1   2   3
    ## 100 100 100
    ##
    ## Number of longitudinal observations:
    ##   1   2   3
    ## 282 279 249

Again, we can see from the time-to-event fixed results in the output (as well as the function call printed at the top of the output) that the baseline hazard is stratified by study.

Extraction of fixed effects, random effects, covariance matrices, boostrapping results, confidence intervals etc. is performed in the same way for models with stratified baseline hazard as the models without stratified baseline hazard, using the functions and methods already shown in this document. However, as study level random effects were included in the model, we can now extract information about them for example:

    rancov(fitted = onestagefit4, type = "individual")

    ##             (Intercept)      time
    ## (Intercept)   1.2190563 0.4582322
    ## time          0.4582322 1.0151515

    rancov(fitted = onestagefit4, type = "study")

    ##           treat1
    ## treat1 0.7694369

As mentioned earlier, the `jointmeta` function has two arguments `longsep` and `survsep` that, if set to `TRUE`, cause the results of separate longitudinal and time-to-event fits to be returned as well as the joint model fit itself. As well as parameter estimates and the like, the entire model fit object is returned as well. If requested in the original joint model fit, this can be extracted in the following way:

    ###CODE NOT RUN
    #to extract the results from a separate longitudinal model
    fitted$sepests$longests$modelfit

    #to extract the results from a separate survival model
    fitted$sepests$survests$modelfit

This allows separate longitudinal and time-to-event models of the same specifications (apart from association structure) to be easily compared to the full joint model fit.

# Conclusion

This vignette has given a brief introduction to the functions available in the `joineRmeta` package. It assumes a working knowledge of joint modelling, and aims to demonstrate the use of the functions rather than the methodology behind them. Additional information about the functions described here can be found in the help files for this package.

# References

Austin, P. C. 2012. “Generating survival times to simulate Cox proportional hazards models with time-varying covariates.” Statistics in Medicine 31 (29): 3946–58. doi:10.1002/sim.5452.

Bender, R., T. Augustin, and M. Blettner. 2005. “Generating survival times to simulate Cox proportional hazards models.” Statistics in Medicine 24: 1713–23.

Gould, L. A., M. E. Boye, M. J. Crowther, J. G. Ibrahim, G. Quartey, S. Micallef, and F. Y. Bois. 2015. “Joint modeling of survival and longitudinal non-survival data: current methods and issues. Report of the DIA Bayesian joint modeling working group.” Statistics in Medicine 34 (14): 2181–95.

Henderson, R., P. J. Diggle, and A. Dobson. 2000. “Joint modelling of longitudinal measurements and event time data.” Biostatistics 1 (4): 465–80. doi:10.1093/biostatistics/1.4.465.

Hsieh, F., Y. K. Tseng, and J. L. Wang. 2006. “Joint modeling of survival and longitudinal data: Likelihood approach revisited.” Biometrics 62 (4): 1037–43.

Ibrahim, J. G., H. Chu, and L. M. Chen. 2010. “Basic concepts and methods for joint models of longitudinal and survival data.” Journal of Clinical Oncology 28 (16): 2796–2801.

Philipson, P., I. Sousa, P. Diggle, P. Williamson, R. Kolamunnage-Dona, and R. Henderson. 2012. JoineR: Joint Modelling of Repeated Measurements and Time-to-Event Data. https://CRAN.R-project.org/package=joineR.

Rizopoulos, D. 2010. “JM: An R package for the joint modelling of longitudinal and time-to-event data.” Journal of Statistical Software 35 (9): 1–32.

Rizopoulos, Dimitris. 2012. Joint Models for Longitudinal and Time-to-Event Data, with Applications in R. Boca Raton, FL: Chapman & Hall/CRC.

“Rizopoulos, D.” 2012. “Fast fitting of joint models for longitudinal and event time data using a pseudo-adaptive Gaussian quadrature rule.” Computational Statistics and Data Analysis 56: 491–501.

Wulfsohn, M.S., and A.A. Tsiatis. 1997. “A joint model for survival and longitudinal data measured with error.” Biometrics 53 (1): 330–39.