library(joineRML) library(knitr)
Tidiers for objects of class mjoint
have been included in latest release of joineRML
package (0.4.5).
The purpose of these tidiers are described in the introductory vignette to broom
:
The broom package takes the messy output of built-in functions in R, such as
lm
,nls
, ort.test
, and turns them into tidy data frames.
There are three distinct tidiers included with broom
:
tidy
: constructs a data frame that summarises the model estimates;
augment
: add columns to the original data that was modeled;
glance
: construct a concise one-row summary of the model.
These methods are specifically useful when plotting results of a joint model or when comparing several models (e.g. in terms of fit).
We use the sample example from the introductory vignette to joineRML
using the heart valve data.
vignette("joineRML", package = "joineRML") help("heart.valve", package = "joineRML")
We analyse only the records with case-complete data for heart valve gradient (grad
) and left ventricular mass index (lvmi
):
data(heart.valve) hvd <- heart.valve[!is.na(heart.valve$grad) & !is.na(heart.valve$lvmi), ]
Further to that, we only select the first 50 individuals to speed up these examples:
hvd <- hvd[hvd$num <= 50, ]
set.seed(12345) fit <- mjoint( formLongFixed = list( "grad" = log.grad ~ time + sex + hs, "lvmi" = log.lvmi ~ time + sex ), formLongRandom = list( "grad" = ~ 1 | num, "lvmi" = ~ time | num ), formSurv = Surv(fuyrs, status) ~ age, data = list(hvd, hvd), timeVar = "time" )
tidy
methodThe tidy
method returns a tidy dataset with model estimates.
tidy(fit)
By default the tidy
method returns the estimated coefficients for the survival component of the joint model; however, it is possible to extract each component by setting the component
argument:
tidy(fit, component = "longitudinal")
It is also possible to require confidence intervals to be calculated by setting conf.int = TRUE
, and modify the confidence level by setting the conf.level
argument:
tidy(fit, ci = TRUE) tidy(fit, ci = TRUE, conf.level = 0.99)
The standard errors reported by default are based on the empirical information matrix, as in mjoint
. It is of course possible to use bootstrapped standard errors as follows:
bSE <- bootSE(fit, nboot = 100, safe.boot = TRUE, progress = FALSE) tidy(fit, boot_se = bSE, conf.int = TRUE)
The results of this example are not included as it would take too long to run for CRAN.
The tidy
method is useful for custom plotting (e.g. forest plots) of results from joineRML
models, all in a tidy framework:
library(ggplot2) out <- tidy(fit, conf.int = TRUE) ggplot(out, aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high)) + geom_point() + geom_errorbar()
augment
methodThe augment
method returns a dataset with added predictions from the joint model. In particular, population-level and individual-level fitted values and residuals are added to the data frame returned by the method:
preds <- augment(fit) head(preds[, c("num", "log.grad", ".fitted_grad_0", ".fitted_grad_1", ".resid_grad_0", ".resid_grad_1")])
head(preds[, c("num", "log.lvmi", ".fitted_lvmi_0", ".fitted_lvmi_1", ".resid_lvmi_0", ".resid_lvmi_1")])
We can plot the resulting predictions for four distinct individuals:
out <- preds[preds$num %in% c(26, 36, 227, 244), ] ggplot(out, aes(x = time, colour = num)) + geom_line(aes(y = log.grad, linetype = "Measured")) + geom_line(aes(y = .fitted_grad_1, linetype = "Fitted")) + labs(linetype = "Type", colour = "ID", y = "Aortic gradient")
glance
methodThe glance
method allows extracting summary statistics from the joint model:
glance(fit)
This allows comparing competing models easily. Say for instance that we fit a second model with only random intercepts:
set.seed(67890) fit2 <- mjoint( formLongFixed = list( "grad" = log.grad ~ time + sex + hs, "lvmi" = log.lvmi ~ time + sex ), formLongRandom = list( "grad" = ~ 1 | num, "lvmi" = ~ 1 | num ), formSurv = Surv(fuyrs, status) ~ age, data = list(hvd, hvd), timeVar = "time" )
We can go ahead and compare the models in terms of AIC and BIC:
glance(fit) glance(fit2)
Several examples of how to use broom
including more details are available on its introductory vignette:
vignette(topic = "broom", package = "broom")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.