Logistic regression with spatial autocorrelation structure
On Mon, Feb 7, 2011 at 1:26 PM, <Dale.W.Steele at gmail.com> wrote:
Suspect I'm not asking the question well. I'm clear on how to get a point estimate for the predicted probability in both cases. How much uncertainty is inherent in these prediction(s)? Not yet clear to me how to obtain something like a 95% confidence interval. Best. --Dale
That is considerably more complex and I'm not quite sure how to answer it in the general case. In fact, I don't really know how the standard error of a predicted probability would be calculated even for a fixed-effects GLM.
On Feb 7, 2011 1:56pm, Douglas Bates <bates at stat.wisc.edu> wrote:
On Mon, Feb 7, 2011 at 10:58 AM, ?Dale.W.Steele at gmail.com> wrote:
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?
Didn't I answer 2 in my previous message in the part about
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 It is certainly possible to automate the process a bit more with a predict method but, for the time being, I think that should be the way to go about it. If you have a simple random effect for doc, as in model m2, then to get a prediction for a specific physician you add the random effect for, say, doc == "18", to the value of eta before transforming it. You can get the conditional modes of the random effects with the ranef extractor.
Thanks. --Dale
On Feb 1, 2011 1:08pm, Douglas Bates bates at stat.wisc.edu> wrote:
On Tue, Feb 1, 2011 at 12:05 PM, Douglas Bates bates at stat.wisc.edu> wrote:
On Mon, Jan 31, 2011 at 11:45 AM, Dale Steele dale.w.steele at 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