Dear Prof. Bates - Thanks for your response. How would I calculate a 95%  
confidence interval for
the probability of appendicitis in a new patient given a particular set of  
covariates for

1) patients enrolled by a specific clinician (eg. doc == "18")?
2) some future clinician?

Thanks. --Dale

On Feb 1, 2011 1:08pm, Douglas Bates <bates@stat.wisc.edu> wrote:
> On Tue, Feb 1, 2011 at 12:05 PM, Douglas Bates bates@stat.wisc.edu> wrote:

> > On Mon, Jan 31, 2011 at 11:45 AM, Dale Steele dale.w.steele@gmail.com>  
> wrote:

> >> Dear mixed-modeling experts,

> >>

> >> I'm interested in modeling the probability of appendicitis in patients

> >> with abdominal pain.

> >>

> >> The R binary data file 'http://www.ped-em.org/appy.rda' contains

> >> the following variables from a pilot study of 138 children with

> >> abdominal pain.

> >

> > Thank you for providing the data.

> >

> >> 'dx' eventual diagnosis: 0=no appendicitis, 1=appendicitis

> >> 'gender' Male/Female

> >> 'wbc' total white blood cell count

> >> 'priorprob' Clinical predicted probability of appendicitis

> >> 'doc' doctor who assigned 'priorprob'

> >>

> >> After taking a history and performing a physical examination, the ER  
> doctor

> >> was asked to make a vertical mark on a 100 mm horizontal line to  
> represent

> >> her estimate of the (percent) probability that the patient had  
> appendicitis.

> >>

> >> My initial thought was to fit a multiple logistic regression model:

> >>

> >> m1
> >>

> >> However, it seems likely that each doctor interpreted the probability  
> scale

> >> differently. The 23 doctors evaluated from 1 to 17 patients each. I'm

> >> not primarily interest in predictions by a specific clinician. Thus,

> >> it seems to make sense to fit a generalized linear mixed model.

> >>

> >> At this point I get muddled. Have I correctly specified a random

> >> intercept model (m2) and a random intercept/random slope model (m3)?

> >> Are there other sensible models?

> >

> > Yes, you have correctly specified them. However, if you check the

> > output from the fits

> >> (m2
> > Generalized linear mixed model fit by the Laplace approximation

> > Formula: dx ~ priorprob + gender + wbc + (1 | doc)

> > Data: appy

> > AIC BIC logLik deviance

> > 108.1 122.7 -49.03 98.07

> > Random effects:

> > Groups Name Variance Std.Dev.

> > doc (Intercept) 0.82922 0.91062

> > Number of obs: 138, groups: doc, 23

> >

> > Fixed effects:

> > Estimate Std. Error z value Pr(>|z|)

> > (Intercept) -8.28862 1.45666 -5.690 1.27e-08

> > priorprob 0.05025 0.01291 3.892 9.96e-05

> > genderMale 1.18705 0.56752 2.092 0.0365

> > wbc 0.34130 0.07090 4.814 1.48e-06

> >

> > Correlation of Fixed Effects:

> > (Intr) prrprb gndrMl

> > priorprob -0.775

> > genderMale -0.293 0.130

> > wbc -0.790 0.353 0.040

> >> (m3
> > Generalized linear mixed model fit by the Laplace approximation

> > Formula: dx ~ priorprob + gender + wbc + (priorprob | doc)

> > Data: appy

> > AIC BIC logLik deviance

> > 111.5 132 -48.76 97.52

> > Random effects:

> > Groups Name Variance Std.Dev. Corr

> > doc (Intercept) 0.02897010 0.170206

> > priorprob 0.00016799 0.012961 1.000

> > Number of obs: 138, groups: doc, 23

> >

> > Fixed effects:

> > Estimate Std. Error z value Pr(>|z|)

> > (Intercept) -8.23977 1.44408 -5.706 1.16e-08

> > priorprob 0.04957 0.01328 3.733 0.000189

> > genderMale 1.16593 0.57481 2.028 0.042521

> > wbc 0.34485 0.07056 4.887 1.02e-06

> >

> > Correlation of Fixed Effects:

> > (Intr) prrprb gndrMl

> > priorprob -0.757

> > genderMale -0.299 0.141

> > wbc -0.800 0.350 0.030

> >> anova(m2, m3)

> > Data: appy

> > Models:

> > m2: dx ~ priorprob + gender + wbc + (1 | doc)

> > m3: dx ~ priorprob + gender + wbc + (priorprob | doc)

> > Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)

> > m2 5 108.07 122.70 -49.034

> > m3 7 111.52 132.01 -48.759 0.5501 2 0.7595

> >> (m4
> > Generalized linear mixed model fit by the Laplace approximation

> > Formula: dx ~ priorprob + wbc + (1 | doc)

> > Data: appy

> > AIC BIC logLik deviance

> > 110.5 122.2 -51.24 102.5

> > Random effects:

> > Groups Name Variance Std.Dev.

> > doc (Intercept) 0.63458 0.7966

> > Number of obs: 138, groups: doc, 23

> >

> > Fixed effects:

> > Estimate Std. Error z value Pr(>|z|)

> > (Intercept) -7.62797 1.33760 -5.703 1.18e-08

> > priorprob 0.04834 0.01229 3.934 8.35e-05

> > wbc 0.34395 0.06771 5.080 3.77e-07

> >

> > Correlation of Fixed Effects:

> > (Intr) prrprb

> > priorprob -0.786

> > wbc -0.813 0.357

> >> anova(m4, m3)

> > Data: appy

> > Models:

> > m4: dx ~ priorprob + wbc + (1 | doc)

> > m3: dx ~ priorprob + gender + wbc + (priorprob | doc)

> > Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)

> > m4 4 110.48 122.19 -51.240

> > m3 7 111.52 132.01 -48.759 4.9639 3 0.1745



> Sorry, that test should have been



> > anova(m4, m2)

> Data: appy

> Models:

> m4: dx ~ priorprob + wbc + (1 | doc)

> m2: dx ~ priorprob + gender + wbc + (1 | doc)

> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)

> m4 4 110.48 122.19 -51.240

> m2 5 108.07 122.70 -49.034 4.4138 1 0.03565



> > you will see that the random intercept/random slope model produces a

> > degenerate fit (estimated correlation of the within-doctor random

> > effects is -1.000), which is not significantly better than the random

> > intercept fit. I would therefore use m2. (Because the fixed-effect

> > for gender had a z-value close to 2 I performed the more reliable

> > likelihood-ratio test, just to check.)

> >

> >

> >> library(lme4)

> >> m2
> >> family=binomial, data=appy)

> >>

> >> m3
> >> family=binomial, data=appy)

> >>

> >> My ultimate goal is to estimate the probability of appendicitis

> >> (and a prediction interval), given a specific 'gender', 'wbc' and

> >> 'priorprob' assigned by a doctor with similar diagnostic ability to

> >> those who participated in our pilot study. I'm stuck on how to code  
> this

> >> prediction.

> >

> > The prediction will be based on the fixed-effects only. Because it

> > applies to a doctor not in this study, we assign the random effect for

> > that doctor to be zero, which is the expected value in the absence of

> > any information on that doctor.

> >

> > The fixef extractor returns the estimates of the fixed-effects

> > parameters. For a female patient with median prior probability (55%)

> > and median white blood cell count (9.45) the estimated linear

> > predictor (eta = -2.299) corresponds to a probability of 9.1%

> >

> >

> >> fixef(m2)

> > (Intercept) priorprob genderMale wbc

> > -8.2886183 0.0502540 1.1870539 0.3412998

> >> eta
> >> eta

> > [,1]

> > [1,] -2.299365

> >> binomial()$linkinv(eta)

> > [,1]

> > [1,] 0.09117557

> >


	[[alternative HTML version deleted]]


