joineRmeta"

In joineRmeta: Joint Modelling for Meta-Analytic (Multi-Study) Data

options(width = 100)

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")

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 @Austin2012, and @Bender2005). 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 @Bender2005). 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))

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)

str(exampledat1$survival)

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))

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))

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)

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)
str(jointdat3.1)

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`
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 @Rizopoulos2012, @Rizopoulos2010). 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 @Gould2015, @Wulfsohn1997 and @Henderson2000). Such models can be fitted using the joineR package (see @Philipson2012).

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 @Hsieh2006 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]])

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")]
joineRmodels2SE <- joineRfits2[c("joineRfit6SE", "joineRfit7SE", "joineRfit8SE")]

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]])

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]])

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]])

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)

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", 
                                           "Study 4", "Study 5"))

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 @Ibrahim2010). 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", 
                                    "Study 4", "Study 5"))

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", 
                                    "Study 4", "Study 5"))

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 @Rizopoulos2012b). The model is fitted using an EM algorithm. Standard errors are estimated via a bootstrapping method, as @Hsieh2006 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 @Wulfsohn1997). 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)

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)

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")

fixef(onestagefit1, type = "Survival")

fixef(onestagefit1, type = "Latent")

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")

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")

formula(onestagefit1, type = "Survival")

formula(onestagefit1, type = "Rand_ind")

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)

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)

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")

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

formula(onestagefit2, type = "Rand_stud")

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)

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")

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

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)

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")

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

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
modelfit$sepests$longests$modelfit

#to extract the results from a separate survival model
modelfit$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

joineRmeta documentation built on Jan. 24, 2020, 5:10 p.m.