Benchmark dose estimation incorporating multiple data sources


Wheeler, M.W.; Bailer, A.John.

Risk Analysis 29(2): 249-256

2008


With the increased availability of toxicological hazard information arising from multiple experimental sources, risk assessors are often confronted with the challenge of synthesizing all available scientific information into an analysis. This analysis is further complicated because significant between-source heterogeneity/lab-to-lab variability is often evident. We estimate benchmark doses using hierarchical models to account for the observed heterogeneity. These models are used to construct source-specific and population-average estimates of the benchmark dose (BMD). This is illustrated with an analysis of the U.S. EPA Region IX's reference toxicity database on the effects of sodium chloride on reproduction in Ceriodaphnia dubia. Results show that such models may effectively account for the lab-source heterogeneity while producing BMD estimates that more truly reflect the variability of the system under study. Failing to account for such heterogeneity may result in estimates having confidence intervals that are overly narrow.

Risk
Analysis,
VoL
29,
No.
2,
2009
DOI:
10.1111/j.1539-6924.2008.01144.x
Benchmark
Dose
Estimation
Incorporating
Multiple
Data
Sources
Matthew
W.
Wheeler
s
and
A.
John
Bailer
2*
With
the
increased
availability
of
toxicological
hazard
information
arising
from
multiple
ex-
perimental
sources,
risk
assessors
are
often
confronted
with
the
challenge
of
synthesizing
all
available
scientific
information
into
an
analysis.
This
analysis
is
further
complicated
because
significant
between-source
heterogeneity/lab-to-lab
variability
is
often
evident.
We
estimate
benchmark
doses
using
hierarchical
models
to
account
for
the
observed
heterogeneity.
These
models
are
used
to
construct
source-specific
and
population-average
estimates
of
the
bench-
mark
dose
(BMD).
This
is
illustrated
with
an
analysis
of
the
U.S.
EPA
Region
IX's
reference
toxicity
database
on
the
effects
of
sodium chloride
on
reproduction
in
Ceriodaphnia
dubia.
Results
show
that
such
models
may
effectively
account
for
the
lab-source
heterogeneity
while
producing
BMD
estimates
that
more
truly
reflect
the
variability
of
the
system
under
study.
Failing
to
account
for
such
heterogeneity
may
result
in
estimates
having
confidence
intervals
that
are
overly
narrow.
KEY
WORDS:
Aquatic
toxicology;
Bayesian
methods;
generalized
linear
mixed
models;
hierarchical
models;
lab-to-lab
variability;
Poisson
responses
1.
INTRODUCTION
With
environmental
hazard
information
com-
piled
from
different
sources,
risk
analysts
are
of-
ten
confronted
with
the
issue
of
combining
data
from
these
multiple
sources
to
better
inform
the
decision-making
process.
Although
the
benefits
of
such
a
meta-analysis
are
obvious,
methodological
considerations
frequently
prevent
such
a
synthesis
since
data
from
multiple
sources
can
contain
source-
specific
variability.
In
this
case,
a
simple
pooled
anal-
ysis
is
ill
advised.
Common
risk
assessment
tools,
often
based
upon
generalized
linear
models,
have
I
National
Institute
for
Occupational
Safety
and
Health,
Risk
Evaluation
Branch,
Cincinnati,
OH,
USA.
2
Department
of
Mathematics
and
Statistics,
Miami
University,
Oxford,
OH,
USA.
*
Address
correspondence
to
A.
John
Bailer,
Dept.
of
Math.
&
Stat.,
Miami
University,
Oxford,
OH
45056,
USA;
tel:
513-529-3538;
fax:
513-529-1493;
.
no
built-in
structure
to
account
for
this
data-source
heterogeneity.
To
illustrate
this
problem,
consider
a
recent
eval-
uation
of
the
U.S.
EPA
Region
IX's
reference
tox-
icity
database,
which
included
data
from
multiple
labs
conducted
using
the
same
organisms
and
the
same
toxicants.
(
'
)
One
set
of
experiments
considered
the
impact
of
sodium
chloride
exposures
on
repro-
duction
of
the
cladoceran,
Ceriodaphnia
dubia.
In
this
analysis,
lab-to-lab
heterogeneity
was
evident
as
analyses
based
upon
individual
labs
yielded
consid-
erably
different
concentration-response
curves.
Fur-
ther,
lab-specific
benchmark
doses
(BMDs)
varied
significantly.
As
these
analyses
show
a
large
amount
of
heterogeneity
between
labs,
the
pooling
of
each
lab
source
into
a
single
dataset
may
result
in
an
anal-
ysis
that
mischaracterizes
the
central
estimate
of
the
concentration-response
curve.
Alternative
modeling
frameworks,
such
as
hierarchical
models
or
gener-
alized
linear
mixed-effects
models
(GLMM)
(2)
that
include
random
effects
for
each
lab,
may
provide
a
249
0272-4332/09/0100-0249$22.00/1®
2008
Society
for
Risk
Analysis
250
Wheeler
and
Bailer
systematic
structure
that
can
readily
account
for
the
added
lab/source
variability.
A
significant
focus
of
previous
work
incorporat-
ing
random
effects
in
risk
assessment
has
focused
on
developmental
toxicity
data
and
the
heterogeneity
often
encountered
between
individual
dams
in
rela-
tion
to
their
litter
size
as
well
as
fetal
weight.
Cat-
alono
and
Ryan
(3)
and
Dunson
et
al.
(4)
among
others
have
proposed
latent
class
models
that
take
into
ac-
count
this
correlative
source
structure
into
risk
anal-
ysis.
These
studies,
however,
have
focused
primarily
on
the
modeling
aspect
of
the
analysis,
with
little
em-
phasis
on
the
estimation
of
specific
endpoints
such
as
the
BMD.
Other
work
estimating
the
BMD
using
mixed
models
has
been
described
in
the
environmental
lit-
erature.
Piegorsch
and
Cox
(5)
reviewed
generalized
linear
models
with
random
effects
in
the
estimation
of
the
effective
dose
in
ordinal
and
continuous
re-
gression.
Simmons
et
al.
(6)
looked
at
modeling
mu-
tagenic
potencies
from
different
Salmonella
strains,
and
Coull
et
al.
(7
conducted
a
meta-analysis
for
methyl-mercury
using
a
Bayesian
hierarchical
model
to
combine
BMD
point
estimates
from
multiple
stud-
ies.
While
these
studies
have
focused
on
risk
assess-
ment,
most
have
looked
at
continuous
response
re-
gression,
with
less
emphasis
on
nonnormal
data.
We
consider
an
extension
of
this
work
and
suggest
a
general
hierarchical
framework
that
can
be
used
by
risk
assessors
to
estimate
the
dose/
concentration-response
curve,
as
well
as
the
BMD,
for
data
assumed
to
be
generated
from
the
exponen-
tial
family
of
distributions.
We
then
use
this
frame-
work
to
model
the
C.
dubia
concentration-response
data,
and
then
estimate
the
concentration
associated
with
specified
decrement
in
reproductive
output,
here
number
of
young
produced
in
three
broods.
2.
HIERARCHICAL
MODELING
FRAMEWORK
As
in
some
of
the
previously
cited
work,
we
use
a
hierarchical
model
to
estimate
adverse
effects
of
multiple-source
toxicology
data.
Here,
we
assume
that
the
mean
response
µ
c
,
at
toxicant
concentra-
tion
c,
can
be
represented
as
a
linear
function
of
the
concentration
through
the
link
function
g(p.
c
),
where
g(p.
c
)
is
a
linear
function
of
the
ambient
concentra-
tion,
and
the
concentration
effect
is
heterogeneous
across
sources
(note
that
an
internal
dose
can
be
sub-
stituted
for
the
concentration).
Given
these
assump-
tions,
a
multilevel
hierarchical
model
can
be
speci-
fled
allowing
for
the
heterogeneity
in
response
across
multiple
sources.
For
example,
the
multiple-source
toxicology
data
described
above
can
be
modeled
us-
ing
a
two-stage
hierarchical
model.
Here,
the
first
level
explains
the
population
mean
response,
given
the
jth
lab,
and
is
modeled
as:
g(µ
=
00j
+
B
li
c
+
0
2
jc
2
+
+
&
J
en.
(1)
The
second
level
models
each
lab's
heterogene-
ity
in
response
as:
Okj
=
Pk
bkj,
(2)
where
k
=
0,
1,
,
m,
13
=[13
0
,
0
1
.
.
n]
represents
the
fixed
effects,
and
bj
=
Poi,
b
1
,
...b
m
i]
^
Nm-Ft
(0,
D).
That
is,
b
;
is
a
random
vector
for
the
jth
lab,
which
is
distributed
multivariate
normal
with
mean
vector
0
and
variance—covariance
matrix
D
=
{a
rc
:
r
=
0,1
,
,m;
c
=
0,1
,
,m1.
Note
that
while
two
levels
are
used
in
this
formulation
to
represent
the
lab-source
variability,
additional
levels
could
be
in-
troduced
to
model
other
source
variability
that
may
be
present
in
the
analysis.
The
hierarchical
model
can
be
rewritten
into
a
single
linear
equation
describing
the
mean
by
substituting
Equation
(2)
into
Equation
(1),
i.e.:
g(pjc)
=
(
Po
+
NJ)
+
(01
+
bii)c
+
+
(P
m
+
b
m
i)cm.
Often
m
<
2
is
sufficiently
complex
to
model
the
aquatic
toxicology
data
we
consider.
Under
this
for-
mulation,
the
adverse
response
is
assumed
to
have
fixed
effects
associated
with
toxicant
exposure
with
additional
source
heterogeneity
accounted
for
by
source-independent
realizations
of
the
specified
ran-
dom
effect.
Given
an
estimate
of
p.
c
,
the
BMD
can
be
esti-
mated.
This
endpoint
is
dependent
upon
the
bench-
mark
response
(BMR)
value
p
as
well
as
the
un-
derlying
data-generating
mechanism.
Crump
(8)
first
discussed
this
quantity
in
terms
of
dichotomous
dose-
response
data,
Gaylor
and
Slikker,
(9)
among
others,
presented
a
formulation
in
terms
of
continuous
dose
response,
and
Bailer
and
Oris
(1
°
)
presented
a
sim-
ilar
formulation
for
concentration-response
studies
where
the
underlying
response
was
based
upon
count
data.
Regardless
of
the
response
distribution
used
to
model
the
response,
the
BMD(p),
the
benchmark
dose
specified
at
a
specified
BMR
=
p,
is
often
of
in-
terest
in
risk
assessment.
By
definition,
the
estimation
of
the
BMD
re-
quires
all
coefficients
of
the
above
model
to
be
Benchmark
Dose
Estimation
Incorporating
Multiple
Data
Sources
251
specified,
which
implies
that
the
estimate
is
source
specific
(i.e.,
it
relies
on
a
specific
realization
of
the
random
variables
b
1
).
This
source-specific
depen-
dence
prohibits
a
single
estimate,
which
includes
the
source
heterogeneity;
consequently,
a
population-
average
estimate
(i.e.,
the
expectation
of
the
BMD
across
all
realizations
of
the
random
effects)
may
be
more
useful
when
making
scientific
reasoning
from
combined
sources.
From
a
Bayesian
perspective
in
which
the
regres-
sion
coefficients
(,8
k
)
and
the
variance
components
(a„)
are
also
random
variables,
the
population-
averaged
estimate
is
an
average
BMD
over
all
ran-
dom
variables
impacting
this
endpoint,
i.e.:
BMDpa
=
f
BMD(p113,
b,
cr)n
-
($,
b,
cr)dfi
db
da,
(3)
where
740,
b,
a)
is
the
joint
posterior
distribution
of
these
quantities.
Thus,
the
population-average
value
integrates
out
the
regression
coefficients
(Os),
the
lab-specific
random
effects
(bs)
and
the
vari-
ance
components
associated
with
the
random
effects
(a
rc)
Multiple
methods
for
estimating
the
BMD
pa
could
be
considered.
While
it
is
possible
to
es-
timate
the
above
model
parameters
using
maxi-
mum
likelihood
methods,
computing
the
BMD
pa
and
subsequent
lower
bound
calculations
may
be
diffi-
cult
as
the
above
formulation
does
not
easily
yield
to
computationally
simple
methods.
A
bootstrap-
based
analysis
(11)
could
be
used
to
determine
this
lower
bound;
since
each
fit
of
a
GLMM
may
take
several
minutes
on
larger
problems,
an
adequate
number
of
bootstrap
replications
may
be
computa-
tionally
intractable.
A
Bayesian
formulation
of
the
problem
may
be
a
more
natural
choice.
In
this
setting
the
BMD
pa
is
a
quantity
that
can
be
readily
sampled
through
the
use
of
the
Markov
chain
Monte
Carlo
(MCMC)
algorithms.
(12)
Consequently,
we
formulate
the
rest
of
the
discussion
in
the
Bayesian
setting.
3.
CERIODAPHNIA
BMD
MODELING
We
return
to
the
C.
dubia
water
toxicity
tests
conducted
by
the
U.S.
EPA
Region
IX
between
July
18,
1989
and
June
16,
1993,
and
analyze
this
dataset
using
a
hierarchical
model/GLMM
in
a
Bayesian
framework.
The
data
were
obtained
from
a
series
of
experiments
conducted
in
seven
labs
with
multi-
ple
replicate
experiments
per
lab,
ranging
from
8
to
46
experiments
in
a
given
lab.
For
a
given
exper-
iment,
different
ambient
concentrations
of
the
ref-
erence
toxicant,
sodium
chloride,
was
introduced,
and
the
number
of
young
from
three
broods
was
recorded
as
the
measured
response.
The
testing
procedure
used
in
all
labs
is
fully
described
by
the
U.S.
EPA.
(13)
Bailer
et
al.
(1)
studied
this
database
and
concluded
that
experiments
exhibited
significant
lab-
to-lab
heterogeneity.
Their
analysis
suggests
that
any
model
based
upon
the
combined
data
from
the
seven
labs
should
account
for
this
interlaboratory
variabil-
ity;
consequently,
we
use
the
hierarchical
model
im-
plied
by
Equations
(1)
and
(2)
to
account
for
this
lab-
source
heterogeneity.
In
this
analysis,
let
Y
ff
i
c
represent
the
observed
number
of
young
given
the
jth
lab,
nth
experiment,
at
the
concentration
c.
Following
Bailer
et
a1.,
(1)
we
assume
that
that
Y
.
*
Poisson(p,
k
)
having
a
loga-
rithmic
link
function
that
relates
the
mean
response
to
a
quadratic
polynomial
of
the
fixed
and
random
effects,
i.e.:
lo
g
oi
l
o=
00+1)0+01+bl*
+02
+b2j)c2,
(4)
where
b
k.
'
is
modeled
as
in
Equation
(2).
Model
(4)
extends
the
model
of
Bailer
et
al.
(1)
by
adding
a
lab
random
effect
to
the
background,
linear,
and
quadratic
terms
of
the
equation.
The
model
can
be
thought
of
as
having
some
central
tendency
repre-
sented
by
the
C
Bs
with
each
lab
having
a
slightly
dif-
ferent
background,
as
well
as
an
unaccounted
for
concentration
effect
in
the
concentration-response
curve,
represented
by
the
bs.
Given
the
concentration-response
curve
(Equa-
tion
(4))
or
the
equivalent
formulation
of
Equations
(1)
and
(2),
the
lab-specific
BMD(p)
is
defined
as
the
concentration
that
decreases
the
expected
brood
size,
below
background
levels
(i.e.,
concentration
=
0),
by
some
proportion
p
for
that
lab.
This
value
is
defined
as
the
concentration
c
that
satisfies
the
equation
/L
k
=
(1
p)p,
j
°.
This
concentration
is
found
by
taking
the
largest
positive
root
of
the
equation
(pi
+
b
J
oc
+
02
+
bp)C
2
=
ln(1
p).
The
RIp
estimator
defined
by
Bailer
and
Oris
(1
°
)
is
a
BMD
estimator
defined
for
a
BMR
corresponding
to
inhibition
in
reproduc-
tion
relative
to
control-group
mean
response.
This
is
a
lab-specific
estimator
and
does
not
take
into
ac-
count
the
observed
heterogeneity
between
labs.
The
population-average
(Equation
(3))
estimator,
which
is
the
average
of
all
lab-source
estimates,
may
be
pre-
ferred
for
risk
management
purposes.
252
Wheeler
and
Bailer
Table
I.
Posterior
Mean
for
the
Population-Average
Parameters
of
Equation
(5)
(i.e.,
fiks),
as
Well
as
the
Mean
and
Median
Estimates
of
the
Lab
Effect
Variance
Fixed
Reg.
Effect
Mean
95%
Credible
Interval
Variance
Component
Mean
Median
2.5%
97.5%
Po
Th
P2
3.06
0.24
—0.81
2.89
—0.31
—1.09
3.23
0.82
—0.55
2
ao
=
Cr
00
2
=
ail
0
2
=
0
-22
0.05
0.51
0.12
0.04
0.38
0.09
Note:
A
large
portion
of
the
variability,
in
the
concentrations
of
interest,
is
accounted
for
in
the
linear
random
effect
term.
4.
HIERARCHICAL
MODELING
OF
THE
SODIUM
CHLORIDE
TOXICITY
DATABASE
Given
Equation
(1)
with
m
=
2
and
the
C.
du-
bia
experimental
data
recorded
from
the
seven
EPA
Region
IX
testing
facilities,
the
posterior
distribution
of
the
model
parameters
was
obtained
using
MCMC
algorithms
implemented
in
WinBugs
version
1.4.(
14
)
In
this
implementation,
all
parameters
in
Equations
(1)
and
(2)
were
assumed
to
have
relatively
flat
prior
information
with
the
following
priors
being
used:
O
f
N3(0,
D),
(
5
)
where
aoo
D=
[
0
0
or, more
simply,
Okj
"'
N(,8k,
a
kk)
for
k
=
0,
1,
2
with
uncertainty
in
these
parameters
reflected
by
,6k
—N(0,
10
6
),
and
cr
kk
^
InvGamma
(0.001,
0.001).
A
full
listing
of
the
WinBugs
code
used
to
es-
timate
the
model
is
described
in
the
Appendix.
In
all,
10,000
burn-in
samples
were
used
on
three
sep-
arate
chains,
and
convergence
was
monitored.
Once
convergence
was
achieved,
20,000
additional
samples
were
taken
from
the
posterior
distribution
for
the
9,
0,
and
a
(i.e.,
entries
from
D)
vectors.
Posterior
mean
estimates
for
the
0
and
a
pa-
rameters,
as
well
as
a
95%
credible
interval
for
the
0
vector,
are
presented
in
Table
I.
Here,
the
population-averaged
concentration-response
curve,
the
curve
defined
by
the
posterior
mean
values
of
(no,
01,
02),
suggests
an
increased
response
asso-
ciated
with
sodium
chloride
at
low
concentrations
(i.e.,
the
expectation
of
the
linear
term
is
positive
but
the
95%
credible
intervals
contain
0),
where
at
higher
concentrations
the
deleterious
effect
of
sodium
chloride
becomes
evident
(i.e.,
the
expecta-
tion
of
the
quadratic
term
is
negative).
Further,
the
linear
lab-source
effect
shows
the
highest
variabil-
ity.
This
reflects
the
observation
that
any
amount
of
sodium
chloride
may
lead
to
a
reduction
in
brood
size
in
some
labs
while
an
increased
brood
size
is
observed
at
low
NaCl
concentrations
in
other
labs.
The
concentration-response
patterns
for
the
differ-
ent
labs
are
displayed
in
Fig.
1.
Here
the
solid
curve
represents
the
population-averaged
concentration-
response
curve
and
the
dashed
lines
represent
the
ex-
pected
lab-specific
responses.
Table
II
presents
the
BMD/RIp
calculations
for
p
=
0.25
and
0.5,
for
the
lab-specific
and
population-
average
estimates.
(Note
that
in
the
original
papers
describing
these
data
and
potency
endpoints,
the
rel-
ative
inhibition
estimator/RIp
was
defined
and
esti-
mated.
This
quantity
corresponds
to
a
BMD
for
a
particular
definition
of
the
BMR.
Thus,
we
refer
to
the
RIp
and
BMD
interchangeably.)
From
this ta-
ble,
it
is
clear
that
there
is
a
significant
lab-to-lab
variability
in
these
estimates,
with
lab-specific
RI25
posterior
means
differing
by
a
factor
of
5
and
RI50
Fig.
I.
Posterior
means
of
the
expected
number
of
offspring
given
the
lab-source
variability.
The
solid
line
represents
the
population-
average
expected
response,
whereas
the
dotted
lines
represent
the
expected
response,
given
an
individual
lab.
0 0
0
a22
0
1
2
3
4
Concentration
Benchmark
Dose
Estimation
Incorporating
Multiple
Data
Sources
253
Table
II.
Posterior
Mean
and
Standard
Deviation
Estimates
for
the
Lab-
and
Population-Average
RIp
Distribution
EPA
Lab
RI25
RI50
Mean
SD
95%
Credible
Interval
Mean
SD
95%
Credible
Interval
2.5%
97.5%
2.5%
97.5%
Lab-specific
estimates
CAAQS
0.58
0.01
0.56
0.61
0.97
0.01
0.95
0.99
CAMEC
0.91
0.02
0.87
0.95
1.13
0.02
1.10
1.16
CAMEC1
0.30
0.02
0.27
0.33
0.60
0.02
0.56
0.64
CAOEE
0.67
0.02
0.64
0.70
1.01
0.01
0.98
1.03
CASRH
0.84
0.02
0.81
0.88
1.20
0.01
1.17
1.22
CAUCD
1.45
0.02
1.41
1.49
1.75
0.02
1.72
1.79
MNEPAD
0.94
0.02
0.91
0.98
1.21
0.01
1.18
1.24
Lab-average
estimates
0.80
0.37
0.43
1.36
1.12
0.38
0.74
1.67
Note:
Here,
for
the
RI25
and
RI50
distribution,
p
represents
a
25%
or
50%
decrement
in
expected
response
in
relation
to
control
group
responses.
Table
HI.
Comparison
of
the
Posterior
Estimates
of
the
flks
and
the
RIp
Estimates
for
the
Pooled-Data/Fixed
and
Hierarchical/Mixed
Effects
Generalized
Linear
Model
Analyses
Pooled-Data/Fixed
Effects
Analysis
Hierarchical/Mixed
Model
Analysis
95%
Credible
Interval
95%
Credible
Interval
Mean
2.5%
97.5%
Mean
2.5%
97.5%
Po
3.12
3.11
3.13
3.06
2.89
3.23
Pi
0.16
0.14
0.19
0.24
-0.31
0.82
P2
-0.68
-0.70
-0.66
-0.81
-1.09
-0.55
RIp25
0.78
0.77
0.80 0.80
0.43
1.36
RIp50
1.14
1.13
1.15
1.12
0.74
1.67
DIC
96800
87766
Note:
The
DIC
estimate
is
shown
as
a
method
of
Bayesian
model
comparison;
here
smaller
is
better.
posterior
means
differing
by
a
factor
of
3.
Individual
credible
intervals
are
tight
around
the
posterior
mean
in
each
lab,
suggesting
a
RIp
posterior
distribution
that
is
quite
narrow
for
a
given
lab;
however,
these
RIp
distributions
were
quite
different
among
labs.
This
suggests
that
inference
from
any
one
lab
may
not
be
generalizable
to
the
population
as
a
whole,
and
may
be
unsuitable
for
risk
management
pur-
poses.
The
table
further
shows
that
the
population-
averaged
RIp/BMD
pa
estimates,
based
upon
Equa-
tion
(2),
tend
to
have
much
larger
credible
intervals
as
they
explicitly
incorporate
the
lab-to-lab
hetero-
geneity.
One
can
look
at
the
effect
of
modeling
the
lab
heterogeneity
on
the
posterior
parameter
estimates,
by
comparing
Equation
(5)
to
a
model
that
does
not
include
random
effects,
i.e.,
essentially
pool
all
lab-
specific
data
into
a
single
dataset
and
then
fit
a
model
with
only
fixed
effects.
Here,
the
posterior
distribu-
tion
of
the
parameters
for
model
log(p,i)
=
8o
+
p
l
c
+
8
2
c
2
was
sampled
using
WinBugs
version
1.4
with
all
pa-
rameters
being
given
the
following
diffuse
priors:
Pk
-
N(0,
a
2
=
10
6
).
Table
III
highlights
this
com-
parison
and
shows
how
an
analysis
that
ignores
the
lab-source
heterogeneity
may
improperly
quantify
the
level
of
precision
in
the
analysis.
Although
the
central
estimates/posterior
means
are
similar
for
all
o
k
parameters,
the
credible
intervals
are
much
wider
for
the
hierarchical/mixed
model
estimates.
Employ-
ing
such
naive/pooled
approaches
might
lead
to
risk
estimates
that
underestimate
the
true
nature
of
the
hazard
involved
or
that
report
a
false
level
of
pre-
cision
for
the
estimate
of
the
BMD.
The
hierarchi-
cal
model
is
superior
to
the
pooled
fixed-effects-only
analysis
in
terms
of
describing
the
experimental
data
as
indicated
by
a
smaller
deviance
information
crite-
rion,
(15)
adding
support
for
the
use
of
a
hierarchical
model.
254
Wheeler
and
Bailer
5.
DISCUSSION
The
BMD
based
on
a
hierarchical
model
allows
the
risk
analyst
to
gain
further
insight
into
the
system
being
studied
by
incorporating
data
source
variability
into
the
estimation
of
a
potency
endpoint.
In
the
example
described
above,
the
model's
linear
term
showed
the
highest
degree
of
variability
among
labs.
This
variability
may
simply
indicate
differences
in
lab
conditions.
For
example,
this
might
indicate
differ-
ences
in
the
lab
populations
of
C.
dubia
used
by
the
different
labs.
It
may
be
that
particular
labs
are
using
subpopulations
that
are
more
sensitive
to
NaCl
expo-
sure.
Such
research
questions
are
interesting
in
their
own
right,
and
this
analysis
might
help
identify
labs
for
additional
investigation.
For
example,
one
may
define
a
lab
or
its
C.
dubia
as
"sensitive"
when
its
linear
term
shows
no
benefi-
cial
effect
(i.e.,
those
in
the
population
whose
01
+
b
1
<
0),
and
seek
to
estimate
the
RIp
from
this
group.
This
calculation
can
be
readily
incorporated
in
a
Bayesian
framework
since
each
MCMC
iteration
produces
a
sample
from
the
posterior
distribution
of
the
parameters.
The
RIp
calculation
could
then
be
restricted
to
those
samples
having
01
+
b1
<
0
(and
02
+
b2
<
0).
If
this
is
done,
then
the
subpopulation's
RI25
posterior
mean
is
estimated
at
0.51
with
sym-
metric
95%
credible
interval
being
(0.32,
0.67).
The
RI50
posterior
mean
is
estimated
to
be
0.84
with
a
symmetric
95%
credible
interval
being
(0.61,
1.05).
Not
surprisingly,
the
RIp
mean
estimates
are
lower
in
this
analysis
than
the
previously
derived population-
average
values,
as
they
represent
RIp
estimates
from
lab
sources
whose
C.
dubia
populations
show
no
low-
concentration
increases
in
response
to
the
exposure
of
sodium
chloride.
Alternatively,
an
investigation
into
labs
where
01+
b
l
>
0
might
be
of
interest
since
these
labs
exhibit
some
increase
in
mean
response
at
low
concentrations.
Arguably
a
subpopulation
analysis
should
only
be
used
if
there
is
reasonable
scientific
evidence
that
such
a
population
exists.
Further,
the
motivation
for
this
additional
analysis
will
come
from
risk
managers,
who
may
want
to
consider
how
these
analyses
sup-
port or
influence
a
risk
management
decision.
This
analysis
was
intended
to
highlight
and
illustrate
other
questions
that
might
be
explored
in
the
context
of
the
modeling
framework
under
discussion.
The
hi-
erarchical
structure
allows
such
broad
questions
to
be
asked,
whereas
the
standard
modeling
framework
may
force
risk
managers
to
make
a
decision
based
upon
one
single
lab.
6.
CONCLUSIONS
The
hierarchical
modeling
framework
describes
and
illustrates
an
analytic
framework
for
model-
ing
experimental
data
from
multiple
sources
and
deriving
BMD
estimates
that
reflect
this
lab-to-lab
variability.
This
reflects
an
integration
of
more
sci-
entific
knowledge
into
the
decision
process
and
for-
mally
incorporates
an
important
source
of
variabil-
ity
into
potency
estimation.
Regardless
of
how
the
hierarchical
model
is
fit,
this
strategy
may
more
re-
alistically
represent
the
true
system
by
capturing
lab
source
variation
explicitly.
While
the
framework
was
demonstrated
by
estimating
the
BMD
based
upon
concentration
(Poisson)-response
data,
it
is
easily
generalized
given
a
suitable
definition
of
the
BMD
for
other
responses.
Finally,
we
note
that
the
afore-
mentioned
methods
assume
the
link
function
that
is
a
linear
function
of
the
parameters.
In
practice,
this
may
not
be
the
case
as
one
might
be
interested
in
a
nonlinear
function
of
the
parameters.
In
such
cases,
minimal
suitable
modifications
can
be
made
to
the
above
framework,
and
the
estimation
would
proceed
as
specified
above.
ACKNOWLEDGMENTS
The
authors
would
like
to
thank
Dr.
Walter
Piegorsch,
Dr.
Robert
Noble,
and
Dr.
Elizabeth
Margosches
along
with
two
anonymous
reviewers
and
editors
for
suggestions
that
led
to
improvements
in
this
article.
Benchmark
Dose
Estimation
Incorporating
Multiple
Data
Sources
255
APPENDIX:
WINBUGS
CODE
Full
listing
of
the
WinBugs
code
used
to
obtain
posterior
estimates
from
Equation
(5)
is
given
be-
low.
The
dataset,
which
is
too
large
to
include
in
the
model{
Appendix,
is
specified
by
the
variables
totyoung[],
conc[],
conc2[],
dlines,
and
labid[].
Further
note
that
the
RIp
estimates
were
calculated
from
the
sampled
parameter
values
in
the
Microsoft
Excel
spreadsheet
package.
##set
up
all
of
the
random
effects
for
(j
in
1:7){
##################################
#specify
the
random
effects
#f
or
individual
labs
##################################
lab
.
int
[j]
—dnorm
(lab
.
mu
.
int
,
1
.
tau
.
int)
lab
.
bone
[j]
—dnorm
(lab
.
mu
.
one
,
1
.
tau
.
one)
lab
.
btwo
[j]
—dnorm
(lab
.
mu
.
two
,
1
.
tau
.
two)
}
######################################################
#Specify
the
likelihood
######################################################
for
(i
in
1:dlines){
totyoung[i]—dpois(mu[i])
#################################################
#specify
the
fixed
effects
models
#################################################
log(mu[i])<-
lab.int[labid[i]]
conc[i]
*lab.bone[labid[i]]
conc2Cirlab.btwo[labid[i]]
}
1.tau.one—dgamma(0.001,0.001)
1.tau.two—dgamma(0.001,0.001)
1.tau.int—dgamma(0.001,0.001)
s.one
<-
1/l.tau.one
s.two
<-
1/l.tau.two
s.int
<-
1/l.tau.int
lab.mu.int—dnorm(0,1.0E-6)
lab.mu.one—dnorm(0,
1
.
0E-6)
lab.mu.two"
,
dnorm(0,1
.0E-6)
}
256
Wheeler
and
Bailer
REFERENCES
1.
Bailer
AJ,
Hughes
MR,
Denton
DL,
Oris
JT.
An
empirical
comparison
of
effective
concentration
estimators
for
evaluat-
ing
aquatic
toxicity
test
responses.
Environmental
Toxicology
and
Chemistry,
2000;
19:141-150.
2.
Clayton
DG.
Generalized
linear
mixed
models.
In
Gilks
WR,
Richardson
S,
Spiegelhalter
DJ
(eds).
Markov
Chain
Monte
Carlo
in
Practice.
London:
Chapman
and
Hall,
1996.
3.
Catalano
PJ,
Ryan
RM.
Bivariate
latent
variable
mod-
els
for
clustered
discrete
and
continuous
outcomes.
Jour-
nal
of
the
American
Statistical
Association,
1992;
87:651-
658.
4.
Dunson
DB,
Chen
Z,
Harry
J.
A
Bayesian
approach
for
joint
modeling
of
cluster
size
and
subunit-specific
outcomes.
Bio-
metrics,
2003;
59:521-530.
5.
Piegorsch
WW,
Cox
LH.
Combining
environmental
informa-
tion:
II.
Environmental
epidemiology
and
toxicology.
Envi-
ronmetrics,
1996;
7:309-324.
6.
Simmons
SJ,
Piegorsch
WW,
Nitcheva
D,
Zeiger
E.
Com-
bining
environmental
information
via
hierarchical
modeling:
An
example
using
mutagenic
potencies.
Environmetrics,
2003;
14:159-168.
7.
Coull
BA,
Mezzetti
M,
Ryan
LM.
A
Bayesian
hierarchical
model
for
risk
assessment
of
methylmercury.
Journal
of
Agri-
cultural,
Biological,
and
Environmental
Statistics,
2003;
8:253-
270.
8.
Crump
KS.
A
new
method
for
determining
allowable
daily
intakes.
Fundamental
and
Applied
Toxicology,
1984;
4:854-
871.
9.
Gaylor
DW,
Slikker
W.
Risk
assessment
for
neurotoxicants.
Pp.
341-353
in
Tilson
H,
Mitchell
C
(eds).
Neurotoxicology.
New
York:
Raven
Press,
1992.
10.
Bailer
AJ,
Oris
J.
Modeling
reproductive
toxicity
in
Ceri-
odaphnia
tests.
Environmental
Toxicology
and
Chemistry,
1993;
12:787-781.
11.
Efron
B,
Tibshirani
RI.
An
Introduction
to
the
Bootstrap.
New
York:
Chapman
and
Hall,
1993.
12.
Gilks
WR,
Richardson
S,
Spiegelhalter
DJ.
Introducing
Markov
chain
Monte
Carlo.
In
Gilks
WR,
Richardson
S,
Spiegelhalter
DJ
(eds).
Markov
Chain
Monte
Carlo
in
Prac-
tice.
London:
Chapman
and
Hall,
1996.
13.
Weber
CI,
Peltier
WH,
Norberg-King
TJ,
Horning
WB,
Kessler
FA,
Menkedick
JR,
Neiheisel
TW,
Lewis
PA,
Kemm
K,
Pickering
QH,
Robinson
EL,
Lazorchack
JM,
Wymer
LJ,
Freyberg
RW.
Short-Term
Methods
for
Estimating
the
Chronic
Toxicity
of
Effluents
and
Receiving
Waters
to
Fresh-
water
Organisms.
Cincinnati,
OH:
U.S.
Environmental
Protec-
tion
Agency,
2nd
edition,
EPA/600/4-89/001A,
1989.
14.
Spiegelhalter
D,
Thomas
A,
Best
N.
(2003)
WinBUGS
User
Manual,
Version
1.4.
Cambridge:
MRC
Biostatistics
Unit,
Institute
of
Public
Health.
Available
at:
http://www.mrc-
bsu.cam.ac.uk/bugs
.
Accessed
on
November
20,
2007.
15.
Spiegelhalter
D,
Best
NG,
Carlin
BP,
Van
Der
Linde
A.
Bayesian
measures
of
model
complexity
and
fit.
Journal
of
the
Royal
Statistical
Society
Series
B,
2002;
64:583-639.