[R] Plot survival curves after coxph() with frailty() random effects terms
I would like to plot a survival curves of a group with different categories after running a Cox model with frailty() random effects terms. I just could display a survival plot of the covariable’s mean. Here an example: library(survival) fit<-coxph(Surv(time, status) ~ sex+ frailty(litter, dist='gamma', method='em'), rats) summary(fit ) suf<-survfit(fit) plot(suf, xscale=1, xlab = "Days", ylab="Survival", conf.int=FALSE) lines(suf[1], lwd=2) Warning message: In `[.survfit`(suf, 1) : survfit object has only a single survival curve #But if I use the next code I get the 2 groups. surv_group <- survfit(Surv(time, status) ~ sex, rats) lines(surv_group[1:2], lwd=2, conf.int=FALSE) However, I am not sure about these 2 curves are well done, appropriate. If any of you could help… -- Andreu Ferrero Gregori [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Can mice() handle crr()? Fine-Gray model, Object has no vcov() method.
Here is an example: # example library(survival) library(mice) library(cmprsk) test1 <- as.data.frame(list(time=c(4,3,1,1,2,2,3,5,2,4,5,1, 4,3,1,1,2,2,3,5,2,4,5,1), status=c(1,1,1,0,2,2,0,0,1,1,2,0, 1,1,1,0,2,2,0,0,1,1,2,0), x=c(0,2,1,1,NA,NA,0,1,1,2,0,1, 0,2,1,1,NA,NA,0,1,1,2,0,1), sex=c(0,0,0,NA,1,1,1,1,NA,1,0,0, 0,0,0,NA,1,1,1,1,NA,1,0,0))) dat <- mice(test1,m=10) #Cox regression: cause 1 models.cox1 <- with(dat,coxph(Surv(time, status==1) ~ x +sex )) summary(pool(models.cox1)) #Cox regression: cause 1 or 2 models.cox <- with(dat,coxph(Surv(time, status==1 | status==2) ~ x +sex )) models.cox summary(pool(models.cox)) # crr() #Fine-Gray model models.FG<- with(dat,crr(ftime=time, fstatus=status, cov1=test1[,c( "x","sex")], failcode=1, cencode=0, variance=TRUE)) summary(pool(models.FG)) #Error in pool(models.FG) : Object has no vcov() method. models.FG -- Andreu Ferrero Gregori [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] frailtypack - Model did not converge
Hey, I am using "frailtyPenal" to fit a general joint model: I got a "Model did not converge" message, and I guess it is because I miss specify something in the command: "Call: frailtyPenal(formula = Surv(Time_final_mes_cor, BD_RE2.Re_IC1_cor) ~ cluster(BD_RE2.e_b_id) + (BD_RE2.X.a_probnp_bnpR) + (BD_RE2.ae_presencia_Cizquierda) + (BD_RE2.CKD_EPI60) + terminal(BD_RE2.death1_cor), formula.terminalEvent = ~(BD_RE2.X.a_probnp_bnpR) + (BD_RE2.ae_presencia_Cizquierda) + (BD_RE2.CKD_EPI60), data = BD_AV, recurrentAG = FALSE, jointGeneral = TRUE, n.knots = 20, kappa = c(1, 1), maxit = 700, LIMlogl = 0.0142) General Joint gamma frailty model for recurrent and a terminal event processes using a Penalized Likelihood on the hazard function Convergence criteria: parameters = 6.81e-09 likelihood = 0.0115 gradient = 1 n= 2473 n recurrent events= 406 n terminal events= 195" Any idea??? Cause PC computing takes hours to "fit/non-fit" this model. Thanks, Andreu Ferrero Gregori [[alternative HTML version deleted]] __ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] glmer() -> corrected AUC optimism by bootstraping technic bootMer() [internal validation of a mixed-effects-model]
>
>
> I would like to do an internal validation of a discriminative ability of a
> mixed effects models.
>
> Here is my scrip:
>
> ###
> bootMer-> boot AUC#
> ###
>
> library(lme4)
> library(lattice)
> data(cbpp)
>
> #fit a model
>
> cbpp$Y<-cbpp$incidence>=1
> glmm<-glmer(Y~period + size + (1|herd), family=binomial, data=cbpp)
> glmm
>
> # funcio: versio 3 - no cal posar endpoint en la funcio
> ##
>
>
>
> AUCFun <- function(fit) {
> library(pROC)
> pred<-predict(fit, type="response")
> AUC<-as.numeric(auc(fit@resp$y, pred))
> }
>
>
> #test
>
> (AUCFun(glmm))
>
> ###run bootMer: AUCFun
>
>
>
> system.time(AUC.boot <- bootMer(glmm,nsim=100,FUN=AUCFun,seed=1982,
> use.u=TRUE,
>type="parametric", parallel="multicore",
> ncpus=2))
>
>
> #...
>
> (boot.ci(AUC.boot, index =c(1,1), type="norm"))
>
> roc(cbpp$Y, predict(glmm, type="response"))
>
>
> #Now it seems more reasonable, bias as "optimism"... but still do not know
> #if I am just doing a AUC with bootstrap CI
> **
>
>
>
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
