1 

DOSE-RESPONSE ANALYSIS OF COMBINED PREPUBERTAL DATA

JOINT MALE AND FEMALE BENCHMARK DOSES

CONTENTS

The Problem .............................................................................................................................................................................. 2

Summary of Results ................................................................................................................................................................ 3

Recommendation ................................................................................................................................................................ 6

Details of the Statistical Analysis ...................................................................................................................................... 7

Statistical Model for the Data ......................................................................................................................................... 7

Dose-Response Shape ....................................................................................................................................................... 9

Approach 1: Estimate BMD for Pooled Blocks 2 and 3 only ........................................................................... 13

Approach 2: Hierarchical Model for All Three Blocks ....................................................................................... 13

References ................................................................................................................................................................................ 14

Appendices ............................................................................................................................................................................... 15

1. Figures Supporting Modeling Decisions ........................................................................................................ 15

residual plots for ML Fit with W fixed to -1 ...................................................................................................... 15

2. Data Preprocessing................................................................................................................................................. 16

3. Code for Bayesian Analyses ................................................................................................................................ 16

i. Fit Separate background, shape, and scale for each block ................................................................ 16

Script (A1script.txt) ............................................................................................................................................... 16

Model (A1.bug)......................................................................................................................................................... 17

Data (data.bug) ........................................................................................................................................................ 17

Initialization .............................................................................................................................................................. 18

ii. Fit just Blocks 2 and 3, separate background but common shape and scale ............................. 19

Script (A2script.txt) ............................................................................................................................................... 19

Model (A2.bug)......................................................................................................................................................... 19

Data (data.bug) ........................................................................................................................................................ 19

Initialization .............................................................................................................................................................. 20

iii. Multilevel Model for Shape and Scale for All Three Blocks, Informative Priors ................. 21
2 

Script (A3script.txt) ............................................................................................................................................... 21

Model (A3.bug)......................................................................................................................................................... 21

Data (data.bug) ........................................................................................................................................................ 22

Initialization .............................................................................................................................................................. 22

4. R Code for Figures and Data Analysis ............................................................................................................ 23

Functions ......................................................................................................................................................................... 23

R Scripts for Analysis .................................................................................................................................................. 27

Preliminaries ............................................................................................................................................................. 27

Dose-Response Analysis ....................................................................................................................................... 30

Benchmark Doses from Blocks 2 and 3.......................................................................................................... 39

Hierarchical Analysis ............................................................................................................................................. 46

THE PROBLEM

There are three blocks of data from the pre-pubertal assays. Two are male-only blocks, and one is 

female-only. The main goal is to estimate benchmark doses for 10% and 20% reduction in T4 

activity for these data. If the data support it, this should be a common set of benchmark doses for 

males and females, otherwise separate sets, each with their own set of profile-likelihood-based 

confidence limits. 

This begs the question, how to determine whether the male-only and female-only data sets support 

a common dose-response. The factor of interest, gender, is almost perfectly confounded with 

blocking. Thus, the magnitude and nature of differences among blocks affect the interpretation of 

any observed difference between the genders. In principle, this should boil down to asking whether 

the difference between male- and female-specific benchmark doses is large relative to the 

difference between benchmark doses computed for the two male-only blocks. However, there is 

the additional complication that one of the male-only blocks has only two non-control dose groups, 

so we cannot fit the preferred hill model (four parameters) to it: 
3 

Figure 1 Doses and sample sizes in the pre-pubertal dataset. Blocks 1 and 2 are male-only, block 3 is female only. 

So, Block 2 (males) has both the smallest and largest non-control doses, while Block 1 (males) has 

no low doses at all. As a consequence, Block 1 data have no information about the behavior of the 

dose-response below 100 mg/kg. Furthermore, only the control dose level appears in more than 

one block, so the three blocks cannot be compared for internal consistency by just comparing the 

responses at the same dose levels. The only way to compare across blocks is to compare doseresponse functions fit to each block. The small number of dose levels in block 1 means that either 

relatively simple models, with no more than three parameters, must fit the data, or that block must 

be left out of comparisons that require more complex models. So, an initial part of the analysis for

these data must be to identify the simplest form of model that will fit all three blocks. Next, it would 

simplify the analysis if the same model shape parameters would yield good fits for both the male 

blocks. If so, then the final step would be to see if the same model shape fits both sexes. If, 

however, the two male blocks differ, then we need to determine why, and evaluate how that affects 

the comparison with the single female block. All code for this analysis, carried out in R, is included 

in the Appendices. 

SUMMARY OF RESULTS

The Hill model fits each block quite well. A lognormal multiplicative error model describes these 

data adequately, and the coefficient of variation is the same across all three blocks. The term in the 

Hill model that controls the asymptotic value of the response as dose increases was fixed, so that 

the level of T4 would go to 0 as dose increases. This was done so that block 1 could be fit by the 

same dose-response model used for the other two blocks. As it turns out, constraining the model in 

this way has little effect on its fit to the data. 

Dose (mg/kg/day)

Block

1

2

3

0 50 100 150 200 250 300

Sample

Size

7

8

9

10

15
4 

The two male blocks are very significantly different from each other, while the dose-response for 

the male block that included low doses (block 2) is very similar to the female block (P = 0.72) 

(Figure 2). 

Figure 2 Estimates of T4 levels in the three pubertal assay blocks, relative to block-specific controls. Error bars 

are 95% confidence intervals for the ratios of dose-specific level to the corresponding block-specific control, 

assuming log-normality and the same coefficient of variation across all blocks and doses (which assumptions is 

supported by the data). 

I considered two interpretations of this result: 

1) “ ‘Something’ went wrong with block 1, causing it to differ from the other two blocks. Since, 

in any case, it only includes two non-control doses, where the response is already quite a bit 

greater than that we are considering for benchmark doses, its data are nearly irrelevant to 

the problem at hand. Since the remaining two blocks are so similar, combine their data, and 

use only them to estimate BMDs.” 

2) “We have no direct information about what, if anything, might have gone wrong with block 

1. The variation among the three blocks may just reflect normal or typical variation among 

experimental blocks, due to minor differences in husbandry, reagents, or experimental 

procedures. In fact, the difference between blocks 1 and the combined blocks 2 and 3, while 

significant, does not appear large at the response level. It makes sense to base benchmark 

dose estimates on characteristics of the “population” of such studies, and base it on a 

statistically valid estimate of the average dose-response.” 

In either case, since gender seems to make no difference to the dose-response, both interpretations 

lead to an analysis that ignores gender. 

Under the first interpretation, the Hill model was fit to the data for blocks 2 and 3 (Figure 3). In this 

model, there is a separate background level of each block, but only a single shape (n), scale (K), and

error standard deviation (sigma) parameter. This yields benchmark doses of: 
5 

Table 1 BMDs and 90% confidence limit (thus, the BMDL corresponds to the conventional 95% one-sided limit 

used for the BMDL) for the analysis of blocks 2 and 3 only, using maximum likelihood methods. Units for BMD 

and confidence limits are mg/kg/day 

Response (BMR) BMD (mg/kg/day) BMDL BMDU 

0.10 3.41 1.49 7.07 

0.20 9.43 4.95 16.67 

Figure 3 Data from blocks 2 and 3 (block 2 contains only males, while block 3, only females), and fitted Hill 

model. Data are expressed here as T4 levels relative to block-specific control (as in Figure 2). 

Under the second interpretation, the appropriate course of action is to fit a multilevel model to the 

data, and base BMD estimates on the parameters for the upper level. It is difficult to estimate 

parameters of such models using maximum likelihood methods, but Bayesian methods lend 

themselves quite handily to multilevel models. Even so, I was unable to reliably estimate model 

parameters, though, without adding additional information to the problem, in the form of 

informative priors. Specifically, I assumed that it was unlikely that the block- or experimentspecific Hill coefficient (n)would fall below 1/3 or above 3 times the population mean value for the 

Hill coefficient in a replicate block, and that the population mean value itself would fall between 0.1 

and 10. Further, it was assumed that the scale parameter, K would mostly fall between 5% and 

95% of the population median, and the population median most likely falls between 2% and 98% of 

the maximum dose. The mean dose-response curves under these assumptions are shown in Figure 

4. The distribution of BMDs is quite skewed (the mean BMD10 is nearly 200000, while the median is 

10.37), so the median is used here as the measure of central tendency. The central estimates of the 

BMDs are larger than for just blocks 2 and 3, as are the BMDLs, while the BMDUs are quite a bit 

greater. 
6 

Table 2 BMDs and confidence limits from the Bayesian analysis of the hierarchical model using all three blocks. 

The BMD is the median of the posterior distribution and BMDL and BMDU are 5 and 95 percentiles of the 

posterior distribution. Units are mg/kg/day. 

Response (BMR) BMD (mg/kg/day) BMDL BMDU 

0.10 10.37 3.47 35.18 

0.20 19.05 8.07 46.85 

Figure 4 The T4 data, as in Figure 2, with mean dose-response curves from a Bayesian analysis. Colored curves 

are block-specific estimates, the solid black curve is from the median and the dashed black curves the 5 and 95 

point-wise percentiles of the posterior distribution of the population mean response. 

The specifics of the informative prior change the answer somewhat; less restrictive priors increase 

the mean estimate, while dropping BMDLs to lower levels and increasing the BMDU. 

RECOMMENDATION

It is difficult to arrive at a definitive recommendation between these two analyses. After discussion

with the investigators who generated the data, it appears there is nothing exceptional about the 

first block other than it was a pilot study. That makes it difficult to disqualify those data on the 

basis of intrinsic characteristics of the experiments themselves. On the other hand, that block has 

no doses with responses close to the benchmark response, so its dose-response shape is totally 

governed by the behavior of two high doses. Small changes in those responses can have large 

effects on the BMD estimates. In effect, taking that block very seriously amounts to saying that the 

Hill model really captures the true dose-response. In balance, I prefer the analysis that drops block

1. This also skirts the potential objections that the Bayesian analysis uses a somewhat ad hoc 

informative prior. That the BMDLs arrived at by either method are not so different makes this 

decision somewhat easier to make. 
7 

DETAILS OF THE STATISTICAL ANALYSIS

STATISTICAL MODEL FOR THE DATA

The first step in the analysis of this or any dose-response data set is to establish a plausible 

statistical model for the data. Practically, initially this comes to deciding how best to make 

variances homogeneous, or if that is not possible, whether a power model for the variance will 

work. First, use the Box-Cox transformation to identify a power transformation to homogenize the 

variances. 

Figure 5 Box-cox estimated power for making variances homogeneous. This suggests a log-transformation 

(corresponding to lambda = 0). 

We can see this by plotting both dose- and block-specific standard deviations and geometric 

standard deviations against dose. The standard deviations decline with dose, while the geometric 

standard deviations are fairly constant (geometric standard deviations are related to coefficients of 

variation and to the standard deviation of the log transformed data). 

-2 -1 0 1 2

-350

-300

-250

-200

λ

log-Likelihood

 95%
8 

Figure 6 Standard deviations (solid symbols) and geometric standard deviations (10sd(log10(T4))) (open symbols) 

by dose and block. This result, constancy of the geometric standard deviation, is consistent with log-normality of 

the data. 

Since this study was conducted in blocks, one likely possibility is that, while coefficients of variation 

may be similar within block, they may differ among blocks. Figure 3 gives scant support for this 

idea; the dose- and block-specific estimates of geometric standard deviation are similar among 

blocks. This question can be formally tested by comparing the log-likelihood of a model fitting a 

separate mean to each dose X block group and a separate coefficient of variation for each block, to 

that log-likelihood of the same model for the means, but a common coefficient of variation. The Pvalue for this comparison is about 0.14, so the data should be well-described by a model with a 

common coefficient of variation. 

Finally, the quantile-quantile plot of residuals from dose- and block-specific means, on the log 

scale. These should form a straight line if log transformation is appropriate. 

Dose (mg/kg/day)

x

0.4

0.6

0.8

1.0

1.2

1.4

0 50 100 150 200 250 300

factor(Block)

1

2

3
9 

Figure 7 Quantile-quantile plot of log10 transformed data, centered on doseXblock means and scaled by a robust 

estimate of its standard deviation, versus theoretical quantiles of the normal distribution. The straight line is 

consistent with a lognormal distribution of the original data. 

All this says that it is reasonable to assume the data are lognormally distributed about dose- and 

block-specific medians, with a constant coefficient of variation. That is 

log ( 4) ~ , 10 , T N (μ σ i j )

The subscripts i and j index dose and block, respectively. 

DOSE-RESPONSE SHAPE

It would certainly be convenient if log10(T4) were linear in dose. That would imply a simple 

exponential model for the dose-response, which could be estimated even for Block 1. However this 

is unlikely: 

theoretical

sample

-2

-1

0

1

2

-2 -1 0 1 2
10 

Figure 8 Means and confidence intervals (computed on the log scale) for T4 in the three blocks, along with fitted 

lines. 

It looks as if, in blocks 2 and 3, the initial drop has a greater slope, and the final (high-dose) drop is 

more gradual, than is consistent with a straight line. All three means for block 1 are pretty much 

right on their lines. Next, look at the residuals from this fit. That should show the patterns more 

clearly. 

Figure 9 Mean residuals from fits shown in Fig. 5. 

The dose-response needs to be concave. In decreasing data, this can be accomplished in both the 

hill model and the exponential model by having a power parameter less than 1, or/and by a 

Dose (mg/kg/day)

log10(T4)

-0.2

0.0

0.2

0.4

0.6

0.8

0 50 100 150 200 250 300

Block

1

2

3

Dose (mg/kg/day)

Residuals from linear fit

-0.2

-0.1

0.0

0.1

0.2

0.3

0 50 100 150 200 250 300

Block

1

2

3
11 

horizontal asymptote greater than 0. Background levels probably differ among blocks (the P-value 

for the comparison is 0.0461), with the block 3 control being about 91% of the block 1 control 

value. 

Now, the modeling switches to maximum likelihood fitting, using a lognormal error model. Since 

there are often only small differences between general exponential model and hill model fits, and 

there is already a history of fitting these data with the hill model, the hill model is used to estimate 

BMDs for these data. However, initially I will turn it into a three-parameter model in two different 

ways: either fix the hill coefficient n to 1 and continue to estimate a horizontal asymptote, or 

estimate n and set the horizontal asymptote to 0. The model with the better log-likelihood from 

these fits is to be preferred, as long as its likelihood compares favorably to the likelihood from a 

saturated model, in which each group is characterized by a separate mean, whose deviance (that is, 

-2 × log-likelihood) is 192.14, with 13 model degrees of freedom. 

The best-fitting model with n fixed to 1 has a deviance of 202.2, and 10 model degrees of freedom, 

yielding a p-value of 0.018, that is, its fit to the data is significantly worse than the fully saturated 

model. In contrast, by estimating n, but setting w to -1 (thus, fitting models that always go to 0 T4

activity at sufficiently high doses), the deviance is 192.97 with 10 model degrees of freedom. The 

p-value for the fit of this model to the data is 0.84, indicating a quite adequate fit to the data (this 

model is the point of comparison of all subsequent models fit to these data in this report). 

Residuals from this fit are shown in figures 12 and 13, Appendix 1. Figure 10 shows a direct 

comparison of fit and data: 

Figure 10 Means with 95% confidence intervals, with the fitted curves superimposed. 

The low-dose shape of the curve for block 1 differs quite a bit from the shapes in blocks 2 and 3. 

The next figure shows the estimate of n (shape) and K (scale) for each block, with associated 

contours characterizing the posterior distribution of each from a Bayesian analysis of this same 

Dose (mg/kg/day)

T4

1

2

3

4

5

6

7

0 50 100 150 200 250 300

Block

1

2

3
12 

model (see Appendix 2, “Bayesian Analysis” for more information about this analysis). Blocks 2 and 

3 have very similar values, while those for block 1 are quite distinct from the other two. 

Figure 11 Block-specific estimates of n and K for the hill model, and contours of the posterior distribution of each 

estimate (see Appendix 2). The heavier contours are 95% confidence contours for the estimates, based on a 

bivariate normal approximation to the posterior density. 

If the three blocks had similar shapes and scales (but possibly different background levels), then a 

common benchmark dose could be estimated from a curve fit to the data with common shape and 

scale parameters. However, such a model has a significantly worse fit to the data than the reference 

model (P = 0.0065). Blocks 1 and 2 differ significantly (P = 0.0019), but blocks 2 and 3 do not (P =

0.72), consistent with Figure 8. Given all this, there seem to be three possible approaches to 

computing BMDs for these data: 

1. Block 1 stands out from the other two, and it was a pilot study. Furthermore, since block 1 

has no lower doses, the BMD estimates from that block are totally determined by responses 

at higher doses. Thus, we should ignore that block, and combine blocks 2 and 3 and use 

them to estimate the BMD, using profile likelihood methods to compute the BMDLs (as does 

the Agency software package BMDS). 

2. Fit a multilevel model to the data, with separate shape and scale parameter for each block, 

and a top-level overall estimate. The top-level estimate corresponds to the pooled estimate. 

This is a statistically supportable approach to estimating a pooled BMD for the whole 

dataset, but would be less valid if it were determined that block 1 was somehow processed 

differently from blocks 2 and 3, and so represented a different dose-response. 

I will work through these options in what follows, using both maximum likelihood and Bayesian 

methods for the first approach, and only Bayesian methods for the second because it is 

computationally substantially more difficult to use maximum likelihood for that. Code for all these 

analyses is given in the Appendices. 

K

n

0.5

1.0

1.5

2.0

2.5

20 40 60 80 100 120

block

1

2

3
13 

APPROACH 1: ESTIMATE BMD FOR POOLED BLOCKS 2 AND 3 ONLY

Both the maximum likelihood method applied in the previous section (and whose code is 

reproduced in Appendix 4) and Bayes approach (whose code is reproduced in Appendix 3) gave 

similar estimates for BMD and BMDL, as shown in the following table (all doses are expressed in 

units of mg/kg/day): 

BMR Maximum Likelihood Bayes 

BMD BMDL BMDU BMD BMDL BMDU 

0.10 3.41 1.49 7.07 4.01 1.57 7.63 

0.20 9.44 4.95 6.70 10.60 5.18 17.80 

It was shown in the previous section that the Hill model gives a very adequate fit to the full dataset. 

Furthermore, Figure 3 shows the fit to blocks 2 and 3 to be quite adequate. The point of the Bayes 

analysis here is to build confidence in the Bayes analysis, since that is the methodology used in 

approach 2. 

APPROACH 2: HIERARCHICAL MODEL FOR ALL THREE BLOCKS

The code for the hierarchical approach is duplicated in Appendix 3.iii. 

To include all three blocks in the analysis, we assume that each block is a separate experiment, and 

that the results of the three blocks in hand are samples from a population of possible blocks, whose 

dose-response curves all differ from each other (due to analytic variability, slightly different 

samples of animals, or other unrecorded differences). We assume that the Hill model characterizes 

the dose-response curves for all the blocks, so the difference among blocks is in the parameters of 

the Hill model: each block has its own set of parameters, sampled from a distribution of parameters. 

In addition to the assumptions made in the maximum likelihood model (the dose-response is a Hill 

model, errors are lognormal and each block has the same coefficient of variation, the background 

level of T4 varies among blocks), we assume that the logarithms of the Hill coefficient, n and scale 

parameter K, are sampled from a multivariate normal distribution with population mean (lnk.mu) 

and population variance matrix (lnk.var). Together these two comprise the population parameters. 

The presence of this second level is the reason for calling this a hierarchical approach. 

Preliminary work with this model indicated that there were not enough blocks to estimate the 

population parameters without additional restrictions. Consequently, priors were set on the 

population parameters to reflect plausible bounds on the distribution of the Hill model parameters. 

It is important to remember that subsequent BMD estimates and, especially, confidence limits are 

sensitive to these priors. First of all, the logarithm (only natural logarithms are used in this 

analysis) of the Hill coefficient was assigned a uniform prior on the interval (log(.1), log(10)). That 

is, an interval centered at 1, with the possibility of ten-fold variation on either side. This is derived 

from the practical observation that Hill coefficients rarely deviate much from 1 in practical 

application. Prior to analysis, the dose variable was rescaled so the maximum dose (300 mg/kg in 

practice) is 1. This allows a reasonably prior on log(K) (on this rescaled level) to be set as uniform 

on the interval (log(0.02), log(0.98)). This bounds the scale parameter (also the dose value at 

which the response is half its maximum value), essentially to fall within the range of the 

administered doses. Priors on the variance parameters were set so that 95% of sampled values of 
14 

the Hill coefficient would fall between .33 and 3, and of K would fall between .05 and .95 (see code 

in Appendix 3 for ‘A3.bug’ for details). 

Three chains were run out for 100000 samples, thinning to every tenth sample, using JAGS version 

2.1.0 (Plummer, 2010). Convergence of the chains was assessed using the Heidelberger and Welch 

diagnostic, as implemented in the R package coda (Plummer et al, 2010), which indicated that all 

three chains had passed the default diagnostic criteria, and Gelman and Rubin’s convergence 

diagnostic, also as implemented in coda, which gives the multivariate potential scale reduction 

factor for the full set of chains as 1.00. 

BMDs and BMDLs were computed respecting the hierarchical structure of the model parameters. 

At the top level, values for the population mean of log n and log K, and the entries for the population 

covariance matrix are sampled from the estimated posterior distribution. Given those values, 

enough samples are drawn for log(n) and log(K) from a multivariate normal distribution that the 

maximum of the two standard errors is less than 0.01. The mean BMD10 and BMD20 over all the 

sampled parameters is then computed, and a new set of population parameters drawn. In this way, 

BMDs for each sampled set of population parameters is computed. The median of this set of BMDs 

was used for the point estimate for BMD, since the distribution is quite skewed, with a few 

extremely large values which distort the mean: the mean BMD10 is 215865, while the median is just 

10.4 (with 5th and 95th percentiles of 3.5 and 35.3). The resulting benchmark doses and confidence 

limits are reported on page 6. 

REFERENCES

Martyn Plummer (2010). JAGS Version 2.1.0 User Manual. http://wwwfis.iarc.fr/~martyn/software/jags/ 

 Martyn Plummer, Nicky Best, Kate Cowles and Karen Vines (2010). coda: Output analysis and 

diagnostics for MCMC. R package version 0.13-5. http://CRAN.R-project.org/package=coda 

R Development Core Team (2010). R: A language and environment for statistical computing. R 

Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.Rproject.org. 
15 

APPENDICES

1. FIGURES SUPPORTING MODELING DECISIONS

RESIDUAL PLOTS FOR ML FIT WITH W FIXED TO -1 

Figure 12 Mean and 1 standard deviation range for residuals (computed on the log scale: log10(sT4) – 

log10(predicted value)) from maximum likelihood fit with w fixed to -1, plotted against predicted value. This is 

for T4 scaled so the grand mean of the control values is 1. 

Figure 13 Mean and 1 standard deviation range for residuals from maximum likelihood fit with w fixed to -1, 

plotted versus dose. 

pred

res

-0.4

-0.2

0.0

0.2

0.4

0.2 0.4 0.6 0.8 1.0

Block

1

2

3

dose

res

-0.4

-0.2

0.0

0.2

0.4

0 50 100 150 200 250 300

Block

1

2

3
16 

2. DATA PREPROCESSING 

Code for preparing the data frames used in the subsequent analyses from the data as supplied in 

spreadsheets: 

library(gdata) 

## First, get the data sets 

## Males: 

infileM <- file.path("..","Triclosan-MalePubertal",

 "T4concentrationsTriclosan.xls") 

dtaM <- read.xls(infileM, sheet=2, header=TRUE) 

## add gender: 

dtaM$sex <- factor(rep("M", nrow(dtaM)), levels=c("F","M")) 

## Females: 

infileF <- file.path("..","T4dataOnly.xls") 

dtaF <- read.xls(infileF, sheet=1, skip=3) 

## rename the variables 

names(dtaF) <- c("animal","fdose","T4") 

## Summary lines are interspersed among the data lines. 

## This time, it is the 'fdose' variable that is missing for 

## the summary rows; this is a "". 

dtaF <- subset(dtaF, !(fdose %in% c(""))) 

dtaF$dose <- as.numeric(gsub("Control","0.0", as.character(dtaF$fdose))) 

dtaF$sex <- factor(rep("F", nrow(dtaF)), levels=c("F","M")) 

dtaF$block <- rep(3, nrow(dtaF)) 

## stack the two datasets (drop the animal variable from dtaF) 

dta <- rbind(dtaM[,c("T4","dose","block", 

"sex")],dtaF[,c("T4","dose","block","sex")]) 

## Save the data frame for later use 

save(dta, file="CombinedPubertal.RData") 

3. CODE FOR BAYESIAN ANALYSES

Both OpenBUGS version 3.1.1 rev 524, and JAGS version 2.1.0 were used for these computations. 

The code reported here is for JAGS, to make it easier to report repeatable results. 

i. FIT SEPARATE BACKGROUND, SHAPE, AND SCALE FOR EACH BLOCK

Script (A1script.txt) 

model in A1.bug 

data in data.bug 

compile, nchains(3) 

parameters in parms1.R, chain(1) 

parameters in parms2.R, chain(2) 

parameters in parms3.R, chain(3) 

initialize 

update 3000 

monitor sigma, thin(10) 

monitor A, thin(10) 

monitor n, thin(10) 

monitor K, thin(10) 

monitor BMD10, thin(10) 

monitor BMD20, thin(10) 
17 

update 100000 

coda *, stem(Full_) 

Model (A1.bug) 

model 

{ 

 sigma ~ dunif(0,10) 

 tau <- 1/(sigma * sigma) 

 for (i in 1:3) { 

 W[i] <- -1 

 lA[i] ~ dt(0,1,2) 

 A[i] <- exp(lA[i]) 

 K[i] ~ dunif(0,2) 

 ln[i] ~ dt(0,1,2) 

 n[i] <- exp(ln[i]) 

 BMD10[i] <- pow(0.1/0.9, 1/n[i]) * K[i] * 300 

 BMD20[i] <- pow(0.2/0.8, 1/n[i]) * K[i] * 300 

 } 

 for (i in 1:N) { 

 g0[i] <- 1 - step(-dose[i]) 

 ddose[i] <- (1 - g0[i]) + dose[i] 

 g1[i] <- step(dose[i] - K[block[i]]) 

 mu[i] <- log(A[block[i]] * ( (1 - g0[i]) + 

 g0[i] * 

 (g1[i] * 

 (1 + W[block[i]]/(pow(K[block[i]]/ddose[i],n[block[i]]) + 1)) 

+ 

 (1 - g1[i]) * 

 (1 + W[block[i]]*pow(ddose[i]/K[block[i]], n[block[i]]) / 

 (1 + pow(ddose[i]/K[block[i]], n[block[i]])))))) 

 T4[i] ~ dlnorm(mu[i], tau) 

 } 

} 

Data (data.bug) 

'T4' <- c(1.06879536664825, 1.05747870982492, 1.02918706776658, 

1.46508792318391, 1.21339709272419, 1.14738326125474, 1.18280020575740, 

1.14298233915678, 0.99984758711349, 0.968621996989845, 1.24650879231839, 

1.03002533864238, 0.98978833660386, 0.875783497494713, 0.92733715635657, 

0.537960334546286, 0.346834574863305, 0.475090018861095, 0.560174512755053, 

0.43401474594677, 0.445540970489055, 0.680047247994818, 0.467545580978872, 

0.244565528015394, 0.261750080969346, 0.170378555506868, 0.145020861513841, 

0.132446798376803, 0.280820743393854, 0.274743279544285, 0.237440225571072, 

0.99942845167559, 0.820248051972794, 1.10044009220980, 1.36868677246661, 

1.06292747051763, 0.837223037207796, 0.675017622740003, 0.75947341347711, 

0.987063956257502, 0.820248051972794, 0.793842519385014, 1.16058602754863, 

1.08451294556955, 1.16813046543085, 0.732648745451428, 0.798033873764027, 

0.614242984244318, 0.703099697079388, 1.13921012021566, 0.788184190973347, 

0.500447712854122, 0.468802987292576, 0.666215778544076, 0.681723789746423, 

0.475299586580045, 0.728457391072415, 0.576730362552154, 0.741450589647355, 

0.586999180780735, 0.176246451637486, 0.168492446036312, 0.176246451637486, 

0.181066509173351, 0.224656594715083, 0.189239650212425, 0.268037112537865, 

0.29926270266151, 0.165348930252053, 0.110442187886986, 0.538169902265237, 

0.770999638019394, 0.827373354417116, 0.709805864085808, 1.51370763398045, 
18 

1.26453161614815, 0.820876755129646, 0.778753643620568, 0.534607251043076, 

1.02688182285812, 0.469012555011526, 1.02310960391701, 0.432547771914115, 

0.711691973556364, 0.485358837089676, 0.598734973041971, 0.404884833012631, 

0.514279182304864, 0.570024195545733, 0.362761721503553, 0.642325058583703, 

0.594753186381909, 0.316866391053364, 0.601459353388329, 0.373449675170036, 

0.508830421612147, 0.327973480157747, 0.321267313151327, 0.419345005620225, 

0.210615557545390, 0.379527139019604, 0.29926270266151, 0.265941435348359, 

0.211453828421193, 0.23681152241422, 0.209777286669588, 0.343691059079045, 

0.348930252052811, 0.254205643087123) 

'dose' <- c(0, 0, 0, 0, 0, 

0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.333333333333333, 0.333333333333333, 

0.333333333333333, 0.333333333333333, 0.333333333333333, 0.333333333333333, 

0.333333333333333, 0.333333333333333, 0.666666666666667, 0.666666666666667, 

0.666666666666667, 0.666666666666667, 0.666666666666667, 0.666666666666667, 

0.666666666666667, 0.666666666666667, 0, 0, 0, 0, 0, 0, 0, 0, 

0, 0, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 

0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 1, 1, 1, 1, 1, 1, 

1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0.03125, 0.03125, 0.03125, 

0.03125, 0.03125, 0.03125, 0.03125, 0.03125, 0.125, 0.125, 0.125, 

0.125, 0.125, 0.125, 0.125, 0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 

0.25, 0.25, 0.25, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5) 

'block' <- c(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, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 

2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 

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

3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 

3, 3, 3) 

'N' <- 109 

Initialization 

Parms1.R 

"sigma" <- 

0.15 

"lA" <- 

c(0.00, -0.18, -0.44) 

"K" <- 

c(0.28, 0.28, 0.28) 

"ln" <- 

c(-0.56, -0.56, -0.56) 

Parms2.R 

sigma <- 0.16 

lA <- c( 0.08, -0.07, -0.18) 

K <- c(0.14, 0.14, 0.14) 

ln <- c(0.15, 0.15, 0.15) 

Parms3.R 

sigma <- 0.30 

lA <- c( 0.15, 0.04, 0.05) 

K <- c(0.06, 0.06, 0.06) 

ln <- c(2.18, 2.18, 2.18) 
19 

ii. FIT JUST BLOCKS 2 AND 3, SEPARATE BACKGROUND BUT COMMON SHAPE AND SCALE

Script (A2script.txt) 

model in A2.bug 

data in data.bug 

compile, nchains(3) 

parameters in parms1.R, chain(1) 

parameters in parms2.R, chain(2) 

parameters in parms3.R, chain(3) 

initialize 

update 3000 

monitor sigma, thin(10) 

monitor A, thin(10) 

monitor n, thin(10) 

monitor K, thin(10) 

monitor BMD10, thin(10) 

monitor BMD20, thin(10) 

update 100000 

coda *, stem(Pooled_) 

Model (A2.bug) 

model 

{ 

 sigma ~ dunif(0,10) 

 tau <- 1/(sigma * sigma) 

 for (i in 1:2) { 

 lA[i] ~ dt(0,1,2) 

 A[i] <- exp(lA[i]) 

 } 

 W <- -1 

 K ~ dunif(0,2) 

 ln ~ dt(0,1,2) 

 n <- exp(ln) 

 BMD10 <- pow(0.1/0.9, 1/n) * K * 300 

 BMD20 <- pow(0.2/0.8, 1/n) * K * 300 

 for (i in 1:N) { 

 g0[i] <- 1 - step(-dose[i]) 

 ddose[i] <- (1 - g0[i]) + dose[i] 

 g1[i] <- step(dose[i] - K) 

 mu[i] <- log(A[block[i]-1] * ( (1 - g0[i]) + 

 g0[i] * 

 (g1[i] * 

 (1 + W/(pow(K/ddose[i],n) + 1)) + 

 (1 - g1[i]) * 

 (1 + W*pow(ddose[i]/K, n) / 

 (1 + pow(ddose[i]/K, n)))))) 

 T4[i] ~ dlnorm(mu[i], tau) 

 } 

} 

Data (data.bug) 

dose <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.01, 

0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.1, 0.1, 

0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
20 

1, 0, 0, 0, 0, 0, 0, 0, 0, 0.03125, 0.03125, 0.03125, 0.03125, 

0.03125, 0.03125, 0.03125, 0.03125, 0.125, 0.125, 0.125, 0.125, 

0.125, 0.125, 0.125, 0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 

0.25, 0.25, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5) 

T4 <- c(0.99942845167559, 

0.820248051972794, 1.10044009220980, 1.36868677246661, 1.06292747051763, 

0.837223037207796, 0.675017622740003, 0.75947341347711, 0.987063956257502, 

0.820248051972794, 0.793842519385014, 1.16058602754863, 1.08451294556955, 

1.16813046543085, 0.732648745451428, 0.798033873764027, 0.614242984244318, 

0.703099697079388, 1.13921012021566, 0.788184190973347, 0.500447712854122, 

0.468802987292576, 0.666215778544076, 0.681723789746423, 0.475299586580045, 

0.728457391072415, 0.576730362552154, 0.741450589647355, 0.586999180780735, 

0.176246451637486, 0.168492446036312, 0.176246451637486, 0.181066509173351, 

0.224656594715083, 0.189239650212425, 0.268037112537865, 0.29926270266151, 

0.165348930252053, 0.110442187886986, 0.538169902265237, 0.770999638019394, 

0.827373354417116, 0.709805864085808, 1.51370763398045, 1.26453161614815, 

0.820876755129646, 0.778753643620568, 0.534607251043076, 1.02688182285812, 

0.469012555011526, 1.02310960391701, 0.432547771914115, 0.7116919733556364, 

0.485358837089676, 0.598734973041971, 0.404884833012631, 0.514279182304864, 

0.570024195545733, 0.362761721503553, 0.642325058583703, 0.594753186381909, 

0.316866391053364, 0.601459353388329, 0.373449675170036, 0.508830421612147, 

0.327973480157747, 0.321267313151327, 0.419345005620225, 0.210615557545390, 

0.379527139019604, 0.29926270266151, 0.265941435348359, 0.211453828421193, 

0.23681152241422, 0.209777286669588, 0.343691059079045, 0.348930252052811, 

0.254205643087123) 

block <- c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 

2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 

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

3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 

3, 3, 3, 3) 

N <- 78 

Initialization 

parms1.R 

"sigma" <- 

0.15 

"lA" <- 

c(-0.18, -0.44) 

"K" <- 

0.28 

"ln" <- 

-0.56 

parms2.R 

sigma <- 0.16 

lA <- c( -0.07, -0.18) 

K <- 0.14 

ln <- 0.15 

parms3.R 

sigma <- 0.30 

lA <- c(0.04, 0.05) 

K <- 0.06 

ln <- 2.18 
21 

iii. MULTILEVEL MODEL FOR SHAPE AND SCALE FOR ALL THREE BLOCKS, INFORMATIVE 

PRIORS

Script (A3script.txt) 

model in A3.bug 

data in data.bug 

compile, nchains(3) 

parameters in parms1.R, chain(1) 

parameters in parms2.R, chain(2) 

parameters in parms3.R, chain(3) 

initialize 

update 3000 

monitor sigma, thin(10) 

monitor A, thin(10) 

monitor n, thin(10) 

monitor K, thin(10) 

monitor lnK.mu, thin(10) 

monitor lnK.var, thin(10) 

monitor BMD10, thin(10) 

monitor BMD20, thin(10) 

update 100000 

coda *, stem(s2_) 

Model (A3.bug) 

model 

{ 

 ## Previous chains, using very broad priors for the variances and 

 ## population means gave surprisingly broad posterior distributions. 

 ## So, see what happens if we use a more informative prior. 

 ## hyperparameters for lnK.tau 

 ## First of all, the variance for the population distribution of 

 ## log(n) (lnK.tau.tau[1,1]) Most likely, this would fall between 

 ## say, .33 and 3 for most (95%) of the experiments. This gives 

 ## a value of .314 for a guess for the variance of log(n) 

 ## For log(K), we have normalized doses so max(K) is 1, so, 

 ## most (95%) of the time, K should fall between, say, .05 and .95 

 ## This gives an estimate of the variance of 0.0428. 

 lnK.tau.tau[1,1] <- 0.314; lnK.tau.tau[1,2] <- 0 

 lnK.tau.tau[2,1] <- 0; lnK.tau.tau[2,2] <- 0.0428 

 ## The population mean of log n is unlikely to fall outside, say 

 ## log(.1) to log(10). Use a uniform prior f 

 lnK.mu[1] ~ dunif(-2.3026, 2.3026) 

 ## The population mean of log(K) is unlikely to fall outside, say 

 ## log 0.02 and log 0.98. Set those as limits for a uniform prior 

 lnK.mu[2] ~ dunif(-3.912023, -0.02020271) 

 sigma ~ dunif(0,10) 

 tau <- 1/(sigma * sigma) 

 n.mu <- exp(lnK.mu[1]) 

 K.mu <- exp(lnK.mu[2]) 

 BMD10 <- pow(0.1/0.9, 1/n.mu) * K.mu * 300 

 BMD20 <- pow(0.2/0.8, 1/n.mu)*K.mu * 300 
22 

 lnK.tau[1:2,1:2] ~ dwish(lnK.tau.tau[,], 2) 

 lnK.var[1:2,1:2] <- inverse(lnK.tau[1:2,1:2]) 

 for (i in 1:3) { 

 W[i] <- -1 

 lA[i] ~ dt(0,1,2) 

 A[i] <- exp(lA[i]) 

 lnK[i,1:2] ~ dmnorm(lnK.mu[], lnK.tau[,]) 

 K[i] <- exp(lnK[i,2]) 

 n[i] <- exp(lnK[i,1]) 

 } 

 for (i in 1:N) { 

 g0[i] <- 1 - step(-dose[i]) 

 ddose[i] <- (1 - g0[i]) + dose[i] 

 g1[i] <- step(dose[i] - K[block[i]]) 

 mu[i] <- log(A[block[i]] * ( (1 - g0[i]) + 

 g0[i] * 

 (g1[i] * 

 (1 + W[block[i]]/(pow(K[block[i]]/ddose[i],n[block[i]]) + 1)) 

+ 

 (1 - g1[i]) * 

 (1 + W[block[i]]*pow(ddose[i]/K[block[i]], n[block[i]]) / 

 (1 + pow(ddose[i]/K[block[i]], n[block[i]])))))) 

 T4[i] ~ dlnorm(mu[i], tau) 

 } 

} 

Data (data.bug) 

Uses the same data as the first entry. 

Initialization 

(parms1.R) 

"sigma" <- 

0.15 

"lA" <- 

c(0.00, -0.18, -0.44) 

"lnK.mu" <- 

c(0.0,-1.61) 

"lnK.tau" <- 

structure(c(11.111111,0,0,100), .Dim=c(2,2)) 

"lnK" <- 

structure(c(0.095, -0.22, -0.11, -1.90, -1.61, -1.51), .Dim=c(3,2)) 

(parms2.R) 

sigma <- 0.16 

lA <- c( 0.08, -0.07, -0.18) 

lnK.mu <- c(-.3, -2.30) 

lnK.tau <- structure(c(4,0,0,1.5625), .Dim=c(2,2)) 

(parms3.R) 

sigma <- 0.30 

lA <- c( 0.15, 0.04, 0.05) 

lnK.mu <- c(-0.5, -0.51) 
23 

lnK.tau <- structure(c(100,0,0,11.111111), .Dim=c(2,2)) 

4. R CODE FOR FIGURES AND DATA ANALYSIS

FUNCTIONS

## ===================================================================== 

## FUNCTIONS 

## MLfunctions.R 

## 

## >>> utility function for general use in the ML functions 

Expand <- function(p, Aux) { 

 unlist(lapply(Aux, function(z) { 

 switch(z[[1]], 

 Ested={p[z[[2]]]}, 

 Fixed=z[[2]], 

 Both={ 

 r1 <- z[[3]] 

 r1[z[[2]]] <- p[z[[3]][z[[2]]]] 

 r1 

 } 

 ) 

 })) 

} 

TEST.Expand <- FALSE 

if (TEST.Expand){ 

## Example, estimate: 

## Parameter Index 

## sigma 1 

## A 2,3,4 

## n 5,6,7 

## Fixed: 

## w value = -1 

## Estimate the first two components of K, with indices 8,9, 

## while fixing the third to -1.9 

## Want the order in the final vector to be sigma, A, w, n, K, so 

Aux <- list(sigma=list("Ested", c(1)), 

 A = list("Ested", c(2,3,4)), 

 w = list("Fixed", c(-1)), 

 n = list("Ested", c(5,6,7)), 

 K = list("Both", c(TRUE,TRUE,FALSE), c(8, 9, -1.9))) 

## Start with 

p <- c(0.1, .8, .9, 1.0, 1.1,1.2,1.3, 2.1,2.2) 

## End with 

Expand(p, Aux) 

## should be 

## 0.1, 0.8, 0.9, 1.0, -1, 1.1, 1.2, 1.3, 2.1, 2.2, -1.9 

} 

## >>> hill log-likelihood model for the conventional parametrization 

## >>> of the hill model. There are some tricks to allow for the 

## >>> possibility of block- or sex-specific values for parameters. 

## >>> Arguments are: 

## >>> p: the vector of parameters 
24 

## >>> Data: data frame within which to evaluate 

## >>> the loglikelihood. 

## >>> It includes sT4, the scaled value of T4, and sdose, 

## >>> the scaled value of dose (dose/max(dose)), for computational 

## >>> stability. 

## >>> It also includes vectors for each parameter 

## >>> indexing into p that gives the value of p to use for each 

## >>> observation in Data. These indices must be as long as 

## >>> dose and sT4, and each index falls in the range 1:k, where 

## >>> k is the maximum level for that parameter. factors converted 

## >>> to numeric are ideal for this. Next, Data has the attribute 

## >>> 'nElts', which is a vector with names 'sigma','A','w','n',and 

## >>> 'K', which gives the largest value for the corresponding index 

## >>> Finally, the argument Aux is either 

## >>> NULL (default) or a list with elements sigma, A, w, n, and K. 

## >>> The value of each element is a list with two components. 

## >>> The first is "Ested", "Fixed", or "Both". If "Ested", the 

## >>> second component gives the indices in the parameter vector 

## >>> p that corresponds to the desired values. If "Fixed", the 

## >>> second component gives the fixed values to use for that 

## >>> parameter in the estimation process. If "Both", the second 

## >>> component is a logical, TRUE if the corresponding level 

## >>> of the parameter is to be estimated, FALSE otherwise. For 

## >>> the "Both" case, there is a third component, as long as the 

## >>> second, which gives the index in p of the corresponding 

## >>> level of the parameter, or the value itself if FALSE. 

## >>> This is implemented in the function Expand(), defined above. 

## >>> 

hillloglikK <- 

 function(p, Data, model=NULL, Aux=NULL) { 

 if (is.null(model)) stop("model must be defined") 

 if (!is.null(Aux)) 

 p <- Expand(p, Aux) 

 nElts <- attr(Data, "nElts") 

 out <- with(Data,{ 

 y <- get(as.character(model[[2]])) 

 x <- get(as.character(model[[3]])) 

 sigma <- p[sigmaI] 

 A <- exp(p[nElts["sigma"] + AI]) 

 w <- p[nElts["sigma"] + nElts["A"] + wI] 

 n <- exp(p[nElts["sigma"] + nElts["A"] + nElts["w"] + nI]) 

 K <- exp(p[nElts["sigma"] + nElts["A"] + nElts["w"] + nElts["n"] + KI]) 

 ## predictions, given parameters. 

 mu <- 

 ifelse(x == 0, 

 A, 

 ifelse(x >= K, 

 A*(1 + w/((K/x)^n + 1)), 

 A*(1 + w*(x/K)^n/((x/K)^n + 1)))) 

 ## Deviance. 

 Dev <- -2*sum(dlnorm(y, meanlog=log(mu), sdlog=sigma, 

 log=TRUE)) 

 list(mu=mu, Dev=Dev) 

 }) 

 out$Dev 
25 

 } 

## >>> hillmodK is just the prediction part of the above. Just returns 

## >>> 'mu' 

hillmodK <- 

 function(x, p, Data, Aux=NULL) { 

 if (!is.null(Aux)) 

 p <- Expand(p, Aux) 

 nElts <- attr(Data, "nElts") 

 with(Data,{ 

 sigma <- p[sigmaI] 

 A <- exp(p[nElts["sigma"] + AI]) 

 w <- p[nElts["sigma"] + nElts["A"] + wI] 

 n <- exp(p[nElts["sigma"] + nElts["A"] + nElts["w"] + nI]) 

 K <- exp(p[nElts["sigma"] + nElts["A"] + nElts["w"] + nElts["n"] + KI]) 

 ## predictions, given parameters. 

 ifelse(x == 0, 

 A, 

 ifelse(x >= K, 

 A*(1 + w/((K/x)^n + 1)), 

 A*(1 + w*(x/K)^n/((x/K)^n + 1)))) 

 }) 

 } 

## >>> hill log-likelihood model for lognormal error model. 

## >>> This formulation is for the form of the model which 

## >>> has BMD for relative change of the mean as a parameter. 

## >>> That is, D is the dose where response is (1 - R) times 

## >>> background, and R is the benchmark response level. 

## >>> So, this is also for decreasing dose-response. 

## >>> This function is set up to be used as on objective function 

## >>> for nlminb or optim, and uses the dataframe Data as does 

## >>> hillloglikK. R must be set as an attribute of Data. 

## >>> arguments are like those of hillloglikK, replacing K with D, 

## >>> the benchmark dose for a 100*R% decrease from control. 

## >>> This formulation is mainly used for computing profile- 

## >>> likelihood confidence intervals for the BMD.

hillloglikD <- 

 function(p, Data, model=NULL, Aux=NULL, R=0.10) {

 if (is.null(model)) stop("model must be defined") 

 if (!is.null(Aux)) 

 p <- Expand(p, Aux) 

 nElts <- attr(Data, "nElts") 

 with(Data,{ 

 y <- get(as.character(model[[2]])) 

 x <- get(as.character(model[[3]])) 

 sigma <- p[sigmaI] 

 A <- exp(p[nElts["sigma"] + AI]) 

 w <- p[nElts["sigma"] + nElts["A"] + wI] 

 n <- exp(p[nElts["sigma"] + nElts["A"] + nElts["w"] + nI]) 

 D <- exp(p[nElts["sigma"] + nElts["A"] + nElts["w"] + nElts["n"] + KI]) 

 ## predictions, given parameters. 

 beta <- -(w/R + 1) 

 mu <- 

 ifelse(x == 0, 
26 

 A, 

 ifelse(x >= D, 

 A*(1 + w/(beta*(D/x)^n + 1)),

 A*(1 + w*(x/D)^n/((x/D)^n - (1 + w/R))))) 

 ## Deviance 

 -2*sum(dlnorm(y, meanlog=log(mu), sdlog=sigma, 

 log=TRUE)) 

 }) 

 } 

hillmodD <- function(x, p, Data, Aux=NULL, R=0.10) { 

 if (!is.null(Aux)) 

 p <- Expand(p, Aux) 

 nElts <- attr(Data, "nElts") 

 with(Data,{ 

 sigma <- p[sigmaI] 

 A <- exp(p[nElts["sigma"] + AI]) 

 w <- p[nElts["sigma"] + nElts["A"] + wI] 

 n <- exp(p[nElts["sigma"] + nElts["A"] + nElts["w"] + nI]) 

 D <- exp(p[nElts["sigma"] + nElts["A"] + nElts["w"] + nElts["n"] + KI]) 

 ## predictions, given parameters. 

 beta <- -(w/R + 1) 

 ifelse(x == 0, 

 A, 

 ifelse(x >= D, 

 A*(1 + w/(beta*(D/x)^n + 1)), 

 A*(1 + w*(x/D)^n/((x/D)^n - (1 + w/R))))) 

 }) 

} 

## Fixup: We are fitting models to scaled parameters, but want 

## to convert back to unscaled versions. This means adding 

## 2 * nrow(dta) * log(controlmean) to the 'objective' value of the 

## returned object from nlminb, adding log(controlmean) to all the 

## values of 'lA', and adding log(maxdose) to all the values of 

## lK. Their indices in par can be found from Aux[["A"]][[2]], for instance 

Fixup <- function(z, n, Aux) { 

 z$objective <- z$objective + 2 * n * log(controlmean) 

 z$par[Aux[["A"]][[2]]] <- z$par[Aux[["A"]][[2]]] + log(controlmean) 

 z$par[Aux[["K"]][[2]]] <- z$par[Aux[["K"]][[2]]] + log(maxdose) 

 z 

} 

## a little function to compute AIC for this analysis 

myAIC <- function(Dev, Df){ 

 Dev + 2*Df 

} 

## and BIC 

myBIC <- function(Dev, Df, n){ 

 Dev + Df * log(n) 

} 

## Compute BMDs from n and K, for Hill model: 

toBMD <- function(n, K, R=0.10) { 

 K * (R / (1 - R)) ^ (1/n) 

} 
27 

R SCRIPTS FOR ANALYSIS

Preliminaries 

## Can I combine the Male and Female data sets. So, we fit 

## a Full model, with separate variance terms for each block and background 

## as well as for each of the three scale and shape parameters for the 

## Hill model. Compare the likelihood to the Reduced model, 

## Which has variance and background terms like the full model, but 

## has only single shape and scale parameters. When I did this in the 

## Bayesian model, the full model had substantially lower DIC than did the 

## reduced model. However, a model with sex-specific shape and scale 

## parameters had slightly greater DIC than the reduced model, which does not 

## really make sense. 

## ----------------------------------------------------------------------- 

## Libraries used: 

## 

library(ggplot2) 

library(MASS) 

library(coda) 

library(ellipse) 

source("MLfunctions.R") 

## ----------------------------------------------------------------------- 

## Constants 

## 

## The golden ratio and its inverse 

## For setting up aspect ratios for plotting 

golden <- 1.61803399 

igolden <- 1/golden 

## ------------------------------------------------------------------------ 

## Data 

## 

## Preliminaries. Loading the data, evaluating a likely error 

## model, with lots of graphics. 

## ... Load the data, and add some scaled versions of the response, T4 

## ... and dose. 

load("CombinedPubertal.RData") 

## data in dta 

## scale dose and T4 in dta 

maxdose <- max(dta$dose) 

dta$sdose <- dta$dose/maxdose 

## Scale T4 so that the overall control mean is 1 

controlmean <- mean(dta$T4[dta$dose == 0]) 

dta$sT4 <- dta$T4 / controlmean 

## ----------------------------------------------------------------------- 

## Plot the design for these data: doses by block, and sample size 

## 

## Figure 1 

win.metafile("Doses.wmf", width=3.18*golden, height=3.18, pointsize=12) 

tmp <- aggregate(dta$T4, by=list(Block=dta$block, Dose=dta$dose), FUN=length) 

sfn <- function(x) c(rep(0,6),1,2,3,4,0,0,0,0,5)[x]

qplot(Dose, factor(Block), data=tmp, asp=igolden, size=sfn(x), 

 ylab="Block") + 
28 

 scale_size("Sample\nSize", breaks=c(1,2,3,4,5), 

 labels=c("7","8","9","10","15"), to=c(1,5)) 

dev.off() 

## ---------------------------------------------------------------------- 

## Distribution of the data, for the statistical model 

## 

## Figure 5 

## ........ Box Cox transformation for the data: 

win.metafile("BoxCox.wmf", width=3.18*golden, height=3.18, pointsize=12) 

par(bty="l", mar=c(4,4,0,0), las=1) 

boxcox(lm(T4 ~ factor(dose):factor(block) - 1, data=dta)) 

dev.off() 

## ------------------------------------------------------------------- 

## 'Saturated' Model 

## Fit a linear model to log(T4) ~ factor(dose):factor(block) 

## and get predicted values for each point. Also, extract 'sigma' 

## and convert it to the mle by squaring, multiplying by df/n, then 

## taking the square root. 

## ........ Do the three blocks share a commong CV? 

## ........ 1) Saturated model with block-specific CVs 

tmp1 <- lm(log(T4) ~ factor(dose), data=subset(dta, block==1)) 

tmp2 <- lm(log(T4) ~ factor(dose), data=subset(dta, block==2)) 

tmp3 <- lm(log(T4) ~ factor(dose), data=subset(dta, block==3)) 

pred <- c(predict(tmp1), predict(tmp2), predict(tmp3)) 

mlsigma <- 

 c(rep(sqrt((summary(tmp1)$sigma)^2 * 

tmp1$df.residual/length(tmp1$residuals)), 

 length(tmp1$residuals)), 

 rep(sqrt((summary(tmp2)$sigma)^2 * 

tmp2$df.residual/length(tmp2$residuals)), 

 length(tmp2$residuals)), 

 rep(sqrt((summary(tmp3)$sigma)^2 * 

tmp3$df.residual/length(tmp3$residuals)), 

 length(tmp3$residuals))) 

Dev0a <- -2*sum(dlnorm(dta$T4, pred, mlsigma, log=TRUE)) 

DF0a <- sum(!is.na(tmp1$coef)) + 1 + 

 sum(!is.na(tmp2$coef)) + 1 + 

 sum(!is.na(tmp3$coef)) + 1 ## model df. 

## ........ 2) Saturated model with single, common CV 

tmp <- lm(log(T4) ~ factor(dose):factor(block), data=dta) 

pred <- predict(tmp) 

mlsigma <- sqrt((summary(tmp)$sigma)^2 * 

tmp$df.residual/length(tmp$residuals)) 

Dev0 <- -2*sum(dlnorm(dta$T4, pred, mlsigma, log=TRUE)) 

DF0 <- sum(!is.na(tmp$coef)) + 1 ## model df. 

## Are the cvs heterogeneous? 

pchisq(Dev0 - Dev0a, df=DF0a - DF0, lower.tail=FALSE) 

## P is about 0.14. MLEs for sdlogs are: 

## Block sdlog 
29 

## 1 0.19 

## 2 0.22 

## 3 0.27 

## Use the homogeneous formulation. 

## ........ Plot the standard deviations and geometric standard deviations 

## ........ (10^sd(log10(T4))) of the data. 

## Figure 6 

win.metafile("standarddeviations.wmf", 

 width=3.18*1.61803399, height=3.1, pointsize=12) 

tmp <- aggregate(dta$T4, by=list(Block=dta$block, Dose=dta$dose), FUN=sd) 

tmp2 <- aggregate(log10(dta$T4), by=list(Block=dta$block, Dose=dta$dose), 

 FUN=function(z) 10^sd(z)) 

p <- ggplot(tmp, aes(Dose, x))+ 

 opts(aspect.ratio=1/1.61803399)+ 

 geom_point(aes(colour=factor(Block))) 

p + geom_point(aes(x=tmp2$Dose, y=tmp2$x, colour=factor(tmp2$Block), 

shape=1)) 

dev.off() 

## ........ Residual plots for log-transformed data

## Figure 7 

win.metafile("residslog10T4.wmf", 

 width=3.18, height=3.18, pointsize=12)

m1 <- lm(log10(T4) ~ factor(block) : factor(dose), data=dta) 

qplot(sample=resid(m1)/mad(resid(m1)), dist=qnorm)+stat_abline(intercept=0, 

slope=1) 

dev.off() 

## It looks to be a reasonable error model, that all the data are lognormal 

## with a common CV. 

## ---------------------------------------------------------------------- 

## Dose-Response Shape 

## 

## Work on the log scale. This has the advantage 

## that, if the response is exponential, we should see straight lines. 

## Plot the log10 transformed data 

dta$l10T4 <- log10(dta$T4) 

## Fit a linear trend for each block, and retain the residuals. 

tmp <- lm(l10T4 ~ factor(block) + factor(block):dose - 1, data=dta) 

tmp2 <- lm(l10T4 ~ 1 + factor(block):dose, data=dta) 

dta$res <- resid(tmp) 

## Figure 8 

win.metafile("lmeansandfits.wmf", width=3.18*golden, height=3.18, 

pointsize=12) 

ggplot(dta, aes(dose, l10T4, colour=factor((block)))) + 

 opts(aspect.ratio=igolden) + 

 scale_colour_hue(name="Block") + 

 ylab(expression(log[10](T4))) + 

 stat_summary(fun.data="mean_cl_normal", geom="pointrange") + 

 geom_abline(intercept=coef(tmp)[1], slope=coef(tmp)[4], colour=2) + 

 geom_abline(intercept=coef(tmp)[2], slope=coef(tmp)[5], colour=3) + 
30 

 geom_abline(intercept=coef(tmp)[3], slope=coef(tmp)[6], colour=4) 

dev.off() 

## Figure 9 

win.metafile("reslinearfit.wmf", width=3.18*golden, height=3.18, 

pointsize=12) 

ggplot(dta, aes(dose, res, colour=factor(block))) +

 opts(aspect.ratio=igolden) + 

 scale_colour_hue(name="Block") + 

 ylab("Residuals from linear fit") + 

 geom_hline(yintercept=0) + 

 stat_summary(fun.y="mean", geom="point") 

dev.off() 

## Test, do the controls differ? 

anova(lm(l10T4 ~ factor(block), subset(dta, dose == 0))) 

## In brief, yes. The P-value for this is 0.0461. 

Dose-Response Analysis 

## =================================================================== 

## Dose-Response Analysis 

## 

## ------------------------------------------------------------------ 

## Fit the n = 1 model to the data in dta. 

## The index for each of these is block. 

## 1) set up Aux 

Aux <- list(sigma=list("Ested", c(1)), 

 A = list("Ested", c(2,3,4)), 

 w = list("Ested",c(5,6,7)), 

 n = list("Fixed", 0), 

 K = list("Ested", c(8, 9, 10))) 

## Set up and populate the input with data and indices 

FitData <- cbind(dta, sigmaI=rep(1, nrow(dta)), 

 AI=dta$block, 

 wI=dta$block, 

 nI=rep(1, nrow(dta)), 

 KI=dta$block) 

attr(FitData, "nElts") <- c(sigma=1, A=3, w=3, n=1, K=3) 

## Initial values. Use the medians from the Bayesian analysis (there's a 

## switch!) 

## Full variance model: different CV for each block

start <- c(0.2, ## sdlog 

 log(c(1.084, 0.9301, 0.8338)), # lA

 c(-0.8984, -0.871, -0.8687), # w

 log(c(0.2791,0.138,0.1249)) # lK

 ) 

## To enforce the bounds on w, use nlminb 

## sigma A w K 

lb <- c(0.01,rep(-Inf, 3),rep(-1,3),rep(-Inf, 3)) 

ub <- c(3, rep(Inf, 3), rep(-0.1,3),rep(Inf, 3)) 

MLfitFull.n <- nlminb(start, hillloglikK, gradient=NULL, hessian=NULL, 

 Data=FitData, 
31 

 model=sT4 ~ sdose, 

 Aux=Aux, 

 control=list(eval.max=1000, iter.max=1000), 

 lower=lb, 

 upper=ub) 

MLfitFull.n <- Fixup(MLfitFull.n, nrow(FitData), Aux) 

## Compare this model to the full one. 

## The fit and resulting deviance is from data scaled by controlmean 

## Add 2 * nrow(dta) * log(controlmean) to get deviance on the original 

## scale 

Dev.n <- MLfitFull.n$objective 

DF.n <- length(MLfitFull.n$par) 

P.value.n <- pchisq(Dev.n - Dev0, DF0 - DF.n, lower.tail=FALSE) 

## about 0.02. not so good. 

## ---------------------------------------------------------------------- 

## Fit the w = -1 model to the data in dta 

## The index for each of these is block. 

## 1) set up Aux 

Aux <- list(sigma=list("Ested", c(1)), 

 A = list("Ested", c(2,3,4)), 

 w = list("Fixed",-1), 

 n = list("Ested",c(5,6,7)), 

 K = list("Ested", c(8, 9, 10))) 

## Set up and populate the input with data and indices 

FitData <- cbind(dta, sigmaI=rep(1, nrow(dta)), 

 AI=dta$block, 

 wI=rep(1, nrow(dta)), 

 nI=dta$block, 

 KI=dta$block) 

attr(FitData, "nElts") <- c(sigma=1, A=3, w=1, n=3, K=3) 

## Initial values. Use the medians from the Bayesian analysis (there's a 

## switch!) 

## Full variance model: different CV for each block

start <- c(0.2, ## sdlog 

 log(c(1.084, 0.9301, 0.8338)), # lA 

 c(0,0,0), # ln

 log(c(0.2791,0.138,0.1249)) # lK

 ) 

## To enforce the bounds on w, use nlminb 

## sigma A n K 

lb <- c(0.01,rep(-Inf, 3),rep(-18,3),rep(-Inf, 3)) 

ub <- c(3, rep(Inf, 3), rep(log(20),3),rep(Inf, 3)) 

MLfitFull.w <- nlminb(start, hillloglikK, gradient=NULL, hessian=NULL, 

 Data=FitData, 

 model=sT4 ~ sdose, 

 Aux=Aux, 

 control=list(eval.max=1000, iter.max=1000), 

 lower=lb, 

 upper=ub) 
32 

MLfitFull.w <- Fixup(MLfitFull.w, nrow(FitData), Aux) 

Dev.w <- MLfitFull.w$objective 

DF.w <- length(MLfitFull.w$par) 

P.value.w <- pchisq(Dev.w - Dev0, DF0 - DF.w, lower.tail=FALSE) 

## about 0.84; good fit. So, these curves go to 0, but with some 

## complications wrt shape. Check residuals and fit before proceeding. 

## ............ Plots 

pred <- 

 hillmodK(FitData$dose, MLfitFull.w$par, Data=FitData, Aux=Aux) 

res <- log(dta$T4) - log(pred) 

## Plot res vs pred, mean(res) vs dose for each block, 

## fitted curve vs data. 

tdta <- FitData 

tdta$res <- res 

tdta$pred <- pred 

## Plot mean residual +/- 1 sd 

win.metafile("wfixed_res_pred.wmf", 

 width=3.18*golden, height=3.18, pointsize=12) 

ggplot(tdta, aes(pred, res, colour=factor(block))) + 

 stat_summary(fun.data="mean_sdl", geom="pointrange", mult=1) + 

 opts(aspect.ratio=igolden) + 

 scale_colour_hue(name="Block") 

dev.off() 

win.metafile("wfixed_res_sdose.wmf", 

 width=3.18*golden, height=3.18, pointsize=12) 

ggplot(tdta, aes(dose, res, colour=factor(block))) + 

 stat_summary(fun.data="mean_sdl", geom="pointrange", mult=1) + 

 opts(aspect.ratio=igolden) + 

 scale_colour_hue(name="Block") 

dev.off() 

## A wrapper around hillmodK to plot the fitted curve from this 

## fit for a specific block. This function is specific to the 

## w = -1 fit, MLfitFull.w 

phillmodK <- function(x, block){ 

 ## Need to point AI, nI, and KI to block 

 tData <- 

 data.frame(AI = rep(block, length(x)), 

 sigmaI = rep(1, length(x)), 

 wI = rep(1, length(x)), 

 nI = rep(block, length(x)), 

 KI = rep(block, length(x))) 

 p <- MLfitFull.w$par 

 attr(tData, "nElts") <- attr(FitData, "nElts") 

 hillmodK(x, p, tData, Aux) 

} 

## figure 10 

win.metafile("wfixed_DRfun.wmf", 
33 

 width=3.18*golden, height=3.18, pointsize=12) 

ggplot(tdta, aes(dose, T4, colour=factor(block))) +

 opts(aspect.ratio=igolden) + 

 scale_colour_hue(name="Block") + 

 stat_summary(fun.data="mean_cl_normal", geom="pointrange") + 

 stat_function(fun=phillmodK, args=list(block=1), colour=2) + 

 stat_function(fun=phillmodK, args=list(block=2), colour=3) + 

 stat_function(fun=phillmodK, args=list(block=3), colour=4) 

dev.off() 

## Get contours for the posterior distribution of n and K from the 

## Bayesian analysis of these data, and plot them by block. The point is 

## the similarity in shape and scale between blocks 2 and 3, and the 

## distinctness of block 1. Plot the MLEs on the same axes as plotting 

## points. The samples are in the subdirectory 'Bayes'. 

full <- 

 mcmc.list(read.coda(file.path(".","Bayes","Model A1","Full_chain1.txt"), 

 file.path(".","Bayes","Model A1","Full_index.txt"), 

 quiet=TRUE), 

 read.coda(file.path(".","Bayes","Model A1","Full_chain2.txt"), 

 file.path(".","Bayes","Model A1","Full_index.txt"), 

 quiet=TRUE), 

 read.coda(file.path(".","Bayes","Model A1","Full_chain3.txt"), 

 file.path(".","Bayes","Model A1","Full_index.txt"), 

 quiet=TRUE)) 

samps <- 

 as.data.frame(rbind(cbind(do.call("rbind", full[,c("K[1]","n[1]")]), 

 block=rep(1,3*nrow(full[[1]]))), 

 cbind(do.call("rbind", full[,c("K[2]","n[2]")]), 

 block=rep(2,3*nrow(full[[1]]))), 

 cbind(do.call("rbind", full[,c("K[3]","n[3]")]), 

 block=rep(3,3*nrow(full[[1]]))) 

 ) 

 ) 

names(samps) <- c("K","n","block") 

samps[,"block"] <- factor(samps[,"block"]) 

## Rescale K by max(dose) == 300 

samps[,"K"] <- 300*samps[,"K"] 

pp <- exp(MLfitFull.w$par[5:10]) 

## Use ellipse.default to plot an approximate 95% confidence ellipse 

s1 <- subset(samps, block==1) 

cov1 <- with(s1, cov(cbind(K,n))) 

e1 <- ellipse(cov1, centre=pp[c(4,1)], npoints=301)

s2 <- subset(samps, block==2) 

cov2 <- with(s2, cov(cbind(K,n))) 

e2 <- ellipse(cov2, centre=pp[c(5,2)], npoints=301)

s3 <- subset(samps, block==3) 

cov3 <- with(s3, cov(cbind(K,n))) 

e3 <- ellipse(cov3, centre=pp[c(6,3)], npoints=301)

## Figure 11. Uses Model A1, in JAGS. 

win.metafile("PlotnvsK.wmf", width=4.75, height=3.18, pointsize=12) 
34 

ggplot(samps, aes(K,n, colour=block)) + geom_density2d() + 

 geom_path(aes(x=e1[,"K"], y=e1[,"n"]), col=2, size=1) + 

 geom_path(aes(x=e2[,"K"], y=e2[,"n"]), col=3, size=1) + 

 geom_path(aes(x=e3[,"K"], y=e3[,"n"]), col=4, size=1) + 

 scale_y_continuous(limits=c(0.2, 2.5)) + 

 scale_x_continuous(limits=c(10, 120)) + 

 geom_point(aes(x=pp[4], y=pp[1]), col=2) + 

 geom_point(aes(x=pp[5], y=pp[2]), col=3) + 

 geom_point(aes(x=pp[6], y=pp[3]), col=4) 

dev.off() 

## --------------------------------------------------------------------- 

## Can we fit combined shape and scale parameters to the 

## two male blocks? 

## The index for each of these is block. 

## 1) set up Aux 

Aux <- list(sigma=list("Ested", c(1)), 

 A = list("Ested", c(2,3,4)), 

 w = list("Fixed",-1), 

 n = list("Ested",c(5,6)), 

 K = list("Ested", c(7, 8))) 

## Set up and populate the input with data and indices 

FitData <- cbind(dta, sigmaI=rep(1, nrow(dta)), 

 AI=dta$block, 

 wI=rep(1, nrow(dta)), 

 nI=as.numeric(dta$sex), 

 KI=as.numeric(dta$sex)) 

attr(FitData, "nElts") <- c(sigma=1, A=3, w=1, n=2, K=2) 

## Initial values. Use the medians from the Bayesian analysis (there's a 

## switch!) 

## Full variance model: different CV for each block

## Use parameters from the last fit to get starting values 

par <- MLfitFull.w$par 

start <- c(par[1], ## sdlog

 par[2:4], # lA

 mean(par[5:6]),par[7], # ln

 mean(par[8:9]),par[10]) # lK

## To enforce the bounds on w, use nlminb 

## sigma A n K 

lb <- c(0.01,rep(-Inf, 3),rep(-18,2),rep(-Inf, 2)) 

ub <- c(3, rep(Inf, 3), rep(log(20),2),rep(Inf, 2)) 

MLfitSex <- nlminb(start, hillloglikK, gradient=NULL, hessian=NULL, 

 Data=FitData, 

 model=sT4 ~ sdose, 

 Aux=Aux, 

 control=list(eval.max=1000, iter.max=1000), 

 lower=lb, 

 upper=ub) 

MLfitSex <- Fixup(MLfitSex, nrow(FitData), Aux) 
35 

Dev.Sex <- MLfitSex$objective 

DF.Sex <- length(MLfitSex$par) 

P.value.Sex <- pchisq(Dev.Sex - Dev.w, DF.w - DF.Sex, lower.tail=FALSE) 

## This is quite a bit worse than the fit to all three blocks: 

## P ~ 0.0017 

## See what a plot looks like 

## ................... Plot 

## A wrapper around hillmodK to plot the fitted curve from this 

## fit for a specific block. This function is specific to the 

## fit with sex-specific shape and scale parameters, MLfitSex 

phillmodK <- function(x, block){ 

 sex <- c(2,2,1)[block] 

 ## Need to point AI, nI, and KI to block 

 tData <- 

 data.frame(AI = rep(block, length(x)), 

 sigmaI = rep(1, length(x)), 

 wI = rep(1, length(x)), 

 nI = rep(sex, length(x)), 

 KI = rep(sex, length(x))) 

 p <- MLfitSex$par 

 attr(tData, "nElts") <- attr(FitData, "nElts") 

 hillmodK(x, p, tData, Aux) 

} 

ggplot(FitData, aes(dose, T4, colour=factor(block))) + 

 opts(aspect.ratio=igolden) + 

 scale_colour_hue(name="Block") + 

 stat_summary(fun.data="mean_cl_normal", geom="pointrange") + 

 stat_function(fun=phillmodK, args=list(block=1), colour=2) + 

 stat_function(fun=phillmodK, args=list(block=2), colour=3) + 

 stat_function(fun=phillmodK, args=list(block=3), colour=4) 

## -------------------------------------------------------------------- 

## Let scale (K) change among blocks, but shape be sex-specific 

## The index for each of these is block. 

## 1) set up Aux 

Aux <- list(sigma=list("Ested", c(1)), 

 A = list("Ested", c(2,3,4)), 

 w = list("Fixed",-1), 

 n = list("Ested",c(5,6)), 

 K = list("Ested", c(7, 8, 9))) 

## Set up and populate the input with data and indices 

FitData <- cbind(dta, sigmaI=rep(1, nrow(dta)), 

 AI=dta$block, 

 wI=rep(1, nrow(dta)), 

 nI=as.numeric(dta$sex), 

 KI=as.numeric(dta$block)) 

attr(FitData, "nElts") <- c(sigma=1, A=3, w=1, n=2, K=3) 

## Use parameters from the last fit to get starting values 

par <- MLfitFull.w$par 

start <- c(par[1], ## sdlog 

 par[2:4], # lA
36 

 mean(par[5:6]),par[7], # ln

 par[8:10]) # lK

## To enforce the bounds on w, use nlminb 

## sigma A n K 

lb <- c(0.01,rep(-Inf, 3),rep(-18,2),rep(-Inf, 3)) 

ub <- c(3, rep(Inf, 3), rep(log(20),2),rep(Inf, 3)) 

## From now on, just number the fits 

MLfit.04 <- nlminb(start, hillloglikK, gradient=NULL, hessian=NULL, 

 Data=FitData, 

 model=sT4 ~ sdose, 

 Aux=Aux, 

 control=list(eval.max=1000, iter.max=1000), 

 lower=lb, 

 upper=ub) 

MLfit.04 <- Fixup(MLfit.04, nrow(FitData), Aux) 

Dev.04 <- MLfit.04$objective 

DF.04 <- length(MLfit.04$par) 

P.value.04 <- pchisq(Dev.04 - Dev.w, DF.w - DF.04, lower.tail=FALSE) 

## P-value is 0.0019 

## Blocks 1 and 2 have different shapes and scales.

## ................. Plot 

## A wrapper around hillmodK to plot the fitted curve from this 

## fit for a specific block. 

phillmodK <- function(x, block){ 

 ## Need to point AI, nI, and KI to block 

 sex <- c(2,2,1)[block] 

 tData <- 

 data.frame(AI = rep(block, length(x)), 

 sigmaI = rep(1, length(x)), 

 wI = rep(1, length(x)), 

 nI = rep(sex, length(x)), 

 KI = rep(block, length(x))) 

 p <- MLfit.04$par 

 attr(tData, "nElts") <- attr(FitData, "nElts") 

 hillmodK(x, p, tData, Aux) 

} 

ggplot(FitData, aes(dose, T4, colour=factor(block))) + 

 opts(aspect.ratio=igolden) + 

 scale_colour_hue(name="Block") + 

 stat_summary(fun.data="mean_cl_normal", geom="pointrange") + 

 stat_function(fun=phillmodK, args=list(block=1), colour=2) + 

 stat_function(fun=phillmodK, args=list(block=2), colour=3) + 

 stat_function(fun=phillmodK, args=list(block=3), colour=4) 

## ----------------------------------------------------------------------- 

## Are blocks 2 and 3 different? If not, maybe something odd happened the 

## first time out (block 1). 

## Let scale (K) change among blocks, but shape be sex-specific 

## The index for each of these is block. 

## 1) set up Aux 

Aux <- list(sigma=list("Ested", c(1)), 
37 

 A = list("Ested", c(2,3,4)), 

 w = list("Fixed",-1), 

 n = list("Ested",c(5,6)), 

 K = list("Ested", c(7, 8))) 

## Set up and populate the input with data and indices 

FitData <- cbind(dta, cblock=c(1,2,2)[dta$block]) 

FitData <- cbind(FitData, sigmaI=rep(1, nrow(dta)),

 AI=FitData$block, 

 wI=rep(1, nrow(dta)), 

 nI=as.numeric(FitData$cblock), 

 KI=as.numeric(FitData$cblock)) 

attr(FitData, "nElts") <- c(sigma=1, A=3, w=1, n=2, K=2) 

## Use parameters from the last fit to get starting values 

par <- MLfitFull.w$par 

start <- c(par[1], ## sdlog 

 par[2:4] - log(controlmean), # lA 

 par[5],mean(par[6:7]), # ln 

 c(par[8],mean(par[9:10])) - log(maxdose)) # lK 

## To enforce the bounds on w, use nlminb 

## sigma A n K 

lb <- c(0.01,rep(-Inf, 3),rep(-18,2),rep(-Inf, 2)) 

ub <- c(3, rep(Inf, 3), rep(log(20),2),rep(Inf, 2)) 

MLfit.05 <- nlminb(start, hillloglikK, gradient=NULL, hessian=NULL, 

 Data=FitData, 

 model=sT4 ~ sdose, 

 Aux=Aux, 

 control=list(eval.max=1000, iter.max=1000), 

 lower=lb, 

 upper=ub) 

MLfit.05 <- Fixup(MLfit.05, nrow(FitData), Aux) 

Dev.05 <- MLfit.05$objective 

DF.05 <- length(MLfit.05$par) 

P.value.05 <- pchisq(Dev.05 - Dev.w, DF.w - DF.05, lower.tail=FALSE) 

## P-value is 0.72. So, Blocks 2 and 3 are not very different; it makes 

## sense 

## to combine them, but Block 1 stands out as different. Again, plots: 

## ........ Plot 

## A wrapper around hillmodK to plot the fitted curve from this 

## fit for a specific block. 

phillmodK <- function(x, block){ 

 ## Need to point AI, nI, and KI to block 

 sex <- c(2,2,1)[block] 

 cblock <- c(1,2,2)[block] 

 tData <- 

 data.frame(AI = rep(block, length(x)), 

 sigmaI = rep(1, length(x)), 

 wI = rep(1, length(x)), 

 nI = rep(cblock, length(x)), 
38 

 KI = rep(cblock, length(x))) 

 p <- MLfit.05$par 

 attr(tData, "nElts") <- attr(FitData, "nElts") 

 hillmodK(x, p, tData, Aux) 

} 

## What if block 1 had the same shape and scale as blocks 2 and 3? 

phillmodK2 <- function(x){ 

 sex <- 2 

 cblock <- 2 

 block <- 1 

 tData <- 

 data.frame(AI = rep(block, length(x)), 

 sigmaI = rep(1, length(x)), 

 wI = rep(1, length(x)), 

 nI = rep(cblock, length(x)), 

 KI = rep(cblock, length(x))) 

 p <- MLfit.05$par 

 attr(tData, "nElts") <- attr(FitData, "nElts") 

 hillmodK(x, p, tData, Aux) 

} 

ggplot(FitData, aes(dose, T4, colour=factor(block))) + 

 opts(aspect.ratio=igolden) + 

 scale_colour_hue(name="Block") + 

 stat_summary(fun.data="mean_cl_normal", geom="pointrange") + 

 stat_function(fun=phillmodK, args=list(block=1), colour=2) + 

 stat_function(fun=phillmodK, args=list(block=2), colour=3) + 

 stat_function(fun=phillmodK, args=list(block=3), colour=4) + 

 stat_function(fun=phillmodK2, colour=2, linetype=2) 

## --------------------------------------------------------------------- 

## Try fitting common shape and scale parameters 

## 1) set up Aux 

Aux <- list(sigma=list("Ested", c(1)), 

 A = list("Ested", c(2,3,4)), 

 w = list("Fixed",-1), 

 n = list("Ested",c(5)), 

 K = list("Ested", c(6))) 

## Set up and populate the input with data and indices 

FitData <- cbind(dta, sigmaI=rep(1, nrow(dta)), 

 AI=FitData$block, 

 wI=rep(1, nrow(dta)), 

 nI=rep(1, nrow(dta)), 

 KI=rep(1, nrow(dta))) 

attr(FitData, "nElts") <- c(sigma=1, A=3, w=1, n=1, K=1) 

## Use parameters from the last fit to get starting values 

par <- MLfitFull.w$par 

start <- c(par[1]*4, ## sdlog 

 par[2:4] - log(controlmean), # lA

 mean(par[5:7]), # ln

 mean(par[8:10]) - log(maxdose)) # lK

## To enforce the bounds on w, use nlminb 
39 

## sigma A n K 

lb <- c(0.01,rep(-Inf, 3),rep(-18,1),rep(-Inf, 1)) 

ub <- c(3, rep(Inf, 3), rep(log(20),1),rep(Inf, 1)) 

MLfit.06 <- nlminb(start, hillloglikK, gradient=NULL, hessian=NULL, 

 Data=FitData, 

 model=sT4 ~ sdose, 

 Aux=Aux, 

 control=list(eval.max=1000, iter.max=1000), 

 lower=lb, 

 upper=ub) 

MLfit.06 <- Fixup(MLfit.06, nrow(FitData), Aux) 

Dev.06 <- MLfit.06$objective 

DF.06 <- length(MLfit.06$par) 

P.value.06 <- pchisq(Dev.06 - Dev.w, DF.w - DF.06, lower.tail=FALSE) 

## Not a good fit: P.value = 0.0065 

Benchmark Doses from Blocks 2 and 3 

## ====================================================================== 

## Benchmark Dose Computations 

## and Confidence Intervals 

## ------ First method; drop block 1, just fit blocks 2 and 3, combined 

## Set up table for BMDs and CLs. 

BMDsave1 <- matrix(0, nrow=2, ncol=6, dimnames=list(c("0.10","0.20"), 

 c("BMD","BMDL","BMDU","BMD.Bayes","BMDL.Bayes", 

 "BMDU.Bayes"))) 

## ...... R = 0.10 

## make rdta, the data with only blocks 2 and 3 

rdta <- subset(dta, block %in% c(2,3)) 

## First, fit the full model, get the BMD estimate for the restricted dataset 

Aux <- list(sigma=list("Ested", c(1)), 

 A = list("Ested", c(2,3)), 

 w = list("Fixed",-1), 

 n = list("Ested",c(4)), 

 K = list("Ested", c(5))) 

## Set up and populate the input with data and indices 

FitData <- rdta 

FitData <- cbind(FitData, sigmaI=rep(1, nrow(rdta)), 

 AI=(FitData$block - 1), 

 wI=rep(1, nrow(rdta)), 

 nI=rep(1, nrow(rdta)), 

 KI=rep(1, nrow(rdta))) 

attr(FitData, "nElts") <- c(sigma=1, A=2, w=1, n=1, K=1) 

## Use parameters from the last fit to get starting values 

## This seems to be robust, use lK for lD initial value 

par <- MLfit.05$par 

start <- c(par[1], ## sdlog

 par[3:4] - log(controlmean), # lA
40 

 mean(par[5:6]), # ln

 mean(par[7:8]) - log(maxdose)) # lK

## To enforce the bounds on w, use nlminb 

## sigma A n D 

lb <- c(0.01,rep(-Inf, 2),rep(-18,1),rep(-Inf, 1)) 

ub <- c(3, rep(Inf, 2), rep(log(20),1),rep(Inf, 1)) 

MLfit.t <- nlminb(start, hillloglikD, gradient=NULL, hessian=NULL, 

 Data=FitData, 

 model=sT4 ~ sdose, 

 Aux=Aux,R=0.10, 

 control=list(eval.max=1000, iter.max=1000), 

 lower=lb, 

 upper=ub) 

## Save the BMD 

BMDsave1[1,"BMD"] <- exp(MLfit.t$par[5]) * maxdose 

## ----- BMD Confidence Intervals 

## Use parameters from the last fit to get starting values 

par <- MLfit.t$par 

start <- c(par[1], ## sdlog

 par[2:3], # lA

 par[4]) # ln

Aux <- list(sigma=list("Ested", c(1)), 

 A = list("Ested", c(2,3)), 

 w = list("Fixed",-1), 

 n = list("Ested",c(4)), 

 K = list("Fixed", par[5])) 

## To enforce the bounds on w, use nlminb 

## sigma A n 

lb <- c(0.01, rep(-Inf, 2),rep(-18,1)) 

ub <- c(3, rep(Inf, 2), rep(log(20),1)) 

## Now, iterate to lower values for D, for the BMDL. 

## Decrement D by 1% each time (add log(0.99) to log(D)) 

## The target value for objective is: 

Dev <- MLfit.t$objective 

TX2 <- Dev + qchisq(0.90, 1) 

niter <- 0 

BMDl <- list(par[5]) 

Deviancel <- list(Dev) 

while (niter < 100 & Dev <= TX2) { 

 Aux[["K"]][[2]][1] <- Aux[["K"]][[2]][1] + log(0.99) 

 MLfit.u <- 

 try(nlminb(start, hillloglikD, gradient=NULL, hessian=NULL, 

 Data=FitData, 

 model=sT4 ~ sdose, 

 Aux=Aux,R=0.10, 

 control=list(eval.max=1000, iter.max=1000), 
41 

 lower=lb, 

 upper=ub)) 

 if (!is(MLfit.u, "try-error") && MLfit.u$convergence == 0) { 

 start <- MLfit.u$par 

 Dev <- MLfit.u$objective 

 BMDl <- c(BMDl, list(Aux[["K"]][[2]][1])) 

 Deviancel <- c(Deviancel, list(Dev)) 

 } 

} 

## Now, turn this sequence of BMDs and Deviances into a spline function 

sf <- splinefun(unlist(Deviancel), unlist(BMDl)) 

BMDsave1[1,"BMDL"] <- exp(sf(TX2))*maxdose 

## Now the upper bound. 

Aux <- list(sigma=list("Ested", c(1)), 

 A = list("Ested", c(2,3)), 

 w = list("Fixed",-1), 

 n = list("Ested",c(4)), 

 K = list("Fixed", par[5])) 

Dev <- MLfit.t$objective 

niter <- 0 

start <- par[1:4] 

BMDu <- list(par[5]) 

Devianceu <- list(Dev) 

while (niter < 100 & Dev <= TX2) { 

 Aux[["K"]][[2]][1] <- Aux[["K"]][[2]][1] + log(1.01) 

 MLfit.u <- 

 try(nlminb(start, hillloglikD, gradient=NULL, hessian=NULL, 

 Data=FitData, 

 model=sT4 ~ sdose, 

 Aux=Aux,R=0.10, 

 control=list(eval.max=1000, iter.max=1000), 

 lower=lb, 

 upper=ub)) 

 if (!is(MLfit.u, "try-error") && MLfit.u$convergence == 0) { 

 start <- MLfit.u$par 

 Dev <- MLfit.u$objective 

 BMDu <- c(BMDu, list(Aux[["K"]][[2]][1])) 

 Devianceu <- c(Devianceu, list(Dev)) 

 } 

} 

## Now, turn this sequence of BMDs and Deviances into a spline function 

sf <- splinefun(unlist(Devianceu), unlist(BMDu)) 

BMDsave1[1,"BMDU"] <- exp(sf(TX2))*maxdose 

## Combine them into one dataframe. 

TESTING <- FALSE 

if (TESTING) { 

lldta <- data.frame(BMD=exp(c(rev(unlist(BMDl)),unlist(BMDu)) + 

log(maxdose)), 

 Deviance=c(rev(unlist(Deviancel)),unlist(Devianceu))) 
42 

ggplot(lldta, aes(BMD, Deviance)) + geom_line() + geom_hline(yintercept=TX2) 

} 

xx <- BMDsave1[1,"BMDL"] 

DRfita1 <- DRfita1 + 

 geom_segment(aes(x=xx, xend=xx, y=0, yend=0.90), color="green") 

DRfita1 

## ...... R = 0.20 

## First, fit the full model, get the BMD estimate for the restricted dataset 

Aux <- list(sigma=list("Ested", c(1)), 

 A = list("Ested", c(2,3)), 

 w = list("Fixed",-1), 

 n = list("Ested",c(4)), 

 K = list("Ested", c(5))) 

## Set up and populate the input with data and indices 

FitData <- rdta 

FitData <- cbind(FitData, sigmaI=rep(1, nrow(rdta)), 

 AI=(FitData$block - 1), 

 wI=rep(1, nrow(rdta)), 

 nI=rep(1, nrow(rdta)), 

 KI=rep(1, nrow(rdta))) 

attr(FitData, "nElts") <- c(sigma=1, A=2, w=1, n=1, K=1) 

## Use parameters from the last fit to get starting values 

## This seems to be robust, use lK for lD initial value 

par <- MLfit.05$par 

start <- c(par[1], ## sdlog

 par[3:4] - log(controlmean), # lA 

 mean(par[5:6]), # ln 

 mean(par[7:8]) - log(maxdose)) # lD

## To enforce the bounds on w, use nlminb 

## sigma A n D 

lb <- c(0.01,rep(-Inf, 2),rep(-18,1),rep(-Inf, 1)) 

ub <- c(3, rep(Inf, 2), rep(log(20),1),rep(Inf, 1)) 

MLfit.t <- nlminb(start, hillloglikD, gradient=NULL, hessian=NULL, 

 Data=FitData, 

 model=sT4 ~ sdose, 

 Aux=Aux,R=0.20, 

 control=list(eval.max=1000, iter.max=1000), 

 lower=lb, 

 upper=ub) 

## Save the BMD 

BMDsave1[2,"BMD"] <- exp(MLfit.t$par[5]) * maxdose 

## ----- BMD Confidence Intervals 

## Use parameters from the last fit to get starting values 
43 

par <- MLfit.t$par 

start <- c(par[1], ## sdlog 

 par[2:3], # lA 

 par[4]) # ln 

Aux <- list(sigma=list("Ested", c(1)), 

 A = list("Ested", c(2,3)), 

 w = list("Fixed",-1), 

 n = list("Ested",c(4)), 

 K = list("Fixed", par[5])) 

## To enforce the bounds on w, use nlminb 

## sigma A n 

lb <- c(0.01, rep(-Inf, 2),rep(-18,1)) 

ub <- c(3, rep(Inf, 2), rep(log(20),1)) 

## Now, iterate to lower values for D, for the BMDL. 

## Decrement D by 1% each time (add log(0.99) to log(D)) 

## The target value for objective is: 

Dev <- MLfit.t$objective 

TX2 <- Dev + qchisq(0.90, 1) 

niter <- 0 

BMDl <- list(par[5]) 

Deviancel <- list(Dev) 

while (niter < 100 & Dev <= TX2) { 

 Aux[["K"]][[2]][1] <- Aux[["K"]][[2]][1] + log(0.99) 

 MLfit.u <- 

 try(nlminb(start, hillloglikD, gradient=NULL, hessian=NULL, 

 Data=FitData, 

 model=sT4 ~ sdose, 

 Aux=Aux,R=0.20, 

 control=list(eval.max=1000, iter.max=1000), 

 lower=lb, 

 upper=ub)) 

 if (!is(MLfit.u, "try-error") && MLfit.u$convergence == 0) { 

 start <- MLfit.u$par 

 Dev <- MLfit.u$objective 

 BMDl <- c(BMDl, list(Aux[["K"]][[2]][1])) 

 Deviancel <- c(Deviancel, list(Dev)) 

 } 

} 

## Now, turn this sequence of BMDs and Deviances into a spline function 

sf <- splinefun(unlist(Deviancel), unlist(BMDl)) 

BMDsave1[2,"BMDL"] <- exp(sf(TX2))*maxdose 

## Now the upper bound. 

Aux <- list(sigma=list("Ested", c(1)), 

 A = list("Ested", c(2,3)), 

 w = list("Fixed",-1), 

 n = list("Ested",c(4)), 

 K = list("Fixed", par[5])) 
44 

Dev <- MLfit.t$objective 

niter <- 0 

start <- par[1:4] 

BMDu <- list(par[5]) 

Devianceu <- list(Dev) 

while (niter < 100 & Dev <= TX2) { 

 Aux[["K"]][[2]][1] <- Aux[["K"]][[2]][1] + log(1.01) 

 MLfit.u <- 

 try(nlminb(start, hillloglikD, gradient=NULL, hessian=NULL, 

 Data=FitData, 

 model=sT4 ~ sdose, 

 Aux=Aux,R=0.20, 

 control=list(eval.max=1000, iter.max=1000), 

 lower=lb, 

 upper=ub)) 

 if (!is(MLfit.u, "try-error") && MLfit.u$convergence == 0) { 

 start <- MLfit.u$par 

 Dev <- MLfit.u$objective 

 BMDu <- c(BMDu, list(Aux[["K"]][[2]][1])) 

 Devianceu <- c(Devianceu, list(Dev)) 

 } 

} 

## Now, turn this sequence of BMDs and Deviances into a spline function 

sf <- splinefun(unlist(Devianceu), unlist(BMDu)) 

BMDsave1[2,"BMDU"] <- exp(sf(TX2))*maxdose 

## Combine them into one dataframe. 

TESTING <- TRUE 

if (TESTING) { 

lldta <- data.frame(BMD=exp(c(rev(unlist(BMDl)),unlist(BMDu)) + 

log(maxdose)), 

 Deviance=c(rev(unlist(Deviancel)),unlist(Devianceu))) 

ggplot(lldta, aes(BMD, Deviance)) + geom_line() + geom_hline(yintercept=TX2) 

} 

## Plot data and curve. Plot everything in terms of measurements 

## relative to the corresponding control value. Start with just 

## a linear model fit of log-transformed data, then transform estimates 

## and confidence intervals back to the original scale. 

t.lm <- lm(l10T4 ~ -1 + factor(block) + factor(dose), data=rdta) 

## Build a matrix of medians and confidence intervals 

## The first two rows are block means; drop them. Computations 

## were based on log10 values 

t.relests <- 10^cbind(coef(t.lm)[-(1:2)], confint(t.lm)[-(1:2),]) 

colnames(t.relests) <- c("Est","lower","upper") 

## Add the blocks back, and dose as a numeric variable 

t.relests <- as.data.frame(cbind(t.relests, 

 block=c(2,3,2,3,3,3,2), 

 dose=as.numeric(sub("^factor\\(dose\\)","", 

 rownames(t.relests))))) 

t.relests <- rbind(data.frame(Est=c(1,1), lower=c(1,1), upper=c(1,1), 
45 

 block=c(2,3), dose=c(0,0)), 

 t.relests) 

## A wrapper around hillmodD to plot the fitted curve from this 

## fit. This function is specific to the 

## w = -1 fit, MLfit.t, R=0.20 

phillmodD <- function(x){ 

 x <- x/maxdose 

 ## Need to point AI, nI, and KI to block 

 tData <- 

 data.frame(AI = rep(1, length(x)), 

 sigmaI = rep(1, length(x)), 

 wI = rep(1, length(x)), 

 nI = rep(1, length(x)), 

 KI = rep(1, length(x))) 

 p <- MLfit.t$par 

 p[2:3] <- 0 

 ## This model is for relative change from control, so 

 ## I modify the parameter vector so that lA is 0, and 

 ## nElts is 1 for every parameter 

 nElts <- attr(FitData, "nElts") 

 nElts[2] <- 1 

 attr(tData, "nElts") <- attr(FitData, "nElts") 

 ## x is on the original scale. internally, though, x is [0,1] 

 hillmodD(x, p, tData, Aux, R=0.20) 

} 

xx <- BMDsave1[1,"BMD"] 

xxl <- BMDsave1[1,"BMDL"] 

xx2 <- BMDsave1[2,"BMD"] 

xx2l <- BMDsave1[2,"BMDL"] 

DRfita1 <- ggplot(t.relests, aes(x=dose, y=Est, colour=factor(block))) + 

 scale_colour_hue(name="Block") + 

 geom_segment(aes(x=xx, xend=xx, 

 y=0,yend=0.90), colour="green") + 

 geom_segment(aes(x=0, xend=xx,y=0.90,yend=0.90), color="green") + 

 geom_segment(aes(x=xxl, xend=xxl, y=0, yend=0.90), color="green") + 

 geom_segment(aes(x=xx2, xend=xx2, y=0,yend=0.80), color="blue") + 

 geom_segment(aes(x=0, xend=xx2, y=0.80, yend=0.80), color="blue") + 

 geom_segment(aes(x=xx2l, xend=xx2l, y=0, yend=0.80), color="blue") + 

 stat_function(fun=phillmodD, colour=1) + 

 geom_pointrange(aes(ymin=lower, ymax=upper, data=t.relests)) + 

 scale_y_continuous(name="T4 Relative to Block Control") + 

 scale_x_continuous(name="Dose") 

## Figure 3 

png(file="A2BMDplot.png", width=3.18*golden, height=3.18, pointsize=12, 

 units="in", res=600) 

print(DRfita1) 

dev.off() 

## The Bayesian estimates come from sampling carried out by JAGS. We 

## just dump the data set here, run the sampling in scripts in the folder 

## 'Model A2', 
46 

## then read in the samples and compute BMDs and confidence limits. 

## The Bayes method has the advantage that we can compute both BMD[LU]10 

## and BMD[LU]20 from the same set of samples, of A, K, and n. 

ldta <- as.list(rdta[,c("sdose","sT4","block")]) 

ldta$N <- nrow(rdta) 

ldta$maxdose <- maxdose 

ldta$controlmean <- controlmean 

dput(ldta, file="rT4data.txt") 

## This needs to be edited to remove some attributes. 

## ... Now run JAGS 

## ... all done 

## Read in the chains 

A2 <- 

 mcmc.list(read.coda(file.path(".","Bayes","Model A2", 

 "Pooled_chain1.txt"), 

 file.path(".","Bayes","Model A2", 

 "Pooled_index.txt"), 

 quiet=TRUE), 

 read.coda(file.path(".","Bayes","Model A2", 

 "Pooled_chain2.txt"), 

 file.path(".","Bayes","Model A2", 

 "Pooled_index.txt"), 

 quiet=TRUE), 

 read.coda(file.path(".","Bayes","Model A2", 

 "Pooled_chain3.txt"), 

 file.path(".","Bayes","Model A2", 

 "Pooled_index.txt"), 

 quiet=TRUE)) 

## Quick test for convergence, using gelman.diag. 

print(gelman.diag(A2)) 

## converged 

## Turn A2 into a matrix 

A2x <- do.call("rbind", lapply(A2, unclass)) 

BMD10 <- A2x[,"BMD10"] 

BMDsave1[1,"BMD.Bayes"] <- mean(BMD10) 

BMDsave1[1,c("BMDL.Bayes", "BMDU.Bayes")] <- 

 quantile(BMD10, prob=c(0.05, 0.95)) 

BMD20 <- A2x[,"BMD20"] 

BMDsave1[2,"BMD.Bayes"] <- mean(BMD20) 

BMDsave1[2,c("BMDL.Bayes", "BMDU.Bayes")] <- 

 quantile(BMD20, prob=c(0.05, 0.95)) 

Hierarchical Analysis 

## ======================================================================= 

## Approach 2: Use a multilevel model in JAGS 

A3.j2 <- mcmc.list(read.coda(file.path(".","Bayes","Model A3", 

 "s2_chain1.txt"), 

 file.path(".","Bayes","Model A3", 
47 

 "s2_index.txt"), 

 thin=10,quiet=TRUE), 

 read.coda(file.path(".","Bayes","Model A3","s2_chain2.txt"), 

 file.path(".","Bayes","Model A3","s2_index.txt"), 

 thin=10,quiet=TRUE), 

 read.coda(file.path(".","Bayes","Model A3","s2_chain3.txt"), 

 file.path(".","Bayes","Model A3","s2_index.txt"), 

 thin=10,quiet=TRUE)) 

vn <- varnames(A3.j2) 

vn2test <- vn[c(1:14,16)] ## The last two are functions of the

## remaining parameters,and lnK.var[2,1] == lnK.var[1,2], which 

## messes up the gelman diagnostic. 

heidel.diag(A3.j2[, vn2test]) 

gelman.diag(A3.j2[, vn2test]) 

A3x.j2 <- do.call("rbind", lapply(A3.j2, unclass)) 

## We compute BMDs respecting the hierarchical structure. So, sample from 

## the uncertainty distribution of lnK.mu[1] == log(n) and 

## lnK.mu[2] == log(K). 

## Then compute the expected BMD given those population parameters by 

## sampling repeatedly (so that the maximum se of the mean parameters 

## is < 0.01) 

## from the population implied, and computing the 

## mean BMD10 and 20. 

popparms <- A3x.j2[, (11:16)] 

BMDdist <- 

 apply(popparms, 1, 

 function(x) { 

 mu <- x[1:2] 

 sigma <- matrix(x[3:6], nrow=2) 

 ss <- ceiling(max(diag(sigma)/0.0001)) 

 lprms <- mvrnorm(ss, mu, sigma) 

 c(mean(toBMD(exp(lprms[,1]), exp(lprms[,2]), R=0.10)), 

 mean(toBMD(exp(lprms[,1]), exp(lprms[,2]), R=0.20))) * maxdose 

 } 

 ) 

 

BMDsave3.j2 <- rbind("0.10"=c(mean(BMDdist[1,]),median(BMDdist[1,]), 

 quantile(BMDdist[1,], prob=c(0.05, 0.95))), 

 "0.20"=c(mean(BMDdist[2,]),median(BMDdist[2,]), 

 quantile(BMDdist[2,], prob=c(0.05, 0.95)))) 

colnames(BMDsave3.j2) <- 

c("mean(BMD).Bayes","median(BMD).Bayes","BMDL.Bayes","BMDU.Bayes") 

BMDsave3.j2 

## Plot data (relative to concurrent background estimate) with CIs, 

## along with block-specific dr curve, population dr curve, and indicators 

## for BMDs and BMDLs 

## First, means etc. for log response. This works because the only dose 

## level that appears in more than one block is 0. 
48 

t.lm2 <- lm(l10T4 ~ -1 + factor(block) + factor(dose), data=dta) 

## Build a matrix of medians and confidence intervals 

## The first two rows are block means; drop them. Computations 

## were based on log10 values 

t.relests2 <- 10^cbind(coef(t.lm2)[-(1:3)], confint(t.lm2)[-(1:3),]) 

colnames(t.relests2) <- c("Est","lower","upper") 

## Add the blocks back, and dose as a numeric variable 

t.relests2 <- as.data.frame(cbind(t.relests2, 

 block=c(2,3,2,3,3,1,3,1,2), 

 dose=as.numeric(sub("^factor\\(dose\\)","", 

 rownames(t.relests2))))) 

t.relests2 <- rbind(data.frame(Est=c(1,1,1), lower=c(1,1,1), upper=c(1,1,1), 

 block=c(1,2,3), dose=c(0,0,0)), 

 t.relests2) 

## Function to fit the model to each parameter value in A3x.j2, over the 

## values in x, then return the means. block 0 is the top-level. 

## x will range over the real range of dose, but estimates are for 

## scaled dose, 

## and for log10 dose 

phillmodK <- function(x, block){ 

 ## x <- (10^x)/maxdose 

 x <- x/maxdose 

 ## cols in Pars are n, K 

 Pars <- A3x.j2[,c(4+block, 7+block)] 

 w <- -1 

 y <- apply(Pars, 1, function(p) { 

 A <- 1 

 n <- p[1] 

 K <- p[2] 

 ifelse(x == 0, 

 A, 

 ifelse(x >= K, 

 A*(1 + w/((K/x)^n + 1)), 

 A*(1 + w*(x/K)^n/((x/K)^n + 1))))

 }) 

 apply(y, 1, mean) 

} 

phillmodK0 <- function(x, pctile){ 

 ## x <- (10^x)/maxdose 

 x <- x/maxdose 

 popparms <- A3x.j2[sample(1:nrow(A3x.j2), 1000, replace=FALSE),(11:16)] 

 w <- -1 

 A <- 1 

 ymx <- 

 apply(popparms, 1, 

 function(z) { 

 mu <- z[1:2] 

 sigma <- matrix(z[3:6], nrow=2) 

 ss <- min(3000,ceiling(max(diag(sigma)/0.0001))) 

 Pars <- exp(mvrnorm(ss, mu, sigma)) 

 y <- apply(Pars, 1, function(p) { 

 n <- p[1] 

 K <- p[2] 

 ifelse(x == 0, 
49 

 A, 

 ifelse(x >= K, 

 A*(1 + w/((K/x)^n + 1)), 

 A*(1 + w*(x/K)^n/((x/K)^n + 1)))) 

 }) 

 apply(y, 1, mean) 

 }) 

 apply(ymx,1, 

 function(z){ 

 quantile(z, prob=pctile/100) 

 }) 

} 

DRfita2 <- ggplot(t.relests2, aes(x=dose, y=Est, colour=factor(block))) + 

 scale_colour_hue(name="Block") + 

 geom_pointrange(aes(ymin=lower, ymax=upper, data=t.relests2)) + 

 scale_x_continuous(name="Dose") + 

 scale_y_continuous(name="T4 Relative to Block Control") 

## Figure 2 

png(file="RelChangefromControl.png", width=3.18*golden, height=3.18, 

 pointsize=12, units="in", res=600) 

print(DRfita2) 

dev.off() 

DRfita2 <- DRfita2 + 

 stat_function(fun=phillmodK, colour=2, args=list(block=1), size=0.25) + 

 stat_function(fun=phillmodK, colour=3, args=list(block=2), size=0.25) + 

 stat_function(fun=phillmodK, colour=4, args=list(block=3), size=0.25) 

## use phillmodK0 to get the median, 5th and 95th percentiles of the 

## destribution of effects 

xx <- seq(0, maxdose, length=101) 

y <- phillmodK0(xx, c(5,50,95)) 

DRfita20 <- 

 DRfita2 + 

 geom_line(aes(x=x, y=y), data=data.frame(x=xx, y=y[2,]), colour=1, size=1) 

+ 

 geom_line(aes(x=x, y=y), data=data.frame(x=xx, y=y[1,]), 

 colour=1, size=0.25, linetype=2) + 

 geom_line(aes(x=x, y=y), data=data.frame(x=xx, y=y[3,]), 

 colour=1, size=0.25, linetype=2) 

 

## Figure 4 

png(file="A3BMDplot.png", width=3.18*golden, height=3.18, pointsize=12, 

 units="in", res=600) 

print(DRfita20) 

dev.off()