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

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

17.3 Poisson Regression 681

0 2 4 6

−2

Pearson Res.

0 2 4 6

−2

Deviance Res.

0 2 4 6

−2

Friedman-Tukey Res.

0 2 4 6

−2

Anscombe Res.

Fig. 17.8 (a) Four types of residual in a Poisson regression ﬁt of IHGA data: Pearson, de-

viance, Friedman–Tukey, and Anscombe. The residuals are plotted against responses y.

%log-likelihood for saturated model

logliksat = sum(y.

log(y+eps)-y-log(factorial(y))) %-338.1663

m2LL = -2

sum( y .

log(yhat./(y + eps)) ) %819.8369

deviance = sum(rdev.^2) %819.8369

dev %819.8369 from glmfit

-2

(loglik - logliksat) % 819.8369

The following is a Bayesian model ﬁt in WinBUGS ( geriatric.odc).

682 17 Regression for Binary and Count Data

model

{

for (i in 1:n)

{

y[i] ~ dpois(lambda[i] )

log( lambda[i]) <- beta.0 + beta.1

x[i]

}

beta.0 ~ dnorm(0, 0.0001)

beta.1 ~ dnorm(0, 0.0001)

lambda.C <- exp(beta.0)

lambda.E <- exp(beta.0 + beta.1 )

diff <- lambda.E - lambda.C

meffect <- exp( beta.1 )

}

DATA

list( y = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

...

2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,

4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 7, 7, 0, 0, 0, 0, 0, 0, 0,

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

...

2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,

4, 4, 4, 5, 6 ),

x = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

...

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1,

1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

...

1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

1, 1, 1, 1, 1 ), n = 572 )

INITS

list( beta.0 = 0.0, beta.1 = 0.0 )

mean sd MC error val2.5pc median val97.5pc start sample

beta.0 -0.05915 0.06093 2.279E-4 -0.1801 -0.05861 0.05881 1001 100000

beta.1 -0.2072 0.09096 3.374E-4 -0.3861 -0.2071 -0.02923 1001 100000

deviance 1498.0 2.012 0.006664 1496.0 1498.0 1504.0 1001 100000

diff -0.1764 0.07737 2.872E-4 -0.3285 -0.1763 -0.02493 1001 100000

lambda.C 0.9443 0.05749 2.137E-4 0.8352 0.9431 1.061 1001 100000

lambda.T 0.7679 0.05188 1.784E-4 0.6693 0.7668 0.8732 1001 100000

meffect 0.8162 0.07437 2.76E-4 0.6797 0.8129 0.9712 1001 100000

17.3 Poisson Regression 683

Example 17.5. Cellular Differentiation Data. In a biomedical study of the

immunoactivating ability of the agents TNF (tumor necrosis factor) and IFN

(interferon) to induce cell differentiation, the number of cells that exhibited

markers of differentiation after exposure to TNF or IFN was recorded (Pier-

gorsch et al., 1988; Fahrmeir and Tutz, 1994). At each of the 16 dose combina-

tions of TNF/IFN, 200 cells were examined. The number y of differentiating

cells corresponding to a TNF/IFN combination are given in Table 17.6.

Table 17.6 Cellular differentiation data.

Number Dose of Dose of Number Dose of Dose of

cells diff TNF (U/ml) IFN (U/ml) cells diff TNF (U/ml) IFN (U/ml)

11 0 0 31 10 0

18 0 4 68 10 4

20 0 20 69 10 20

39 0 100 128 10 100

22 1 0 102 100 0

38 1 4 171 100 4

52 1 20 180 100 20

69 1 100 193 100 100

The suggested model is Poisson with the form

λ =E(y|TNF, IFN) =exp

+β

TNF +β

IFN +β

TNF×IFN

From

load ’celular.dat’

number = celular(:,1);

TNF = celular(:,2);

IFN = celular(:,3);

[b, dev, stats] = glmfit([TNF IFN TNF.

IFN],number,...

’poisson’,’link’,’log’)

the estimators for β

–β

are

3.43463 0.01553 0.00895 −0.0000567

Since b

<0.0001, it is tempting to drop the interaction term. However, since

the standard error of b

is s.e.(b

) =0.000013484, Wald’s Z = b

/s.e.(b

) statis-

tic is

−4.2050, suggesting that the term TNF×IFN might be signiﬁcant. How-

ever, the overdispersion parameter, which theoretically should be

stats.s=1,

is estimated as

stats.s=3.42566, and the Wald statistic should be adjusted

to Z

=−0.0000567/(3.42566 ×0.000013484) = 1.2275. Since the p-value is 2

normcdf(-1.2275) = 0.2196

, after all, the interaction term turns out not to be

signiﬁcant and the additive model could be ﬁt:

[b, dev, stats] = glmfit([TNF IFN],number,...

’poisson’,’link’,’log’)

684 17 Regression for Binary and Count Data

which gives estimates b

= 3.57311, b

= 0.01314, and b

= 0.00585. Details

can be found in

celular.m.

The additive model was also ﬁt in a Bayesian manner.

model

{

for (i in 1:n)

{

numbercells[i] ~ dpois(lambda[i])

lambda[i] <- exp(beta0 + beta1

tnf[i] + beta2

ifn[i])

}

beta0 ~ dnorm(0, 0.00001)

beta1 ~ dnorm(0, 0.00001)

beta2 ~ dnorm(0, 0.00001)

}

DATA

list(n=16,

numbercells = c(11,18,20,39,22,38,52,69,31,68,69,128,102,171,180,193),

tnf = c(0,0,0,0, 1,1,1,1, 10,10,10,10, 100,100,100,100),

ifn = c(0,4,20,100, 0,2,20,100, 0,4,20,100, 0,4,20,100 ) )

INITS

list(b0=0, b1=0, b2=0)

mean sd MC error val2.5pc median val97.5pc start sample

b0 3.573 0.05139 0.001267 3.473 3.574 3.672 1001 100000

b1 0.01313 5.921E-4 1.18E-5 0.01197 0.01314 0.01431 1001 100000

b2 0.00585 6.399E-4 1.142E-5 0.004585 0.005855 0.007086 1001 100000

Note that, because of the noninformative priors on β

, β

, and β

, the

Bayesian estimators almost coincide with the MLEs from

glmfit.

17.4 Log-linear Models

Poisson regression can be employed in the context of contingency tables

(Chap. 14). The logarithm of the table cell count is modeled in a linear fashion.

Let in an r

× c contingency table the probability of a cell (i, j) be p

i j

. If

n subjects are cross-tabulated, let n

i j

be the number of subjects classiﬁed in

the cell (i, j). The count n

i j

is realization of random variable N

i j

. We assume

that the sample size n is random, since in that case the cell frequency N

i j

a Poisson random variable with the intensity

i j

. If the sample size n is ﬁxed,

then the N

i j

s are realizations of a multinomial random variable. In Fisher’s

exact test context (p. 546), we saw that if in addition the marginal counts are

ﬁxed, the N

i j

s are hypergeometric random variables.

Then the expected table is

17.4 Log-linear Models 685

1 2 ··· c Total

1 np

1·

2 np

2·

r np

r·

Total np

·1

·2

·c

Note that both n

i j

and e

i j

= n

i·

×n

·j

/n involve observations and empirical

marginal probabilities deﬁned by the observed marginal frequencies. If we de-

note

i j

=EN

i j

= np

i j

, then both n

i j

and e

i j

are estimators of µ

i j

; the ﬁrst is

unconstrained and the second is constrained when the independence of factors

is assumed, p

i j

= p

i·

·j

. This point is important since e

i j

are “expected” under

independence once the table is observed and ﬁxed, while

i j

are expectations

of the random variables N

i j

The log-linear model for the expected frequency

i j

is given as

log

i j

=λ

+λ

+(λ

)

i j

, i =1, . .., r; j =1, . .., c,

where

and λ

are contributions by row and column, respectively, and λ

i j

is a row–column interaction term.

This model is called a saturated log-linear model and is similar to the two-

factor ANOVA model. Note that an unconstrained model has (1

+r +c +r ×c)

free parameters but only r

×c observations, n

i j

. To make this overparameter-

ized model identiﬁable, constraints on parameters are imposed. A standard

choice is STZ:

= 0,

(λ

)

i j

(λ

)

i j

= 0. Thus, with the

constraints we have 1

−(r −1) +(c +1) +(r −1)(c −1) = r × c free parameters,

which is equal to the number of observations, and the model gives a perfect

ﬁt for the observed frequencies. This is the reason why this model is called

saturated.

The hypothesis of independence in a contingency table that was discussed

Chap. 14 has simple form:

: λ

i j

=0, i =2,... , r; j =2,..., c.

Under H

the log-linear model becomes additive:

log

i j

=λ

+λ

, i =1, . .., r; j =1, . .., c.

The MLEs of components in the log-linear model [not derived here, but see

Agresti (2002)] are

686 17 Regression for Binary and Count Data

i, j

log n

i j

log n

i j

−

log n

i j

−

i j

= log n

i j

−(

If any n

i j

is equal to 0, then all entries in the table are replaced by n

i j

+0.5.

Traditional analysis involves testing that particular

λ components are

equal to 0. One approach, often implemented in statistical software, would

be to ﬁnd the variance of

λ (using a hypergeometric model), Var (

λ), and use

the statistic (

λ)

/Var (

λ) that has a χ

-distribution with one degree of freedom

(

λ here is any of λ

, λ

, or λ

i j

). Large values of this statistic are critical

for H

Next, we focus on the Bayesian analysis of a log-linear model.

Example 17.6. In Exercise 14.9 concerning a psychological experiment, seek-

ing assistance (help) was dependent on a subject’s perception of the number of

bystanders. The resulting chi-square statistic

= 7.908 was signiﬁcant with

a p-value of about 2%. We revisit this exercise and provide a Bayesian solution

using a loglinear approach. The WinBUGS program

bystanders.odc is used

to conduct statistical inference.

model{

for (i in 1:r) for (j in 1:c) n[i,j] ~ dpois(mu[i,j]);

log(mu[i,j]) <- lambda0+ lambdaR[i]+lambdaC[j]+lambdaRC[i,j]

}

lambda0 ~ dnorm(0,prec)

lambdaR[1] <- 0

lambdaC[1] <- 0

for (i in 2:r) { lambdaR[i] ~ dnorm(0,prec) }

for (i in 2:c) { lambdaC[i] ~ dnorm(0,prec) }

for (j in 1 : c) { lambdaRC[1, j] <- 0 }

for (i in 2 : r) { lambdaRC[i, 1] <- 0;

for (j in 2 : c) { lambdaRC[i, j] ~ dnorm(0, prec)}

}

DATA

list(r=3,c=2,prec=0.0001,

n=structure(.Data=c(11,2,16,10,4,9),.Dim=c(3,2)))

17.4 Log-linear Models 687

INITS

list(lambda0 = 0,lambdaR = c(NA,0,0), lambdaC=c(NA,0),

lambdaRC = structure(.Data = c(NA, NA, NA,0,NA,0),

.Dim = c(3,2)) )

mean sd MC error val2.5pc median val97.5pc start sample

lambda0 2.351 0.3075 0.003347 1.707 2.366 2.908 2001 100000

lambdaC[2] −1.922 0.8446 0.01376 −3.796 −1.85 −0.481 2001 100000

lambdaR[2] 0.3895 0.3977 0.004085 −0.372 0.3842 1.186 2001 100000

lambdaR[3] −1.089 0.6145 0.006468 −2.375 −1.061 0.03655 2001 100000

lambdaRC[2,2] 1.435 0.9372 0.01451 −0.2453 1.38 3.465 2001 100000

lambdaRC[3,2] 2.801 1.051 0.01653 0.9045 2.744 5.041 2001 100000

The hypothesis of independence is assessed by testing that all λ

are

equal to 0. In this case

lambdaRC[1,1], lambdaRC[1,2], and lambdaRC[2,1] are

set to 0 because of identiﬁability constraints. The 95% credible set for inter-

action

lambdaRC[2,2] contains 0 – therefore, this interaction is not signiﬁcant.

On the other hand,

lambdaRC[3,2] is signiﬁcantly positive since the credible

set [0.9045, 5.041] does not contain 0.

Example 17.7. Upton–Fingleton Square. An example from Upton and Fin-

gleton (1985) concerns ﬁnding directional trends in spatial count data. The

simple count data set is given as follows:

0 0 0 1 0

3 3 5 2 0

2 3 6 2 5

1 5 4 6 7

6 2 4 3 4

We are interested in testing if there are any north–south or east–west trends

present in the spatial pattern. The idea of Upton and Fingleton was to estab-

lish a Poisson regression where the response was the intensity n

i j

in the loca-

tion (i, j) and the covariates are x-coordinate, related to the east–west trend,

and y-coordinate, related to the north–south trend.

Thus, the observed frequencies n

i j

will be modeled as

E(n

i j

) =exp{β

+β

}, i, j =1,..., 5,

where x

= i and y

= j. Signiﬁcant β

or β

in a well-ﬁtting Poisson regression

will indicate the presence of corresponding trends (see

UptonFingleton.m

for details).

688 17 Regression for Binary and Count Data

17.5 Exercises

17.1. Blood Pressure and Heart Disease. This example is based on data

given by Cornﬁeld (1962). A sample of male residents of Framingham, Mas-

sachusetts, aged 40–59, were chosen. During a 6-year follow-up period they

were classiﬁed according to several factors, including blood pressure and

whether they had developed coronary heart disease. The results for these

two variables are given in the table below. The covariate for blood pressure

represents the interval, for example, 122 stands for blood pressure in the

interval 117–126.

Blood pressure No disease Disease Total

112 (<117) 156 3 159

122 (117–126) 252 17 269

132 (127–136) 284 12 296

142 (137–146) 271 16 287

152 (147–156) 139 12 151

162 (157–166) 85 8 93

177 (167–186) 99 16 115

197 (

> 186) 43 8 51

Using logistic regression, estimate the probability of disease for a person

with an average blood pressure equal to 158.

17.2. Blood Pressure and Heart Disease in WinBUGS. Use data from the

previous exercise. Parameter p represents the probability of development

of coronary disease and can be estimated as

1 +e

where b

and b

are Bayes estimators of β

and β

obtained by WinBUGS.

Use noninfromative priors on

and β

(a) What are values b

and b

? Provide plots of posterior densities for b

and b

(b) What are the 95% credible intervals for

and β

pressure of 158.

17.3. Sex of Turtles and Incubation Temperature. Below are data on the

relationship between the proportion of male turtles and incubation temper-

ature for turtle eggs from New Mexico.

Temp in

◦

C Male Female

27.2 0 10

28.3 8 4

29.9 8 2

17.5 Exercises 689

Predict the probability of female turtle for a temperature of 28

◦

17.4. Health Promotion. Students at the University of the Best in England

(UBE) investigated the use of a health promotion video in a doctor’s surgery.

Covariates Age and Amount of Weekly Exercise for a sample of 30 men

were obtained, and each man was asked a series of questions on the video.

On the basis of the responses to these questions the psychologist simply

recorded whether the promotion video was Effective or Not Effective.

The data collected are as follows:

Number Age Exercise Code of response Response

1 27 3 0 Not Effective

2 26 5 0 Not Effective

3 28 10 1 Effective

4 40 4 1 Effective

5 19 3 0 Not Effective

6 36 14 0 Not Effective

7 41 5 1 Effective

8 27 5 0 Not Effective

9 33 1 0 Not Effective

10 34 2 1 Effective

11 51 8 1 Effective

12 33 21 1 Effective

13 23 12 0 Not Effective

14 41 19 1 Effective

15 38 2 0 Not Effective

16 25 9 0 Not Effective

17 40 6 0 Not Effective

18 36 25 1 Effective

19 40 13 0 Not Effective

20 39 3 0 Not Effective

21 45 10 1 Effective

22 39 5 0 Not Effective

23 40 2 1 Effective

24 20 1 0 Not Effective

25 47 9 1 Effective

26 31 17 1 Effective

27 37 7 1 Effective

28 30 0 0 Not Effective

29 32 13 1 Effective

30 25 10 0 Not Effective

(a) Use the ﬁtted model to estimate the probability that a male of age 40

who exercises 10 hours a week would ﬁnd the video “effective.”

Comment: The linear predictor is

−7.498 +0.17465 Age +0.16324 Exercise.

690 17 Regression for Binary and Count Data

The positive coefﬁcient with covariate Age means that older subjects tend to

respond “video is effective” with higher probabilities. Similarly, the positive

coefﬁcient with predictor Exercise indicates that increasing the values of

Exercise also increases the probabilities for the response “video is effective.”

(b) Comment on the ﬁt of the model based on deviance and residuals.

17.5. PONV. Despite advances over the past decade, including the advent of 5-

receptor antagonists, combination therapy, and multimodal strategies,

postoperative nausea and vomiting (PONV) remains a serious and frequent

adverse event associated with surgery and anesthesia. PONV can be very

distressing for patients, can lead to medical complications, and impose eco-

nomic burdens. A meta-analysis of several studies gives rates of 37% for

nausea and 20% for vomiting in patients undergoing general anesthesia.

However, indiscriminate prophylaxis is not recommended (the “prevent-or-

cure” dilemma).

There are considerable variations in the reported incidence of PONV, which

can be attributed to a number of factors. Risk factors for PONV can be di-

vided into patient risk factors, procedural risk factors, anesthetic risk fac-

tors, and postoperative risk factors. The main and well-understood risk fac-

tors are gender, history of motion sickness/PONV, smoking status, and use

of postoperative opioids.

A data set

PONV.xls or PONV.mat (courtesy of Jelena Velickovic, MD

anesthesiologist from Belgrade) contains records of 916 patients consisting

of some demographic, anamnetic, clinical, and procedural variables. Sev-

eral variables are of interest to be modeled and controlled: manifestation

of PONV from 0 to 2 hours after surgery, PONV from 2 to 24 hours, and

PONV from 0 to 24 hours (

PONV0to2, PONV2to24, and PONV0to24).

Three score variables (

SinclairScore, ApfelScore, and LelaScore)

summarize the relevant demographic and clinical information prior to

surgery with the goal of predicting PONV.

The starter ﬁle

ponv.m gives a basic explanation of variables and helps to

read all the variables into MATLAB.

Fit a logistic model for predicting

PONV0to24 based on a modiﬁed Sin-

clair Score deﬁned as

MSS = SinclairScore + 1/20

LelaScore. What is

the probability of

PONV0to24 for a person with MSS = 1.3?

17.6. Mannose-6-phosphate Isomerase. McDonald (1985) counted allele fre-

quencies at the mannose-6-phosphate isomerase (Mpi) locus in the amphi-

pod crustacean Megalorchestia californiana, which lives on sandy beaches

of the Paciﬁc coast of North America. There were two common alleles,

Mpi90 and Mpi100. The latitude of each collection location, the number

of each allele, and the proportion of the Mpi100 allele are shown here: