1

I am trying to plot the dose-response relation between pm2.5 exposure and hypertension incidence. I fitted a random effects cox model, in which a random effect was added for the study centers. Then I used the plotHR function to plot the dose-response relation. But I met an error. Below is the example R code I used.

library(survival)
library(coxme) 
library(splines)
library(Greg)

data("eortc")
eortc$age<-rnorm(2323,40,10)

efit1 <- coxph(Surv(y, uncens) ~ ns(age) + trt , eortc)
plotHR(efit1,term = 1,xlim = c(30,70))

efit2 <- coxme(Surv(y, uncens) ~ ns(age) + trt + (1|center), eortc)
plotHR(efit2,term = 1,xlim = c(30,70))

I can plot efit1 using plotHR, but I met an error when plot efit2 in which I added a random effect in the cox model. Anyone knows how to solve this problem? Thanks!

epi
  • 25
  • 3
  • It appears that plotHR is not designed for `coxme`-models. I don't see any "A" or "B" in the data or output. You should use [edit] facilities to clarify. Do NOT use comments for that purpose. – IRTFM Sep 04 '18 at 19:17
  • yeah. But I haven't found other functions which can solve this problem. – epi Sep 04 '18 at 19:21
  • The code errors out when it hits a request for `predict( ..., type="terms")`. Since `predict.coxme` allows only one of `c("lp", "risk")` there can be no success with that type of model using plotHR. – IRTFM Sep 04 '18 at 19:29
  • thank you very much for your reply. I am new to the R software. Do you mean I can not use the predict function for the coxme model? – epi Sep 04 '18 at 19:34
  • That's not at all what I said. Please read the help page: `?predict.coxme`. It doesn't allow the use of the "type" of "prediction" required by `plotHR`. – IRTFM Sep 04 '18 at 19:43
  • Persons using R should realize that R-Core package authors are generally sticklers for statistical precision, and as a result may not always try to emulate a SAS-like-style of "giving you everything you might want". For many years there were plaintive pleas from new users to provide predictions for mixed models. and these were resisted on theoretic grounds. – IRTFM Sep 06 '18 at 22:31

1 Answers1

1

This show what predict.coxme can return for such a model. I'm guessing it won't be totally satisfying, but I claim that it is useful because it demonstrates what mixed models are like "underneath the hood". Each center has a separate prediction but since there is no claim that these have any particular distribution, efforts to aggregate them at their means is considered statistically suspect. This plots the exponentiated linear predictors for the model involving a "spline" fit with no knots (which gives you a bunch of lines). I'm coloring by trt status: black for "0" and red for "1"

plot( eortc$age, exp( predict(efit2)), col=1+(eortc$trt==1) )

efit

The "average" difference between the trt==0 and trt==1 groups does show up and is consistent with the measured treatment effect of exp(0.705) -> [1] 2.023847, and the "age" effect which in the model was not significant is a very shallow linear rise.

IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • Thank you very much for your answer. But why did I still cannot run the code : plot( eortc$age, exp( predict(efit2)), col=1+(eortc$trt==1) )? – epi Sep 04 '18 at 20:29
  • I have no idea. Generally, people will need to present full text of any error messages to get an informed response. – IRTFM Sep 04 '18 at 21:03
  • Thank you for reply. Could you please show me the code for your "efit2"? When I am using my "efit2" model, R gives an error because "predict" function cannot be used for coxme. – epi Sep 05 '18 at 00:30
  • There is a `predict.coxme` function in package 'coxme', so I'm unable to understand unless you have a really old installation? – IRTFM Sep 05 '18 at 15:56
  • Looked at `news(pac="coxme")`. Two instance of hte word appearred: 1) "Changes in version 2.2-8 Add predict.coxme to the NAMESPACE file;: and 2) "Changes in version 2.2: Fix two bugs in the return values, caused by incorrect mapping of the b coefficients into the final output, when there were multiple random terms in the model. (The fit was correct. One bug the returned ranef(fit) to be incorrectly ordered, the other led to NA in the linear.predictor.)" Perhaps the reason `predict.coxme` could not be found is that you were using version below 2.2-8. – IRTFM Sep 06 '18 at 22:25