1 / 29

Workshop 2 Frailty models

Workshop 2 Frailty models. Contents. Inference for shared gamma frailty models Theoretical considerations Fitting parametric and semiparametric models More flexible shared gamma frailty models Time-varying covariates Smoothing splines Other approaches for frailty models

zudora
Download Presentation

Workshop 2 Frailty models

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Workshop 2Frailty models

  2. Contents • Inference for shared gamma frailty models • Theoretical considerations • Fitting parametric and semiparametric models • More flexible shared gamma frailty models • Time-varying covariates • Smoothing splines • Other approaches for frailty models • Choice of the frailty density • A Bayesian approach for frailty models

  3. The data set: Time to first insemination • Database of regional Dairy Herd Improvement Association (DHIA) • Milk recording service • Artificial insemination • Select sample • Subset of 2567 cows from 49 dairy farms

  4. Fixed covariates data setinsemfix.dat

  5. Time-varying covariates data setinsemtvc.dat

  6. 1 1 0 . . . 1 1 0 1 1 . . . 1 1 Fitting parametric frailty models (1) • Read the data • insemfix<-read.table("c://insemfix.dat",header=T) • Create vectors heifer<- herd stat timeto insem$heifer 1 1 . . 49 .. 49 53 32 201 . . . 103 38

  7. 1 2 3 … 49 Fitting parametric frailty models (2) • Derive quantities n, Di and e • n<-length(levels(as.factor(herd)) 49 • Di<-aggregate(stat,by=list(herd),FUN=sum)[,2] • e<-sum(Di) 11 23 … 34

  8. Fitting parametric frailty models (3) • Observable likelihood for constant hazard

  9. Fitting parametric frailty models (4) • Calculate observable likelihood for parameters • l=log(p[1]), q=log(p[2]), b=p[3] • Cumulative hazard • Hazard Timeto*exp(heifer*p[3])*exp(p[1]) cumhaz cumhaz • aggregate(cumhaz,by=list(herd),FUN=sum)[,2] stat*log(exp(heifer*p[3])*exp(p[1])) lnhaz lnhaz • aggregate(lnhaz,by=list(herd),FUN=sum)[,2]

  10. Fitting parametric frailty models (5) likelihood.exponential<-function(p){ cumhaz<-(timeto*exp(heifer*p[3]))*exp(p[1]) cumhaz<-aggregate(cumhaz,by=list(herd),FUN=sum)[,2] lnhaz<-stat*log(exp(heifer*p[3])*exp(p[1])) lnhaz<-aggregate(lnhaz,by=list(herd),FUN=sum)[,2] lik<-e*log(exp(p[2]))-n*log(gamma(1/exp(p[2]))) +sum(log(gamma(Di+1/exp(p[2]))))- sum((Di+1/(exp(p[2])))*log(1+cumhaz*exp(p[2])))+sum(Di)+sum(lnhaz) -lik}

  11. Fitting parametric frailty models (6) initial<-c(log(0.015),log(0.03),-0.2) t<-nlm(likelihood.exponential,initial,print.level=2) lambda<-exp(t$estimate[1]) theta<-exp(t$estimate[2]) beta<-t$estimate[3]

  12. Interpretation q (1) • From q, the heterogeneity parameter, to density for median time to event calcm<-function(m){ (log(2)/(theta*lambda))^(1/theta) *(1/m)^(1+1/theta)* (1/gamma(1/theta)) *exp(-log(2)/(theta*lambda*m)) } med<-seq(1,100) fmed<-sapply(med,calcm) plot(med,fmed,type="l",xlab="Median time",ylab="Density")

  13. Interpretation q (2)

  14. HR=exp(-0.24)=0.79 Fitting semiparametric models library(survival) coxph(Surv(timeto,stat)~heifer+frailty(herd)) Call: coxph(formula = Surv(timeto, stat) ~ heifer + frailty(herd)) coef se(coef) se2 Chisq DF p heifer -0.24 0.0432 0.0430 30.9 1.0 2.7e-08 frailty(herd) 205.3 40.9 0.0e+00 Iterations: 10 outer, 23 Newton-Raphson Variance of random effect= 0.123 I-likelihood = -16953 Degrees of freedom for terms= 1.0 40.9 Likelihood ratio test=281 on 41.9 df, p=0 n= 2579

  15. Time-varying covariates data Contribution to the denominator: tij=10 : 2.29 tij=55 : 2.61

  16. Fitting Cox models withtime-varying covariates (1) #Read data insemtvc<-read.table("c://insemtvc.dat",header=T) #Restricted data set insemtvcr<-insemtvc[1:10000,] #Create four column vectors, four different variables herd<-insemtvcr$herd begin<-insemtvcr$begin end<-insemtvcr$end score<-insemtvcr$score ureum<-insemtvcr$ureum coxph(Surv(begin,end,score)~ureum+frailty(herd)) ureumtv.frail<-coxph(Surv(begin,end,score)~ureum+ frailty(herd,eps=0.0000001), outer.max=100) summary(ureumtv.frail)

  17. Fitting Cox models withtime-varying covariates (2) Call: coxph(formula = Surv(begin, end, score) ~ ureum + frailty(herd, eps = 1e-07), outer.max = 100) coef se(coef) se2 Chisq DF p ureum -0.0539 0.0683 0.0674 0.62 1.0 4.3e-01 frailty(herd, eps = 1e-07 104.69 15.5 2.8e-15 Iterations: 38 outer, 51 Newton-Raphson Variance of random effect= 0.815 I-likelihood = -4517.5 Degrees of freedom for terms= 1.0 15.5 Likelihood ratio test=372 on 16.5 df, p=0 n= 10000

  18. Stratified model & time varying var • Semiparametric unadjusted model • Semiparametric model stratified for herd • Semiparametric model with random herd effect

  19. Semiparametric unadjusted model insemtvc<-read.table('c://insemtvc.dat',header=T) ureumtv.unadjust <-coxph(Surv(begin,end,score)~ureum,data=insemtvc) > summary(ureumtv.unadjust) Call: coxph(formula = Surv(begin, end, score) ~ ureum, data = insemtvc) n= 93608 coef exp(coef) se(coef) z p l.95 u .95 ureum -0.0273 0.973 0.0162 -1.68 0.093 0.943 1.00 Rsquare= 0 (max possible= 0.779 ) Likelihood ratio test= 2.83 on 1 df, p=0.0926 Wald test = 2.82 on 1 df, p=0.0929 Score (logrank) test = 2.82 on 1 df, p=0.0929

  20. Semiparametric marginal model > ureumtv.marg<-coxph(Surv(begin,end,score)~ureum+cluster(herd),data=insemtvc) > summary(ureumtv.marg) Call: coxph(formula = Surv(begin, end, score) ~ ureum + cluster(herd), data = insemtvc) n= 93608 coef exp(coef) se(coef) robust se z p lower .95 upper .95 ureum -0.0273 0.973 0.0162 0.0285 -0.957 0.34 0.92 1.03 Rsquare= 0 (max possible= 0.779 ) Likelihood ratio test= 2.83 on 1 df, p=0.0926 Wald test = 0.92 on 1 df, p=0.338 Score (logrank) test = 2.82 on 1 df, p=0.0929, Robust = 0.91 p=0.339

  21. Semiparametric stratified model ureumtv.strat<-coxph(Surv(begin,end,score)~ureum+strata(herd),data=insemtvc) > summary(ureumtv.strat) Call: coxph(formula = Surv(begin, end, score) ~ ureum + strata(herd), data = insemtvc) n= 93608 coef exp(coef) se(coef) z p lower .95 upper .95 ureum -0.0588 0.943 0.0198 -2.97 0.003 0.907 0.98 Rsquare= 0 (max possible= 0.444 ) Likelihood ratio test= 8.86 on 1 df, p=0.00291 Wald test = 8.83 on 1 df, p=0.00296 Score (logrank) test = 8.83 on 1 df, p=0.00296

  22. Semiparametric frailty model > ureumtv.frail<-coxph(Surv(begin,end,score)~ureum+frailty(herd,eps=0.0000001),outer.max=100,data=insemtvc) > summary(ureumtv.frail) Call: coxph(formula = Surv(begin, end, score) ~ ureum + frailty(herd, eps = 1e-07), data = insemtvc, outer.max = 100) n= 93608 coef se(coef) se2 Chisq DF p lower .95 upper .95 ureum -0.0525 0.0187 0.0186 7.85 1 0.0051 0.915 0.984 frailty(herd,...) 1238.50 165 0.0000 Iterations: 35 outer, 60 Newton-Raphson Variance of random effect= 0.334 I-likelihood = -69868.9 Degrees of freedom for terms= 1 165 Rsquare= 0.022 (max possible= 0.779 ) Likelihood ratio test= 2066 on 166 df, p=0 Wald test = 7.85 on 166 df, p=1

  23. Unadjusted versus adjusted #Average ureum concentration in herd avureum.cow<-tapply(insemtvc$ureum,list(insemtvc$cowid),mean) herd.cow<-tapply(insemtvc$herd,list(insemtvc$cowid),mean) avureum.herd<-tapply(avureum.cow,list(herd.cow),mean) #Number of inseminations per herd score.herd<-tapply(insemtvc$score,list(insemtvc$herd),sum) results<-data.frame(frailty=ureumtv.frail$frail,ureum=avureum.herd,score=score.herd) results0<-results[results$score==0,] plot(results$frailty,results$ureum,xlab="Frailty",ylab="Mean ureum concentration") points(results0$frailty,results0$ureum,pch=3)

  24. Unadjusted versus adjustedgraphical

  25. Unadjusted model deleting herds without inseminations > insemtvcr<-insemtvc[insemtvc$herd!=1302100819 & insemtvc$herd!=1302310377 & insemtvc$herd!=3103311094 & insemtvc$herd!=3402102723 & insemtvc$herd!=3403801839 & insemtvc$herd!=4300110030,] > ureumtvcr.unadjust<-coxph(Surv(begin,end,score)~ureum,data=insemtvcr) > summary(ureumtvcr.unadjust) Call:coxph(formula = Surv(begin, end, score) ~ ureum, data = insemtvcr) n= 90232 coef exp(coef) se(coef) z p lower .95 upper .95 ureum -0.0433 0.958 0.0163 -2.65 0.008 0.927 0.989 Rsquare= 0 (max possible= 0.789 ) Likelihood ratio test= 7.04 on 1 df, p=0.00795 Wald test = 7.02 on 1 df, p=0.00804 Score (logrank) test = 7.02 on 1 df, p=0.00805

  26. Fitting Cox models withsmoothing splines (1) #Read data insemtvc<-read.table("c://insemtvc.dat",header=T) #Create four column vectors, four different variables herd<-insemtvc$herd begin<-insemtvc$begin end<-insemtvc$end score<-insemtvc$score ureum<-insemtvc$ureum library(survival) fit<-coxph(Surv(begin,end,score)~pspline(ureum,nterm=3,theta=0.5)) temp<-order(ureum) plot(ureum[temp],predict(fit,type='terms')[temp,1],type='l',xlab='Ureum', ylab='Risk score')

  27. Fitting Cox models withsmoothing splines (2)

  28. Fitting the shared normal frailty model (1) insemfix<-read.table("c://docs//onderwijs//workshopICP2//Data//insemfix.dat", header=T) #Create four column vectors, four different variables herd<-insemfix$herd timeto<-insemfix$timeto-min(insemfix$timeto) stat<-insemfix$stat ureum<-insemfix$ureum heifer<-insemfix$heifer coxph(Surv(timeto,stat)~ureum+frailty(herd,distribution="gaussian"))

  29. Fitting the shared normal frailty model (2) Call: coxph(formula = Surv(timeto, stat) ~ ureum + frailty(herd, distribution = "gaussian")) coef se(coef) se2 Chisq DF p ureum -0.0289 0.0290 0.0286 1 1.0 0.32 frailty(herd, distributio 205 38.8 0.00 Iterations: 6 outer, 14 Newton-Raphson Variance of random effect= 0.0948 Degrees of freedom for terms= 1.0 38.8 Likelihood ratio test=246 on 39.8 df, p=0 n= 2579

More Related