0
votes

I am trying to replicate this plot, here:

enter image description here

Here is the source of this plot, slide 89:

http://www.drizopoulos.com/courses/Int/JMwithR_CEN-ISBS_2017.pdf

The top of the plot is the hazard function over time, whereas the bottom green curve is the fitted linear mixed effects model over time.

I have been able to plot both of these separately, however, cannot seem to combine them using either par(mfrow=c(2,1)) or the gridExtra package (because only one is a ggplot object).

I am using the aids and aids.id datasets (as a part of the JM package) in R.


# Load packages JM and lattice
library("JM")
library("lattice")
library("ggplot2")

#Fit models
lmeFit.aids <- lme(CD4 ~ obstime + obstime:drug,
    random = ~ obstime | patient, data = aids)

coxFit.aids <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)

#Plot longitudinal process
p1<-ggplot(data=aids,aes(x=obstime,y=fitted(lmeFit.aids)))
p1<-p1+geom_smooth(se=FALSE)
p1

#Plot survival process 
library(rms)
p2<-psm(Surv(Time,death)~1,data=aids.id)
survplot(p2,what='hazard')


Thank you!

2
Without your data, there's not much we can suggest without making up our own data that many not adequately represent your structure. The following links have several great suggestions for making a question reproducible, namely in the use of dput(x) to provide unambiguous sample data. Thanks! stackoverflow.com/q/5963269, minimal reproducible example, and stackoverflow.com/tags/r/infor2evans
@r2evans, thanks for your comment. I have updated my question accordingly. I am using the 'aids' and 'aids.id' data as a part of the JM package.statistics123
(I'm a little surprised JM passes CRAN's tests: it does not list any of its imported packages in the DESCRIPTION, something I thought was a strict requirement.)r2evans
I loaded the nlme and survival packages to get lme and coxph, respectively, but I cannot produce the p1 plot, complaining: x has insufficient unique values to support 10 knots: reduce k. I can kludge it with jitter (stackoverflow.com/questions/30208670/…), but that comes nowhere close to either of the plots you have.r2evans
However, you asked about combining ggplot and base graphics, perhaps my answer helps?r2evans

2 Answers

1
votes

Up front, patchwork allows you to combine ggplot2 and base graphics into one plot. Adapted from ?wrap_elements,

library(ggplot2)
library(patchwork)
gg <- ggplot(mtcars, aes(mpg, disp)) + geom_point()
gg / wrap_elements(full = ~ plot(mtcars$mpg, mtcars$disp))

ggplot2 over base graphics

1
votes

I was able to extract the values of the hazard at various time points using the survest() function. Then, I was able to plot this using ggplot, meaning I could use grid.arrange().

est<-survest(p2,,what='hazard')
hazard<-data.frame(time=est$time, hazard=est$surv)