1 / 25

Workshop frailty models

Workshop frailty models. Luc Duchateau, Rosemary Nguti and Paul Janssen. Contents. The statistical package R The data set: time to first insemination Inference for shared gamma frailty models Theoretical considerations Fitting parametric and semiparametric models

candice-bowers
Download Presentation

Workshop 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 frailty models Luc Duchateau, Rosemary Nguti and Paul Janssen

  2. Contents • The statistical package R • The data set: time to first insemination • 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 statistical package R • Freeware (unlike Splus) • Powerful statistical package • Data handling a bit less powerful • Downloadable from the Internet • Base • Libraries • From within R • Download library yourself

  4. Installing R via Internet • Go to www.r-project.org • Choose option download • Choose closest mirrorsite http://cran.za.r-project.org/ • Choose operating system Windows (95 and later) • Choose subdirectory base • Choose rw1090.exe

  5. Installing additional packages • Automatically within R • Menu item « Packages» • « Install Packages from CRAN … » • Choose « survival » • New session, specify « library(survival) » • Download zipped package from CRAN • Unzip it under directory « c:\Program Files\R\rw1091\library »

  6. 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

  7. Fixed covariates data setinsemfix.dat

  8. Time-varying covariates data setinsemtvc.dat

  9. 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

  10. 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

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

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

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

  14. Fitting parametric frailty models (6) initial<-c(100,17,0) t<-nlm(likelihood.exponential,initial,print.level=2) h<-1/t$estimate[1] h theta<-1/t$estimate[2] theta beta<-t$estimate[3] beta

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

  16. Interpretation q (2)

  17. 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

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

  19. Fitting Cox models withtime-varying covariates (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 stat<-insemtvc$stat ureum<-insemtvc$ureum heifer<-insemtvc$heifer library(survival) coxph(Surv(begin,end,stat)~heifer+ureum+frailty(herd))

  20. HR=exp(-0.00349)=0.997 Call: coxph(formula = Surv(begin, end, stat) ~ heifer + ureum + frailty(herd)) coef se(coef) se2 Chisq DF p heifer -0.22820 0.0466 0.0464 23.99 1.0 9.7e-07 ureum -0.00349 0.0366 0.0358 0.01 1.0 9.2e-01 frailty(herd) 177.64 39.4 0.0e+00 Iterations: 10 outer, 23 Newton-Raphson Variance of random effect= 0.116 I-likelihood = -14401.2 Degrees of freedom for terms= 1.0 1.0 39.4 Likelihood ratio test=251 on 41.3 df, p=0 n= 23076 Fitting Cox models withtime-varying covariates (2)

  21. 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 stat<-insemtvc$stat ureum<-insemtvc$ureum heifer<-insemtvc$heifer library(survival) fit<-coxph(Surv(begin,end,stat)~heifer+pspline(ureum,nterm=3,theta=0.5)) temp<-order(ureum) plot(ureum[temp],predict(fit,type='terms')[temp,2],type='l',xlab='Ureum', ylab='Risk score')

  22. Fitting Cox models withsmoothing splines (2) Call: coxph(formula = Surv(begin, end, stat) ~ heifer + pspline(ureum, nterm = 3, theta = 0.5)) coef se(coef) se2 Chisq DF p heifer -0.2284 0.0457 0.0457 24.97 1.00 5.8e-07 pspline(ureum, nterm = 3, 0.0094 0.0379 0.0349 0.06 1.00 8.0e-01 pspline(ureum, nterm = 3, 3.31 1.18 8.7e-02 Iterations: 1 outer, 3 Newton-Raphson Theta= 0.5 Degrees of freedom for terms= 1.0 2.2 Likelihood ratio test=30.2 on 3.18 df, p=1.56e-06 n= 23076

  23. Fitting Cox models withsmoothing splines (3)

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

  25. Fitting the shared normal frailty model (2) Call: coxph(formula = Surv(timeto, stat) ~ heifer + frailty(herd, distribution = "gaussian")) coef se(coef) se2 Chisq DF p heifer -0.241 0.0431 0.043 31.1 1.0 2.4e-08 frailty(herd, distributio 199.9 38.6 0.0e+00 Iterations: 6 outer, 14 Newton-Raphson Variance of random effect= 0.0916 Degrees of freedom for terms= 1.0 38.6 Likelihood ratio test=276 on 39.5 df, p=0 n= 2579

More Related