Vidakovic B. Statistics for Bioengineering Sciences: With Matlab and WinBugs Support

Подождите немного. Документ загружается.

18.5 Bayesian Approach 721

}

for(i in 1:n2) {

ImmPeroxPos[i] ~ dexp(lam2) I(CensorIPP[i], )

}

lam1 ~ dgamma(0.001, 0.001)

lam2 ~ dgamma(0.001, 0.001)

}

DATA

list( n1 = 36, n2 = 9,

ImmPeroxNeg=c(19, 25, 30, 34, 37, 46, 47, 51,

56, 57, 61, 66, 67, 74, 78, 86,

NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,

NA, NA, NA, NA, NA, NA, NA, NA, NA, NA),

CensorIPN = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,

122, 123, 130, 130, 133, 134, 136, 141, 143, 148,

151, 152, 153, 154, 156, 162, 164, 165, 182, 189),

ImmPeroxPos = c(22, 23, 38, 42, 73, 77, 89, 115, NA),

CensorIPP= c(0,0,0,0,0,0,0,0,144))

INITS

list(lam1=1, lam2 = 1,

ImmPeroxNeg=c(NA, NA, NA, NA, NA, NA, NA, NA,

NA, NA, NA, NA, NA, NA, NA, NA,

200, 200, 200, 200, 200, 200, 200, 200, 200, 200,

200, 200, 200, 200, 200, 200, 200, 200, 200, 200),

ImmPeroxPos = c(NA, NA, NA, NA, NA, NA, NA, NA, 200) )

mean sd MCrror val2.5pc median val97.5pc start sample

ImmPeroxNeg[17] 374.8 270.8 0.4261 128.0 289.7 1107.0 1001 500000

ImmPeroxNeg[18] 376.2 270.6 0.4568 129.0 291.2 1107.0 1001 500000

ImmPeroxNeg[19] 383.2 270.6 0.4342 136.0 298.0 1115.0 1001 500000

···

ImmPeroxNeg[35] 435.1 270.2 0.4319 188.1 349.9 1166.0 1001 500000

ImmPeroxNeg[36] 442.3 270.6 0.4272 195.0 357.1 1173.0 1001 500000

ImmPeroxPos[9] 233.1 102.8 0.1730 146.0 200.5 508.2 1001 500000

lam1 0.0042 0.0010 2.849E-6 0.0024 0.0041 0.0065 1001 500000

lam2 0.0128 0.0045 7.496E-6 0.0055 0.0123 0.0231 1001 500000

The Bayes estimator of λ

1,B

= 0.004211, and a 95% credible set is

[0.002404, 0.006505]. The Bayes estimator and interval are close to classical

(Exercise 18.2).

Example 18.13. Smoking Cessation Experiment. The data set for this ex-

ample comes from a clinical trial discussed in Banerjee and Carlin (2004). A

number of smokers entered into a smoking cessation study, and 263 of them

quit. These 263 quitters were monitored and checked to see if and when they

relapsed.

RelapseT is the time to relapse; it is either observed or censored,

722 18 Inference for Censored Data and Survival Analysis

with the censoring indicator contained in the vector censored.time. The inde-

pendent covariates are

Age (age of the individual), AgeStart (age when he/she

started smoking),

SexF (Female=1, Male=0), SIUC (whether the individual re-

ceived an intervention or not), and

F10Cigs (the average number of cigarettes

smoked per day).

A logistic distribution is constrained to a nonnegative domain to model

RelapseT. The parameters of dlogis(mu,tau) are the mean mu, which depends

on the linear combination of covariates, and

tau, which is a rate parameter.

The standard deviation is

π/(

3τ) ≈1.8138/τ.

model {

for (i in 1:N)

{

RelapseT[i] ~ dlogis(mu[i],tau) I(censored.time[i],)

mu[i] <- beta[1] + beta[2]

Age[i] + beta[3]

AgeStart[i] +

beta[4]

SexF[i] + beta[5]

SIUC[i] + beta[6]

F10Cigs[i]

}

for( j in 1:6){

beta[j] ~ dnorm(0, 0.01)

}

tau ~ dgamma(1,0.01)

meanT <- mean(mu[])

sigma <- 1.8138/tau #1.8138 ~ pi/sqrt(3)

# Evaluate Survival Curve for a Subject with covariates:

Ag <- 50; AgSt <- 18; SxF <- 0; S <- 1; Cigs <- 20;

fmu <- beta[1] + beta[2]

Ag + beta[3]

AgSt +

beta[4]

SxF + beta[5]

S + beta[6]

Cigs

for(i in 1:100) {

time[i] <- i/10

Surv[i] <- 1/(1 + exp(tau

(time[i] - fmu)))

}

# Data and Inits omitted (see Smoking.odc)

mean sd MCrror val2.5pc median val97.5pc start sample

beta[1] 2.686 2.289 0.1254 –2.026 2.701 7.284 1001 100000

beta[2] 0.07817 0.03844 0.002078 0.009544 0.07815 0.1509 1001 100000

beta[3] –0.07785 0.07222 0.003745 –0.2181 –0.07965 0.06802 1001 100000

beta[4] –0.9555 0.5382 0.009515 –2.02 –0.9494 0.08981 1001 100000

beta[5] 1.666 0.5993 0.01482 0.4886 1.672 2.817 1001 100000

beta[6] –0.02347 0.02459 9.21E-4 –0.07072 –0.02394 0.0281 1001 100000

meanT 5.314 0.3382 0.00604 4.712 5.292 6.034 1001 100000

sigma 3.445 0.345 0.005367 2.841 3.421 4.187 1001 100000

tau 0.5317 0.05256 8.15E-4 0.4332 0.5302 0.6384 1001 100000

Surv[1] 0.9628 0.01129 5.199E-4 0.9366 0.9642 0.9808 1001 100000

Surv[2] 0.9609 0.01173 5.459E-4 0.9336 0.9623 0.9797 1001 100000

···

Surv[99] 0.1402 0.05249 0.003402 0.06123 0.1334 0.2645 1001 100000

Surv[100] 0.1342 0.05118 0.003311 0.05765 0.1274 0.2557 1001 100000

18.5 Bayesian Approach 723

The Surv values (ordinate, mean, std, median, and quantiles) are exported

to MATLAB as data ﬁle

smokingoutbugs.mat. The ﬁle smokingbugs.m reads

in the data and plots the posterior estimator of the survival curve (Fig. 18.6).

Note that the survival curve is S(t

|µ,τ) = 1/(1 +exp{τ(t −µ)}). The posterior

distribution of S(t

|µ,τ) is understood as a distribution of a function of µ and τ

for t ﬁxed.

2 4 6 8 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Time

S(t)

Fig. 18.6 Bayesian estimator of survival curve. The green bands are the 0.025 and 0.975

percentiles of the posterior distribution for S(t), while the blue errorbars have a size of

posterior standard deviations.

Example 18.14. Duration of Remissions in Acute Leukemia. A data set

analyzed by Freireich et al. (1963), and subsequently by many authors, comes

from a trial of 42 leukemia patients (under age 20) treated in 11 US hospitals.

The effect of 6-mercaptopurine (6-MP) therapy on the duration of remis-

sions induced by adrenal corticosteroids has been studied as a model for test-

ing new agents. Some patients were treated with 6-MP and the rest were con-

trols. The trial was designed as matched pairs. The matching was done with

respect to remission status (partial = 1, complete = 2). Randomization to 6-MP

or control arms was done within a pair. Patients were followed until leukemia

relapsed or until the end of the study. Overall survival was not signiﬁcantly

different for the two treatment programs since patients maintained on placebo

were treated with 6-MP when relapse occurred.

# Status Contr 6-MP # Status Contr 6-MP # Status Contr 6-MP

1 1 1 10 8 2 11 34+ 15 2 8 17+

2 2 22 7 9 2 8 32+ 16 1 23 35+

3 2 3 32+ 10 2 12 25+ 17 1 5 6

4 2 12 23 11 2 2 11+ 18 2 11 13

5 2 8 22 12 1 5 20+ 19 2 4 9+

6 1 17 6 13 2 4 19+ 20 2 1 6+

7 2 2 16 14 2 15 6 21 2 8 10+

724 18 Inference for Censored Data and Survival Analysis

#Duration of Steroid-induced Remissions in Acute Leukemia

model {

for (i in 1:n) {

log(mu[i]) <- b0+ b1[treat[i]] + b2[status[i]]

t[i] ~ dweib(r,mu[i]) I(t.cen[i],)

S[i] <- exp(-mu[i]

pow(t[i],r));

f[i] <- mu[i]

pow(t[i],r-1)

S[i]

Lik[i] <- pow(f[i],1-delta[i])

pow(S[i],delta[i]);

logLik[i] <- log(Lik[i])

}

b0 ~ dnorm(0,0.00001)

b1[1] <- 0

b1[2] ~ dnorm(0,0.001)

b2[1] <- 0

b2[2] ~ dnorm(0,0.001)

r ~ dgamma(0.01,0.01)

Dev <- -2

sum(logLik[]) #deviance

}

DATA

list(n=42,

t = c(1, 22, 3, 12, 8, 17, 2, 11, 8, 12,

2, 5, 4, 15, 8, 23, 5, 11, 4, 1, 8,

10, 7, NA, 23, 22, 6, 16, NA, NA, NA,

NA, NA, NA, 6, NA, NA, 6, 13, NA, NA, NA),

t.cen = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,0,32,0,0,0,0,

34,32,25,11,20,19,0,17,35,0,0,9,6,10),

treat = c(1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1,1,

2,2,2,2,2, 2,2,2,2,2, 2,2,2,2,2, 2,2,2,2,2,2),

status = c(1,2,2,2,2,1,2,2,2,2,2,1,2,2,2,1,1,2,2,2,2,

1,2,2,2,2,1,2,2,2,2,2,1,2,2,2,1,1,2,2,2,2),

delta = c(0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0,0,

0,0,1,0,0, 0,0,1,1,1, 1,1,1,0,1, 1,0,0,1,1,1) )

INITS

list(r=1,b0 = 0, b1=c(NA,0),b2=c(NA,0),

t = c(NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,

NA,NA,NA,NA,NA,NA,NA,NA,32,NA,NA,NA,NA,34,32,

25,11,20,19,NA,17,35,NA,NA,9,6,10))

18.5 Bayesian Approach 725

mean sd MCrror val2.5pc median val97.5pc start sample

Dev 241.2 7.527 0.02825 228.5 240.6 257.8 1001 200000

b0 –3.209 0.6867 0.009883 –4.636 –3.181 –1.949 1001 200000

b1[2] –1.772 0.4241 0.003255 –2.629 –1.76 –0.9696 1001 200000

b2[2] 0.1084 0.4306 0.002543 –0.6967 0.09381 0.9995 1001 200000

r 1.37 0.2069 0.003338 0.9877 1.361 1.798 1001 200000

t[24] 57.34 29.86 0.1341 32.62 48.55 132.9 1001 200000

t[29] 58.96 29.15 0.1345 34.6 50.17 134.7 1001 200000

···

t[41] 37.21 30.96 0.1128 7.161 29.06 115.5 1001 200000

t[42] 39.68 31.09 0.1177 10.95 31.45 117.3 1001 200000

Example 18.15. Photocarcinogenicity. Grieve (1987) and Dellaportas and

Smith (1993) explored photocarcinogenicity in four treatment groups with 20

rats in each treatment group. Treatment 3 is the test drug and the others are

some type of control. Response is time to death or censoring time.

We will ﬁnd Bayes estimators for median survival times for the four treat-

ments. The survival times are modeled as Weibull

W ei(r,µ

), i =1,... , 4. The

WinBUGS code speciﬁes the censoring and priors on the Weibull model. The

data part provides observed and censored times. Posterior densities for me-

dian survival times are given in Fig. 18.7.

model{

for(i in 1 : M) {

for(j in 1 : N) {

t[i, j] ~ dweib(r, mu[i])I(t.cen[i, j],)

}

mu[i] <- exp(beta[i])

beta[i] ~ dnorm(0.0, 0.001)

median[i] <- pow(log(2)

exp(-beta[i]), 1/r)

}

r ~ dexp(0.001)

}

DATA

list( t = structure(.Data = c(12, 1, 21, 25, 11, 26, 27, 30, 13,

12, 21, 20, 23, 25, 23, 29, 35, NA, 31, 36, 32, 27, 23, 12,

18, NA, NA, 38, 29, 30, NA, 32, NA, NA, NA, NA, 25, 30,

37, 27, 22, 26, NA, 28, 19, 15, 12, 35, 35, 10, 22, 18, NA,

12, NA, NA, 31, 24, 37, 29, 27, 18, 22, 13, 18, 29, 28, NA,

16, 22, 26, 19, NA, NA, 17, 28, 26, 12, 17, 26), .Dim = c(4,20)),

t.cen = structure(.Data = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

0, 0, 0, 0, 0, 40, 0, 0, 0, 0, 0, 0, 0, 40, 40, 0, 0, 0, 40, 0, 40,

40, 40, 40, 0, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 24, 0,

40, 40, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 20, 0, 0, 0, 0, 29, 10, 0,

0, 0, 0, 0, 0), .Dim = c(4, 20)), M = 4, N = 20)

INITS

726 18 Inference for Censored Data and Survival Analysis

list( r=1, beta=c(0,0,0,0) ) #generate the rest

mean sd MCrror val2.5pc median val97.5pc start sample

median[1] 23.82 1.977 0.02301 20.20 23.74 27.96 1001 100000

median[2] 35.11 3.459 0.01650 29.23 34.79 42.82 1001 100000

median[3] 26.80 2.401 0.02169 22.46 26.66 31.94 1001 100000

median[4] 21.38 1.845 0.01446 18.09 21.26 25.34 1001 100000

Fig. 18.7 Posterior densities of median survival times median[1]–median[4].

18.6 Exercises

18.1. Simulation of Censoring.

y = exprnd(10,50,1); % Random failure times exponential(10)

d = exprnd(20,50,1); % Drop-out times exponential(20)

t = min(y,d); % Observe the minimum of these times

censored = (y>d); % Observe whether the subject failed

Using MATLAB’s ecdf calculate and plot empirical CDF and conﬁdence

bounds for arguments

t and censored.

18.2. Immunoperoxidase. In the context of Example 18.5 ﬁnd a conﬁdence in-

terval for

, the rate parameter for immunoperoxidase-negative patients.

Use the fact that for MLE

−λ

λ/

(log

−log λ

18.6 Exercises 727

both have an approximately standard normal distribution.

18.3. Massachusetts Data. This exercise will illustrate the survival analy-

sis of outcome predictions from DNA microarrays. Such survival analy-

ses are often found in DNA microarray papers that concern the diagnosis

or prognosis of human cancers. A description is given in Bhattacharjee et

al. (2001) which is available at:

http://www.pnas.org/content/98/24/

13790.full

A publicly available data set of 125 adenocarcinomas can be downloaded

from:

http://www.genome.wi.mit.edu/MPR/lung.

18.4. Expected Lifetime. Let T be a lifetime with survival function S(t). Using

integration by parts in the deﬁnition of

ET show

ET =

∞

S(t)dt.

18.5. Censored Rayleigh. The lifetime (in hours) of a certain sensor has

Rayleigh distribution, with survival function

S(t)

=exp

−

λt

, λ >0.

Twelve sensors are placed under test for 100 hours, and the following fail-

ure times are recorded 23, 40, 41, 67, 69, 72, 84, 84, 88, 100+, 100+. Here +

denotes a censored time.

(a) If failure times t

,... , t

are observed, and t

r+1

,... , t

are censored, show

that the MLE of

λ is

λ =

=r+1

Evaluate the MLE for the given data. Consult Example 18.4.

(b) Calculate and plot Kaplan-Meier estimator and superimpose S(t) eval-

uated at

λ.

18.6. MLE for Equally Censored Data. A cohort of n subjects are monitored in

the time interval [0,T], where T is ﬁxed in advance. Suppose that r failures

are observed (r can be any number from 0 to n) at times t

, t

,... , t

≤ T.

There are (n

−r) subjects that survived the entire period [0, T], and their

failure times are not observed.

Suppose that f (t) is the density of a lifetime. The likelihood is

i=1

f (t

)

(

1 −F(T)

)

n−r

for some normalizing constant C.

728 18 Inference for Censored Data and Survival Analysis

(a) Express the likelihood L for the exponential lifetime distribution, that

is, f (t)

=λe

−λt

, t ≥0, and F(t) =1 − e

−λt

, t ≥0.

(b) Take the log of the likelihood,

` =log(L).

` with respect to λ. Set the derivative to 0 and

solve for

λ. If the solution

λ maximizes `, (`

(

λ) <0), then

λ is the MLE of

λ.

(d) Show that the MLE is

mle

+(n −r)T

(e) If in the interval [0,8] four subjects failed at times t

=2, t

=5/2, t

=4,

and t

= 5 and two subjects survived without failure, ﬁnd the estimator

λ, assuming that the lifetime has an exponential distribution. Check your

calculations with

data = [2 2.5 4 5 8 8];

cens = [0 0 0 0 1 1];

lamrec = mle(data, ’distribution’,’exponential’,’censoring’,cens)

lammle = 1/lamrec %MATLAB uses reciprocal parametrization

(f) What is the MLE of λ if the unobserved failure times are ignored, that

is, only four observed failure times, t

, t

, and t

, are used.

18.7. Malignant Melanoma. Survival times of 256 males with malignant

melanoma who had metastases on admission to the MD Anderson Clinic

are reported (McDonald, 1963). The period of admission was between 1944

and 1960. The table shows survival times (no censoring) of the patients on

a year scale.

Find the mean survival time.

Survival time Number of patients Number of patients

(in years) admitted in time interval dying in time interval

0–1 256 167

1–2 89 48

2–3 41 23

3–4 18 6

4–5 12 3

5–6 9 6

6–7 3 1

7–8 2 1

8–9 1 1

9+ 0 0

18.8. Rayleigh Survival Times. It was observed that in clinical studies dealing

with cancer survival times follow Rayleigh distribution with pdf

f (x)

=2λte

−λt

, t ≥0, λ >0.

18.6 Exercises 729

(a) Show that the hazard function is linear.

(b) Find the mean survival time as a function of

λ.

,... , t

observed and t

k+1

,... , t

censored times, show that the like-

lihood is proportional to

∝λ

exp{−λ

i=1

If the prior on

λ is gamma G (α,β), show that the posterior is gamma G (α +

k, β +

(d) Show that the Bayes estimators of the hazard and survival functions

are

(t) =

2(k +α)t

β +

and

(t) =

1 +

β +

−(k+α)

18.9. Western White Clematis. Muenchow (1986) tested whether male or fe-

male ﬂowers (of Clematis ligusticifolia) were equally attractive to insects.

The data in table represent waiting times (in minutes), which includes cen-

sored data (observations with +).

Male ﬂowers Female ﬂowers

1 9 27 1 19 57

1 9 27

2 23 59

2 9 30

4 23 67

2 11 31

4 26 71

4 11 35

5 28 75

4 14 36

6 29 75+

5 14 40

7 29 78+

5 14 43

7 29 81

6 16 54

8 30 90+

6 16 61

8 32 94+

6 17 68

8 35 96

7 17 69

9 35 96+

7 18 70

14 37 100+

8 19 83

15 39 102+

8 19 95

18 43 105+

8 19 102+

18 56

104+

Compare survival functions for Male and Female ﬂowers.

730 18 Inference for Censored Data and Survival Analysis

MATLAB AND WINBUGS FILES AND DATA SETS USED IN THIS CHAPTER

http://springer.bme.gatech.edu/Ch18.Survival/

chada.m, cordband.m, ImmunoPerox.m, kmcdfsm.m, KMmuenchow.m,

kmplot.m, limphomaLR.m, logrank.m, logrankdata.mat, logrankmod.m,

Muenchow.m, PBC.m, simulation1.m, simulation2.m, smokingbugs.m,

weibullsim.m

ibrahim1.odc, ibrahim2.odc, Immunoperoxidase.odc, Leukemia.odc,

photocar.odc, Smoking.odc

gehan.dat KMmuenchow.txt limphoma.mat pbc.xls pbcdata.dat

prostatecanc.dat smokingoutbugs.mat

CHAPTER REFERENCES

Armitage, P. and Berry, G. (1994). Statistical Methods in Medical Research, 3rd edn. Black-

well, London.

Banerjee, S. and Carlin, B. P. (2004). Parametric spatial cure rate models for interval-

censored time-to-relapse data. Biometrics, 60, 268–275.

Bhattacharjee, A., Richards, W. G., Staunton, J., Li, C., Monti, S., Vasa, P., Ladd, C., Be-

heshti, J., Bueno, R., Gillette, M., Loda, M., Weber, G., Mark, E. J., Lander, E. S.,

Wong, W., Johnson, B. E., Golub, T. R., Sugarbaker, D. J., and Meyerson, M. (2001).

Classiﬁcation of human lung carcinomas by mRNA expression proﬁling reveals dis-

tinct adenocarcinoma subclasses. Proc. Natl. Acad. Sci. U.S.A., 98, 24, 13790–13795.

Cardillo, G. (2008). LogRank: comparing survival curves of two groups using the log rank

test,

http://www.mathworks.com/matlabcentral/fileexchange/22317.

Chadha, V., Warady, B. A., Blowey, D. L., Simckes, A. M., and Alon, U. S. (2000). Tenckhoff

catheters prove superior to Cook catheters in pediatric acute peritoneal dialysis. Am.

J. Kidney Dis., 35, 6, 1111–1116.

Cox, D. R. (1972). Regression models and life tables (with discussion). J. R. Stat. Soc. Ser. B,

34, 187–220.

Cox, D. R. and Oakes, D. (1984). Analysis of Survival Data. Chapman & Hall, New York.

Crowder, M. J., Kimber, A. C., Smith, R. L. and Sweeting, T. J. (1991). Analysis of Reliability

Data. Chapman and Hall, London

Dellaportas, P. and Smith, A. F. M. (1993). Bayesian inference for generalized linear and

proportional hazards models via Gibbs sampling. Appl. Stat., 42, 443–460.