310 likes | 680 Views
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
E N D
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
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
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
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
Fitting parametric frailty models (3) • Observable likelihood for constant hazard
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]
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}
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]
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")
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
Time-varying covariates data Contribution to the denominator: tij=10 : 2.29 tij=55 : 2.61
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)
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
Stratified model & time varying var • Semiparametric unadjusted model • Semiparametric model stratified for herd • Semiparametric model with random herd effect
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
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
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
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
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)
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
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')
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"))
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