Luc Duchateau Ghent University Belgium Overview Frailty distributions The parametric gamma frailty model The parametric positive stable frailty model The parametric lognormal frailty model Frailty distributions ID: 289592
Download Presentation The PPT/PDF document "Basics of the parametric frailty model" is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
Slide1
Basics of the parametric frailty model
Luc Duchateau
Ghent University, BelgiumSlide2
Overview
Frailty distributions
The parametric gamma frailty model
The parametric positive stable frailty model
The parametric lognormal frailty modelSlide3
Frailty distributions
Power variance function family
Gamma
Inverse Gaussian
Positive stable
General PVF
Compound Poisson
LognormalSlide4
Parametric gamma frailty model
Frailty density function (1)
Two-parameter gamma density
One-parameter gamma densitySlide5
Parametric gamma frailty model
Frailty density function examples
One-parameter gamma densitySlide6
Parametric gamma frailty model
Laplace transform of frailty density
Characteristic function
Moment generating function
Laplace transform for positive r.v.Slide7
Parametric gamma frailty model
Laplace transf. generates moments
Generate
n
th
moment
Use
nth derivative of Laplace transform
Evaluate at s=0Slide8
Parametric gamma frailty model
Gamma Laplace transform
Gamma Laplace transformSlide9
Parametric gamma frailty model
Joint survival function (1)
Joint survival function in conditional model
Now use notation
For cluster with covariates Slide10
Parametric gamma frailty model
Joint survival function (2)
Applied
to
Laplace
transform
of gamma distribution we
obtainSlide11
Parametric gamma frailty model
Population survival function (1)
Integrate
conditional
survival
function
Population
density functionPopulation hazard
function Slide12
Parametric gamma frailty model
Population survival function (2)
Applied
to
gamma
distribution
we have
Population hazard function Slide13
Graphically
#Set parameters
condHR
<-2;Ktau.list<-c(0.05,0.1,0.25,0.5,0.75)
Theta.list
<-2*
Ktau.list
/(1-Ktau.list);
Sij
<-
seq
(0.999,0.001,-0.001);
Fij
<-1-Sij
#Plot
population
/
conditional
hazard
plot(
Fij
,(
Sij
)^(
Theta.list
[1]),
xlab
="
Sx,f
(t)",type="n",
ylab
="
Population
/
conditional
hazard",
axes
=
F,ylim
=c(0,1.7))
box();
axis
(1,at=
seq
(0,1,0.2),
labels
=
seq
(1,0,-0.2),
lwd
=0.5);
axis
(2,lwd=0.5)
lines
(c(Fij,1),c((
Sij
)^(
Theta.list
[1]),0),
lty
=1,lwd=1)
lines
(c(Fij,1),c((
Sij
)^(
Theta.list
[2]),0),
lty
=2,lwd=1)
lines
(c(Fij,1),c((
Sij
)^(
Theta.list
[3]),0),
lty
=3,lwd=1)
lines
(c(Fij,1),c((
Sij
)^(
Theta.list
[4]),0),
lty
=4,lwd=1)
lines
(c(Fij,1),c((
Sij
)^(
Theta.list
[5]),0),
lty
=5,lwd=1)
legend(0,1.75,legend=c(
expression
(paste(
tau
,"=0.05, ",
theta
,"=0.105")),
expression
(paste(
tau
,"=0.10, ",
theta
,"=0.222")),
expression
(paste(
tau
,"=0.25, ",
theta
,"=0.500")),
expression
(paste(
tau
,"=0.50, ",
theta
,"=2.000")),
expression
(paste(
tau
,"=0.75, ",
theta
,"=6.000"))),
ncol
=2,lty=c(1,2,3,4,5))Slide14
Parametric gamma frailty model
Population vs conditional hazard
Slide15
Parametric gamma frailty model
Population hazard ratio
Using population hazard functions
For the gamma frailty distributionSlide16
Graphically
plot(
Fij
,(
Sij
^(-
Theta.list
[1]))/(1/
condHR+Sij
^(-
Theta.list
[1])-1)
,
xlab
="
Sx,f
(t)",type="n",
ylab
="
Population
hazard ratio",
axes
=
F,ylim
=c(1,2.5))
box()
axis
(1,at=
seq
(0,1,0.2),
labels
=
seq
(1,0,-0.2),
lwd
=0.5)
axis
(2,at=
seq
(1,2.5,0.5),
labels
=
seq
(1,2.5,0.5),
srt
=90,lwd=0.5)
lines
(c(Fij,1),c((
Sij
^(-
Theta.list
[1]))/(1/
condHR+Sij
^(-
Theta.list
[1])-1),1),
lty
=1,lwd=1)
lines
(c(Fij,1),c((
Sij
^(-
Theta.list
[2]))/(1/
condHR+Sij
^(-
Theta.list
[2])-1),1),
lty
=2,lwd=1)
lines
(c(Fij,1),c((
Sij
^(-
Theta.list
[3]))/(1/
condHR+Sij
^(-
Theta.list
[3])-1),1),
lty
=3,lwd=1)
lines
(c(Fij,1),c((
Sij
^(-
Theta.list
[4]))/(1/
condHR+Sij
^(-
Theta.list
[4])-1),1),
lty
=4,lwd=1)
lines
(c(Fij,1),c((
Sij
^(-
Theta.list
[5]))/(1/
condHR+Sij
^(-
Theta.list
[5])-1),1),
lty
=5,lwd=1)
legend(0,2.5,legend=c(
expression
(paste(
tau
,"=0.05, ",
theta
,"=0.105")),
expression
(paste(
tau
,"=0.10, ",
theta
,"=0.222")),
expression
(paste(
tau
,"=0.25, ",
theta
,"=0.500")),
expression
(paste(
tau
,"=0.50, ",
theta
,"=2.000")),
expression
(paste(
tau
,"=0.75, ",
theta
,"=6.000"))),
ncol
=2,lty=c(1,2,3,4,5))Slide17
Parametric gamma frailty model
Population hazard ratio example
Slide18
Parametric gamma frailty model
The conditional frailty density (1)
Assuming no covariate information
which corresponds for gamma with
~Gamma( , )Slide19
Parametric gamma frailty model
The conditional frailty density (1)
~Gamma( , )Slide20
Quadruples of correlated event times
Cluster of fixed size 4
Example: Correlated infection times in 4 udder quartersSlide21
Exercise
Fit gamma
frailty
model
with
Weibull
baseline hazard to time to infection data at udder quarter
levelSlide22
R Program gamma frailty model
setwd
("c://docs//onderwijs//survival//Flames//notas//")
udder <-
read.table
("udderinfect.dat", header =
T,skip
=2)library(parfm)
cowid<-as.factor(udder$cowid);timeto<-udder$timekstat<-udder$censor;heifer<-udder$LAKTNR
udder<-data.frame(cowid=cowid,timeto=timeto,stat=stat,heifer=heifer)parfm
(Surv(timeto,stat)~heifer,cluster="cowid",data=udder,frailty="gamma")Slide23
Parametric gamma frailty model
Example – parameter estimates
Udder quarter infection dataSlide24
Population and
con
ditional
hazards
Exercise
Depict
the
population hazard together with the
conditional hazards for frailties equal to the mean, median
and the 25th and 95th percentile of the frailty densitySlide25
Population and
conditional
hazards
R
program
lambda
<-0.838;theta<-1.793;alpha<-1.979;beta<-0.317time<-seq(0,4,0.1)condhaz
<-function(t){frail*alpha*lambda*t^(alpha-1)}marghaz<-function(t){(alpha*lambda*t^(alpha-1))/(1+theta*lambda*t^(alpha))}frail<-1;condhaz.frailm<-sapply(time,condhaz);marghaz.marg<-sapply
(time,marghaz);lowfrail<-qgamma(0.25,shape=1/theta,rate=1/theta);upfrail<-qgamma(0.75,shape=1/theta,rate
=1/theta)frail<-lowfrail;condhaz.fraill<-sapply(time,condhaz)frail<-upfrail;condhaz.frailu<-sapply
(time,condhaz)Slide26
Population and
conditional
hazards
Graph
miny
<-min(
condhaz.frailm,condhaz.fraill,condhaz.frailu)maxy<-max(condhaz.frailm,condhaz.fraill,condhaz.frailu)
par(cex=1.2,mfrow=c(1,2))plot(c(min(time),max(time)),c(miny,maxy),type='n',xlab='Time (year quarters)',ylab='hazard function')
lines(time,condhaz.frailm,lty=1);lines(time,marghaz.marg,lty=1,lwd=3)lines(time,condhaz.fraill,lty=2);lines(time,condhaz.frailu,lty=3)Slide27
Population
and
conditional
hazards
Plot
Udder quarter infection data
Heifer
Multiparous
cowSlide28
Population and
con
ditional
hazard ratio -
Exercise
Depict
the population and conditional hazard ratio as a
function of the poulation survival functionSlide29
Population and
con
ditional
hazard ratio -
R-program
#
Set parameterscondHR<-exp(0.317);theta<-1.793;Sij<-
seq(0.999,0.001,-0.001);Fij<-1-Sijpar(mfrow=c(1,1))#Plot population/conditional hazard ratioplot(Fij,(Sij
^(-theta))/(1/condHR+Sij^(-theta)-1),xlab="Sx,f(t)",ylab="Population hazard ratio",type="n",axes=F,ylim
=c(1,2.5))box()axis(1,at=seq(0,1,0.2),labels=seq(1,0,-0.2),lwd=0.5)axis(2,at=seq(1,2.5,0.5),labels=
seq(1,2.5,0.5),srt=90,lwd=0.5)lines(Fij,(Sij^(-theta))/(1/condHR+Sij^(-theta)-1))segments(0,condHR,1,condHR)Slide30
Population
and
conditional
hazard ratio
- plot
Udder quarter infection data
Multiparous cow versus heiferSlide31
The frailty density
mean
and
variance
time evolution - ExerciseDepict the frailty
density mean and variance time evolutionSlide32
The frailty density
mean
and
variance
time evolution - R-program#Plot E(u)plot(
Fij,(Sij^(theta)),xlab="Sx,f(t)",ylab="Conditional mean",type="l")Slide33
The frailty
density
mean
and
variance time evolution – plotUdder quarter infection dataSlide34
Parametric gamma frailty model
Kendall’s tau
Dependence measures developed for binary data. Take two random clusters
i
,
k
with event times
Position gives also covariate informationKendall’s tau is
or alternativelySlide35
Parametric gamma frailty model
Kendall’s tau for gamma frailties
Kendall’s tau can be expressed in terms of the Laplace transform (without proof)
Using the Laplace transform of the gamma frailty, we obtain
Slide36
The cross ratio
function
a
local
version
Kendall’s tau
We only consider bivariate data such as time to
reconstitutionConsider the bivariate risk set for two pairs and
This bivariate risk set takes values between its
maximal size s (number of clusters) and 2
Slide37
The cross ratio
function
definition
We
then
define
the local
measure as We can now
consider this local dependence measure for different
values of r, where r/s is a proxy for time in terms of
survival for uncensored data Slide38
The cross ratio
function
estimation
Consider
all
pairs
with
particular value r = ra and take ratio of concordant
and discordant pairs Often, we rather take a group of adjacent r
a’s due to low sample size We will work
this out of uncensored data, otherwise we need som further approximations
Slide39
The cross ratio
function
R
programme
#Read data
timetodiag
<-
read.table
("timetodiag.csv",header=T,sep=";")timetodiag<-timetodiag[timetodiag$c2!=0,];t1<- timetodiag$t1;t2<- timetodiag$t2numobs<-
length(t1);limit.low<-(seq(0,10)*10)+1;limit.up<- limit.low+9numpairs<-choose(numobs,2)
res<-cbind(limit.low,limit.up,NA);results<-matrix(NA,nrow=numpairs,ncol
=8) Slide40
The cross ratio
function
R
program
#Put
values
pairwise in
results sectioniter<-0for (i in 1:(numobs-1)){ for
(j in (i+1):numobs){ iter<-iter+1 results[iter,1]<-t1[i]
results[iter,2]<-t2[i] results[iter,3]<-t1[j] results[iter,4]<-t2[j]
}} Slide41
The cross ratio
function
R
program
#
determine
the
size of the risk set
for each pairfor (iter in 1:numpairs){ minval1<-min(results
[iter,1],results[iter,3])minval2<-min(results[iter,2],results[iter,4])temp<-timetodiag
[t1>=minval1 & t2>=minval2,]m<-length(temp$t1)results[iter,6]<-m}
Slide42
The cross ratio
function
R
program
#
determine
the cross ratio
function
for each group of ra valuesfor (i in 1:10){low<-
limit.low[i]up<- limit.up[i]temp<-results[results[,6]>= low &
results[,6]<= up,]conc<-0;discord<-0for (j in 1:length(temp[,1])){ signcomp<-
sign((temp[j,1]-temp[j,3])* (temp[j,2]-temp[j,4])) if (signcomp==1) conc<-conc+1
if (signcomp==-1) discord<-discord+1}res[i,3] <-conc
/
discord
}
Slide43
The cross ratio
function
Plot
and
add
model
based g(r)
resrplot((resr[,1]+ resr[,2])/(2*numobs
),resr[,3],xlim=c(1,0),xlab="Estimated survival function",ylab
="Cross ratio function")theta<-1.793cr<-theta+1segments(0,cr,1,cr)
Slide44
Parametric
gamma
frailty
model
Cross ratio
function
from modelThe cross ratio
function, a local measure:Interpretation: time
to recovery from mastitis Positive experience: constitution at time t
2For positively correlated data, we assume that hazard in numerator
>hazard in denominator Slide45
Parametric gamma frailty model
Cross ratio function example
Cross ratio for gamma density is constant
For the reconstitution data, we have
=0.47
=2.793Slide46
Parametric positive stable (PS) frailty model
The positive stable distribtion
Laplace transform
Infinite mean!Slide47
Parametric PS frailty model
Frailty density function examples
Positive stable density functionsSlide48
Parametric PS frailty model
Joint survival function
Joint survival functionSlide49
Parametric PS frailty model
Marginal likelihood (1)
Example: Udder quarter infections, quadruples, clusters of size 4
Five different types of contributions, according to number of events in cluster
Order subjects, first uncensored (1, …, l)
Contribution of cluster
i
is equal to
=0
>0Slide50
Parametric PS frailty model
Marginal likelihood (2)
Derivatives of Laplace transformsSlide51
Parametric PS frailty model
Marginal likelihood (3)
Marginal likelihood expression cluster
iSlide52
Parametric PS frailty model
Population survival function
Integrate conditional survival function
Population density functionSlide53
Parametric PS frailty model
Population hazard function
Population hazard function
Ratio population/conditional hazard Slide54
Parametric PS frailty model
Population vs conditional hazard
Slide55
Parametric PS frailty model
Population hazard ratio
Using population hazard functions
For the PS frailty distributionSlide56
Parametric PS frailty model
Population hazard ratio example
Slide57
Parametric
PS
frailty
model
R
programmeSlide58
Parametric PS frailty model
Example – parameter estimates
Udder quarter infection data
Cond. HR=
Pop. HR=Slide59
Parametric PS frailty model
The conditional frailty density
Assuming no covariate information, conditional density not PS, still PVF
Slide60
Parametric PS frailty model
Dependence measures
Kendall’s tau is given by
Cross ratio function
=0.47Slide61
Parametric PS frailty model
Dependence measures
Cross ratio function – two dimensional
Slide62
Parametric lognormal frailty model
Introduced by McGilchrist (1993) as
Therefore, for frailty we have
Slide63
Parametric lognormal frailty model
Frailty density function examples
Lognormal density functionsSlide64
Parametric lognormal frailty model
Laplace transform
No explicit expression for Laplace transform … difficult to compare
Maximisation of the likelihood is based on numerical integration of the normally distributed frailtiesSlide65
Parametric lognormal frailty model
Example udder quarter infection (1)
Numerical integration using Gaussian quadrature (nlmixed procedure)
Difficult to compare with previous results as mean of frailty no longer 1
Convert results to density function of median event time Slide66
Parametric frailty model
udder infection: lognormal/gamma