### 1 

# D OSE 

# -R ESPONSE 

# A NALYSIS OF 

# C OMBINED 

# P RE



# P UBERTAL 

# D ATA

J OINT 

## M

ALE AND 

## F

EMALE 

## B ENCHMARK 

## D OSES



# C ONTENTS

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 co mmon 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



# T HE 

# P ROBLEM

There are three blocks of data from the pre-puberta l assays. Two are male-only blocks, and one is female-only. The main goal is to estimate benchmar k doses for 10% and 20% reduction in T4 activity for these data. If the data support it, t his should be a common set of benchmark doses for males and females, otherwise separate sets, each wi th their own set of profile-likelihood-based

confidence limits.

This begs the question, how to determine whether th e male-only and female-only data sets support a common dose-response. The factor of interest, ge nder, is almost perfectly confounded with blocking. Thus, the magnitude and nature of differ ences among blocks affect the interpretation of any observed difference between the genders. In pr inciple, this should boil down to asking whether the difference between male- and female-specific be nchmark 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-on ly blocks has only two non-control dose groups, so we cannot fit the preferred hill model (four par ameters) to it:

### 3

Figure 1 Doses and sample sizes in the pre-puberta l dataset. Blocks 1 and 2 are male-only, block 3 i s female only.



### So, Block 2 (males) has both the smallest and large st non-control doses, while Block 1 (males) has no low doses at all. As a consequence, Block 1 dat a 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 f or 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 nu mber of dose levels in block 1 means that either relatively simple models, with no more than three p arameters, must fit the data, or that block must be left out of comparisons that require more comple x 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 wo uld simplify the analysis if the same model shape param eters 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 t o determine why, and evaluate how that affects the comparison with the single female block. All c ode for this analysis, carried out in R, is include d in the Appendices.



# S UMMARY OF 

# R ESULTS

The Hill model fits each block quite well. A logno rmal multiplicative error model describes these data adequately, and the coefficient of variation i s the same across all three blocks. The term in t he Hill model that controls the asymptotic value of th e response as dose increases was fixed, so that the level of T4 would go to 0 as dose increases. T his was done so that block 1 could be fit by the same dose-response model used for the other two blo cks. 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 differen t 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 puber tal assay blocks, relative to block-specific contro ls. 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 assump tions is supported by the data).



### I considered two interpretations of this result:

1) “ ‘Something’ went wrong with block 1, causing it t o 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 block s are so similar, combine their data, and use only them to estimate BMDs.” 2) “We have no direct information about what, if anyth ing, might have gone wrong with block 1. The variation among the three blocks may just r eflect normal or typical variation among experimental blocks, due to minor differences in hu sbandry, 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 “populatio n” of such studies, and base it on a statistically valid estimate of the average dose-re sponse.”

In either case, since gender seems to make no diffe rence 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 y ields benchmark doses of:

### 5

Table 1 BMDs and 90% confidence limit (thus, the B MDL corresponds to the conventional 95% one-sided l imit used for the BMDL) for the analysis of blocks 2 and 3 only, using maximum likelihood methods. Units f or 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 contain s only males, while block 3, only females), and fit ted Hill model. Data are expressed here as T4 levels relati ve to block-specific control (as in Figure 2).



### Under the second interpretation, the appropriate co urse 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. Eve n so, I was unable to reliably estimate model parameters, though, without adding additional infor mation 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 populati on 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 par ameter, K would mostly fall between 5% and 95% of the population median, and the population me dian most likely falls between 2% and 98% of the maximum dose. The mean dose-response curves und er these assumptions are shown in Figure 4. The distribution of BMDs is quite skewed (the m ean BMD 10 

### is nearly 200000, while the median is 10.37), so the median is used here as the measure o f central tendency. The central estimates of the BMDs are larger than for just blocks 2 and 3, as ar e the BMDLs, while the BMDUs are quite a bit

greater.

### 6

Table 2 BMDs and confidence limits from the Bayesi an analysis of the hierarchical model using all thr ee 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 do se-response curves from a Bayesian analysis. Color ed 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 distributio n of the population mean response.



### The specifics of the informative prior change the a nswer somewhat; less restrictive priors increase the mean estimate, while dropping BMDLs to lower le vels and increasing the BMDU.



## R ECOMMENDATION

It is difficult to arrive at a definitive recommend ation between these two analyses. After discussion with the investigators who generated the data, it a ppears 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 experimen ts themselves. On the other hand, that block has no doses with responses close to the benchmark resp onse, 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 th at 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 eit her method are not so different makes this

decision somewhat easier to make.

### 7



# D ETAILS OF THE 

# S TATISTICAL 

# A NALYSIS

S TATISTICAL 

## M ODEL FOR THE 

## D ATA 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, initi ally 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 ide ntify a power transformation to homogenize the

variances.

Figure 5 Box-cox estimated power for making varian ces homogeneous. This suggests a log-transformatio n (corresponding to lambda = 0).



### We can see this by plotting both dose- and block-sp ecific standard deviations and geometric standard deviations against dose. The standard dev iations 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 ge ometric standard deviations (10 sd(log10(T4)) ) (open symbols) by dose and block. This result, constancy of the g eometric standard deviation, is consistent with log -normality of
the data.



### Since this study was conducted in blocks, one likel y possibility is that, while coefficients of variat ion 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 geo metric standard deviation are similar among blocks. This question can be formally tested by co mparing the log-likelihood of a model fitting a separate mean to each dose X block group and a sepa rate 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 dat a should be well-described by a model with a

common coefficient of variation.

Finally, the quantile-quantile plot of residuals f rom dose- and block-specific means, on the log scale. These should form a straight line if log tr ansformation 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 log 10 transformed data, centered on doseXblock means and scaled by a robust estimate of its standard deviation, versus theoreti cal quantiles of the normal distribution. The stra ight line is consistent with a lognormal distribution of the ori ginal data.



### All this says that it is reasonable to assume the d ata are lognormally distributed about dose- and block-specific medians, with a constant coefficient of variation. That is 

# ( )
10 , 

### log ( 4) ~ , i j 

### T 

### μ σ N

The subscripts i

and j index dose and block, respectively.



## D OSE 

## -R ESPONSE 

## S HAPE It would certainly be convenient if log 10 

### (T4) were linear in dose. That would imply a simpl e exponential model for the dose-response, which coul d be estimated even for Block 1. However this is unlikely: 

### theoreticalsample

-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) dro p is more gradual, than is consistent with a straight li ne. 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 decreasi ng data, this can be accomplished in both the hill model and the exponential model by having a po wer 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 le vels probably differ among blocks (the P-value for the comparison is 0.0461), with the block 3 con trol being about 91% of the block 1 control

value.

Now, the modeling switches to maximum likelihood fi tting, using a lognormal error model. Since there are often only small differences between gene ral exponential model and hill model fits, and there is already a history of fitting these data wi th the hill model, the hill model is used to estima te BMDs for these data. However, initially I will tur n it into a three-parameter model in two different ways: either fix the hill coefficient n to 1 and continue to estimate a horizontal asympto te, 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 likel ihood compares favorably to the likelihood from a saturated model, in which each group is characteriz ed by a separate mean, whose deviance (that is,

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

The best-fitting model with n fixed to 1 has a deviance of 202.2, and 10 model d egrees of freedom, yielding a p-value of 0.018, that is, its fit to th e data is significantly worse than the fully satura ted 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 (th is 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, wit h 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 ar e 95% confidence contours for the estimates, based on a bivariate normal approximation to the posterior den sity.



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

common benchmark dose could be estimated from a cur ve fit to the data with common shape and scale parameters. However, such a model has a sign ificantly worse fit to the data than the reference model (P = 0.0065). Blocks 1 and 2 differ signifi cantly (P = 0.0019), but blocks 2 and 3 do not (P = 0.72), consistent with Figure 8. Given all this, t here 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 blo ck 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 s hape and scale parameter for each block, and a top-level overall estimate. The top-level es timate corresponds to the pooled estimate. This is a statistically supportable approach to est imating a pooled BMD for the whole
dataset, but would be less valid if it were determi ned 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 m ethods 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. 

### Kn

0.5

1.0 1.5 2.0 2.5

20 40 60 80 100 120 block

1

2

3

### 13



## A PPROACH 

## 1: E STIMATE 

## BMD FOR 

## P OOLED 

## B LOCKS 

## 2 AND 

## 3 ONLY Both the maximum likelihood method applied in the p revious 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 a nalysis, since that is the methodology used in approach 2.



## A PPROACH 

## 2: H IERARCHICAL 

## M ODEL FOR 

## A LL 

## T HREE 

## B LOCKS The code for the hierarchical approach is duplicate d in Appendix 3.iii.

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

that the results of the three blocks in hand are sa mples from a population of possible blocks, whose dose-response curves all differ from each other (du e 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 param eters, 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 th e logarithms of the Hill coefficient, n and scale

parameter K

, are sampled from a multivariate normal distribut ion with population mean ( lnk.mu ) and population variance matrix ( lnk.var ). Together these two comprise the population para meters. The presence of this second level is the reason for calling this a hierarchical approach.

Preliminary work with this model indicated that the re were not enough blocks to estimate the population parameters without additional restrictio ns. Consequently, priors were set on the population parameters to reflect plausible bounds o n the distribution of the Hill model parameters. It is important to remember that subsequent BMD est imates and, especially, confidence limits are sensitive to these priors. First of all, the logar ithm (only natural logarithms are used in this analysis) of the Hill coefficient was assigned a un iform prior on the interval (log(.1), log(10)). Th at is, an interval centered at 1, with the possibility of ten-fold variation on either side. This is der ived from the practical observation that Hill coefficien ts 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 unifor m on the interval (log(0.02), log(0.98)). This bound s the scale parameter (also the dose value at which the response is half its maximum value), esse ntially to fall within the range of the administered doses. Priors on the variance paramet ers 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, thinn ing 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 (P lummer et al, 2010), which indicated that all three chains had passed the default diagnostic crit eria, and Gelman and Rubin’s convergence diagnostic, also as implemented in coda, which give s the multivariate potential scale reduction

factor for the full set of chains as 1.00.

BMDs and BMDLs were computed respecting the hierarc hical 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 populatio n covariance matrix are sampled from the estimated po sterior 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.0 1. The mean BMD 10 

### and BMD 20 

### 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 BMD 10 

### is 215865, while the median is just 10.4 (with 5 th 

### and 95 th 

### percentiles of 3.5 and 35.3). The resulting bench mark doses and confidence limits are reported on page 6.



# R EFERENCES

Martyn Plummer (2010). JAGS Version 2.1.0 User Man ual. 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. ht tp://CRAN.R-project.org/package=coda

R Development Core Team (2010). R: A language and e nvironment for statistical computing. R Foundation for Statistical Computing, Vienna, Austr ia. ISBN 3-900051-07-0, URL http://www.Rproject.org.

### 15



# A PPENDICES

1. F IGURES 

## S UPPORTING 

## M ODELING 

## D ECISIONS

RESIDUAL PLOTS FOR 

## ML F IT 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 valu e. This is for T4 scaled so the grand mean of the control valu es is 1.

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

### predres

-0.4

-0.2 0.0 0.2 0.4

0.2 0.4 0.6 0.8 1.0 Block

1

2

3



### doseres

-0.4 -0.2 0.0 0.2 0.4

0 50 100 150 200 250 300 Block

1

2

3

### 16



## 2. D ATA 

## P REPROCESSING Code for preparing the data frames used in the subs equent 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 li nes. ## This time, it is the 'fdose' variable that is mi ssing for ## the summary rows; this is a "". dtaF <- subset(dtaF, !(fdose %in% c(""))) dtaF$dose <- as.numeric(gsub("Control","0.0", as.ch aracter(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. C ODE FOR 

## B AYESIAN 

## A NALYSES Both OpenBUGS version 3.1.1 rev 524, and JAGS versi on 2.1.0 were used for these computations. The code reported here is for JAGS, to make it easi er to report repeatable results.



## i. F

IT 

## S EPARATE 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[b lock[i]]) + 1))

+ (1 - g1[i]) * (1 + W[block[i]]*pow(ddose[i]/K[block[i]], n[b lock[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.029 18706776658, 1.46508792318391, 1.21339709272419, 1.1473832612547 4, 1.18280020575740, 1.14298233915678, 0.99984758711349, 0.9686219969898 45, 1.24650879231839, 1.03002533864238, 0.98978833660386, 0.8757834974947 13, 0.92733715635657, 0.537960334546286, 0.346834574863305, 0.47509001886 1095, 0.560174512755053, 0.43401474594677, 0.445540970489055, 0.680047247994 818, 0.467545580978872, 0.244565528015394, 0.261750080969346, 0.17037855550 6868, 0.145020861513841, 0.132446798376803, 0.280820743393854, 0.27474327954 4285, 0.237440225571072, 0.99942845167559, 0.820248051972794, 1.100440092209 80, 1.36868677246661, 1.06292747051763, 0.837223037207796, 0.675017622740 003, 0.75947341347711, 0.987063956257502, 0.820248051972794, 0.79384251938 5014, 1.16058602754863, 1.08451294556955, 1.16813046543085, 0.7326487454514 28, 0.798033873764027, 0.614242984244318, 0.703099697079388, 1.13921012021 566, 0.788184190973347, 0.500447712854122, 0.468802987292576, 0.66621577854 4076, 0.681723789746423, 0.475299586580045, 0.728457391072415, 0.57673036255 2154, 0.741450589647355, 0.586999180780735, 0.176246451637486, 0.16849244603 6312, 0.176246451637486, 0.181066509173351, 0.224656594715083, 0.18923965021 2425, 0.268037112537865, 0.29926270266151, 0.165348930252053, 0.110442187886 986, 0.538169902265237, 0.770999638019394, 0.827373354417116, 0.70980586408 5808, 1.51370763398045,

### 18

1.26453161614815, 0.820876755129646, 0.778753643620 568, 0.534607251043076, 1.02688182285812, 0.469012555011526, 1.023109603917 01, 0.432547771914115, 0.711691973556364, 0.485358837089676, 0.59873497304 1971, 0.404884833012631, 0.514279182304864, 0.570024195545733, 0.36276172150 3553, 0.642325058583703, 0.594753186381909, 0.316866391053364, 0.60145935338 8329, 0.373449675170036, 0.508830421612147, 0.327973480157747, 0.32126731315 1327, 0.419345005620225, 0.210615557545390, 0.379527139019604, 0.29926270266 151, 0.265941435348359, 0.211453828421193, 0.23681152241422, 0.209777286669 588, 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.33333333333 3333, 0.333333333333333, 0.333333333333333, 0.333333333333333, 0.66666666666 6667, 0.666666666666667, 0.666666666666667, 0.666666666666667, 0.66666666666 6667, 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.0 1, 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.0312 5, 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. F

IT JUST 

## B LOCKS 

## 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.0 1, 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.0312 5, 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.368686772466 61, 1.06292747051763, 0.837223037207796, 0.675017622740003, 0.75947341347 711, 0.987063956257502, 0.820248051972794, 0.793842519385014, 1.16058602754 863, 1.08451294556955, 1.16813046543085, 0.732648745451428, 0.798033873764 027, 0.614242984244318, 0.703099697079388, 1.13921012021566, 0.788184190973 347, 0.500447712854122, 0.468802987292576, 0.666215778544076, 0.68172378974 6423, 0.475299586580045, 0.728457391072415, 0.576730362552154, 0.74145058964 7355, 0.586999180780735, 0.176246451637486, 0.168492446036312, 0.17624645163 7486, 0.181066509173351, 0.224656594715083, 0.189239650212425, 0.26803711253 7865, 0.29926270266151, 0.165348930252053, 0.110442187886986, 0.53816990226 5237, 0.770999638019394, 0.827373354417116, 0.709805864085808, 1.51370763398 045, 1.26453161614815, 0.820876755129646, 0.778753643620568, 0.53460725104 3076, 1.02688182285812, 0.469012555011526, 1.02310960391701, 0.432547771914 115, 0.7116919733556364, 0.485358837089676, 0.598734973041971, 0.40488483301 2631, 0.514279182304864, 0.570024195545733, 0.362761721503553, 0.64232505858 3703, 0.594753186381909, 0.316866391053364, 0.601459353388329, 0.37344967517 0036, 0.508830421612147, 0.327973480157747, 0.321267313151327, 0.41934500562 0225, 0.210615557545390, 0.379527139019604, 0.29926270266151, 0.265941435348 359, 0.211453828421193, 0.23681152241422, 0.209777286669588, 0.343691059079 045, 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. M ULTILEVEL 

## M ODEL FOR 

## S HAPE AND 

## S CALE FOR 

## A LL 

## T HREE 

## B LOCKS



## , I NFORMATIVE



## P RIORS 

### 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 t he variances and ## population means gave surprisingly broad poste rior distributions. ## So, see what happens if we use a more informat ive 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 w ould fall between ## say, .33 and 3 for most (95%) of the experimen ts. This gives ## a value of .314 for a guess for the variance o f 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.04 28. 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.0 428

## The population mean of log n is unlikely to fa ll 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 f all 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[b lock[i]]) + 1))

+ (1 - g1[i]) * (1 + W[block[i]]*pow(ddose[i]/K[block[i]], n[b lock[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.5 1), .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 C ODE FOR 

## F IGURES AND 

## D ATA 

## A NALYSIS



## F UNCTIONS

## ================================================ ===================== ## FUNCTIONS ## MLfunctions.R
## ## >>> utility function for general use in the ML f unctions 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 ind ices 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 convention al 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(do se)), for computational ## >>> stability. ## >>> It also includes vectors for each par ameter ## >>> indexing into p that gives the value of p to use for each ## >>> observation in Data. These indices m ust be as long as ## >>> dose and sT4, and each index falls in the range 1:k, where ## >>> k is the maximum level for that param eter. 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 fo r the corresponding index ## >>> Finally, the argument Aux is either ## >>> NULL (default) or a list with element s sigma, A, w, n, and K. ## >>> The value of each element is a list w ith two components. ## >>> The first is "Ested", "Fixed", or "Bo th". If "Ested", the ## >>> second component gives the indices in the parameter vector ## >>> p that corresponds to the desired val ues. If "Fixed", the ## >>> second component gives the fixed valu es to use for that ## >>> parameter in the estimation process. If "Both", the second ## >>> component is a logical, TRUE if the c orresponding level ## >>> of the parameter is to be estimated, FALSE otherwise. For ## >>> the "Both" case, there is a third com ponent, 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 E xpand(), 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"] + nElt s["w"] + nI]) K <- exp(p[nElts["sigma"] + nElts["A"] + nElt s["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), sdlo g=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"] + nElt s["w"] + nI]) K <- exp(p[nElts["sigma"] + nElts["A"] + nElt s["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 erro r model. ## >>> This formulation is for the form of the mode l 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-respons e. ## >>> This function is set up to be used as on obj ective 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, rep lacing K with D, ## >>> the benchmark dose for a 100*R% decrease fro m control. ## >>> This formulation is mainly used for computin g 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"] + nElt s["w"] + nI]) D <- exp(p[nElts["sigma"] + nElts["A"] + nElt s["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 parameter s, but want ## to convert back to unscaled versions. This mean s adding ## 2 * nrow(dta) * log(controlmean) to the 'objecti ve' value of the ## returned object from nlminb, adding log(controlm ean) to all the ## values of 'lA', and adding log(maxdose) to all t he 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(controlm ean) 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 analys is 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 S CRIPTS FOR 

## A NALYSIS


### Preliminaries

## Can I combine the Male and Female data sets. So , we fit ## a Full model, with separate variance terms for e ach block and background ## as well as for each of the three scale and shape parameters for the ## Hill model. Compare the likelihood to the Reduc ed model, ## Which has variance and background terms like the full model, but ## has only single shape and scale parameters. Whe n I did this in the ## Bayesian model, the full model had substantially lower DIC than did the ## reduced model. However, a model with sex-specif ic shape and scale ## parameters had slightly greater DIC than the red uced 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, evaluatin g 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 blo ck, 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, D ose=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, s ize=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 sta tistical model
## ## Figure 5 ## ........ Box Cox transformation for the data: win.metafile("BoxCox.wmf", width=3.18*golden, heigh t=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):fac tor(block) ## and get predicted values for each point. Also, extract 'sigma' ## and convert it to the mle by squaring, multiplyi ng 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(tmp 3)) 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=T RUE)) 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), dat a=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=TR UE)) DF0 <- sum(!is.na(tmp$coef)) + 1 ## model df.

## Are the cvs heterogeneous? pchisq(Dev0 - Dev0a, df=DF0a - DF0, lower.tail=FALS E)

## 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 geomet ric standard deviations ## ........ (10^sd(log10(T4))) of the data. ## Figure 6 win.metafile("standarddeviations.wmf", width=3.18*1.61803399, height=3.1, poi ntsize=12) tmp <- aggregate(dta$T4, by=list(Block=dta$block, D ose=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=fa ctor(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 al l 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 th e residuals. tmp <- lm(l10T4 ~ factor(block) + factor(block):dos e - 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="poi ntrange") + geom_abline(intercept=coef(tmp)[1], slope=coef(tm p)[4], colour=2) + geom_abline(intercept=coef(tmp)[2], slope=coef(tm p)[5], colour=3) +

### 30

geom_abline(intercept=coef(tmp)[3], slope=coef(tm p)[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 dt a. ## 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 indi ces 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 Bayesi an 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, ite r.max=1000), lower=lb, upper=ub) MLfitFull.n <- Fixup(MLfitFull.n, nrow(FitData), Au x) ## Compare this model to the full one. ## The fit and resulting deviance is from data scal ed by controlmean ## Add 2 * nrow(dta) * log(controlmean) to get devi ance 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 dt a ## 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 indi ces 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 Bayesi an 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, ite r.max=1000), lower=lb, upper=ub)

### 32

MLfitFull.w <- Fixup(MLfitFull.w, nrow(FitData), Au x) 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 fi t before proceeding. ## ............ Plots

pred <- hillmodK(FitData$dose, MLfitFull.w$par, Data=FitD ata, Aux=Aux) res <- log(dta$T4) - log(pred) ## Plot res vs pred, mean(res) vs dose for each blo ck, ## 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, points ize=12) ggplot(tdta, aes(pred, res, colour=factor(block))) + stat_summary(fun.data="mean_sdl", geom="poin trange", 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, points ize=12) ggplot(tdta, aes(dose, res, colour=factor(block))) + stat_summary(fun.data="mean_sdl", geom="poin trange", mult=1) + opts(aspect.ratio=igolden) + scale_colour_hue(name="Block") dev.off()

## A wrapper around hillmodK to plot the fitted cur ve from this ## fit for a specific block. This function is spec ific 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, points ize=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="poi ntrange") + 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 b y 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 s ame axes as plotting ## points. The samples are in the subdirectory 'Ba yes'.

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% c onfidence 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.1 8, pointsize=12)

### 34

ggplot(samps, aes(K,n, colour=block)) + geom_densit y2d() + geom_path(aes(x=e1[,"K"], y=e1[,"n"]), col=2, siz e=1) + geom_path(aes(x=e2[,"K"], y=e2[,"n"]), col=3, siz e=1) + geom_path(aes(x=e3[,"K"], y=e3[,"n"]), col=4, siz e=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 param eters 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 indi ces 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 Bayesi an 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=NUL L, hessian=NULL, Data=FitData, model=sT4 ~ sdose, Aux=Aux, control=list(eval.max=1000, ite r.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.Se x, lower.tail=FALSE)

## This is quite a bit worse than the fit to all th ree blocks: ## P ~ 0.0017 ## See what a plot looks like ## ................... Plot ## A wrapper around hillmodK to plot the fitted cur ve from this ## fit for a specific block. This function is spec ific 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="poi ntrange") + 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 indi ces 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=NUL L, hessian=NULL, Data=FitData, model=sT4 ~ sdose, Aux=Aux, control=list(eval.max=1000, ite r.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 cur ve 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="poi ntrange") + 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 som ething 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 indi ces 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=NUL L, hessian=NULL, Data=FitData, model=sT4 ~ sdose, Aux=Aux, control=list(eval.max=1000, ite r.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 ver y different; it makes
## sense ## to combine them, but Block 1 stands out as diffe rent. Again, plots: ## ........ Plot ## A wrapper around hillmodK to plot the fitted cur ve 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="poi ntrange") + 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 indi ces 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=NUL L, hessian=NULL, Data=FitData, model=sT4 ~ sdose, Aux=Aux, control=list(eval.max=1000, ite r.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 bloc ks 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 indi ces 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 v alue 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, ite r.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 lo g(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, h essian=NULL, Data=FitData, model=sT4 ~ sdose, Aux=Aux,R=0.10, control=list(eval.max=1000, ite r.max=1000),

### 41

lower=lb, upper=ub)) if (!is(MLfit.u, "try-error") && MLfit.u$converge nce == 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 in to 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, h essian=NULL, Data=FitData, model=sT4 ~ sdose, Aux=Aux,R=0.10, control=list(eval.max=1000, ite r.max=1000), lower=lb, upper=ub)) if (!is(MLfit.u, "try-error") && MLfit.u$converge nce == 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 in to 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)),unl ist(BMDu)) + log(maxdose)), Deviance=c(rev(unlist(Deviancel )),unlist(Devianceu)))

### 42

ggplot(lldta, aes(BMD, Deviance)) + geom_line() + g eom_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 indi ces 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 v alue 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, ite r.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 lo g(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, h essian=NULL, Data=FitData, model=sT4 ~ sdose, Aux=Aux,R=0.20, control=list(eval.max=1000, ite r.max=1000), lower=lb, upper=ub)) if (!is(MLfit.u, "try-error") && MLfit.u$converge nce == 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 in to 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, h essian=NULL, Data=FitData, model=sT4 ~ sdose, Aux=Aux,R=0.20, control=list(eval.max=1000, ite r.max=1000), lower=lb, upper=ub)) if (!is(MLfit.u, "try-error") && MLfit.u$converge nce == 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 in to 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)),unl ist(BMDu)) + log(maxdose)), Deviance=c(rev(unlist(Deviancel )),unlist(Devianceu)))

ggplot(lldta, aes(BMD, Deviance)) + geom_line() + g eom_hline(yintercept=TX2) }

## Plot data and curve. Plot everything in terms o f measurements ## relative to the corresponding control value. St art with just ## a linear model fit of log-transformed data, then transform estimates ## and confidence intervals back to the original sc ale.

t.lm <- lm(l10T4 ~ -1 + factor(block) + factor(dose ), data=rdta) ## Build a matrix of medians and confidence interva ls ## 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 varia ble t.relests <- as.data.frame(cbind(t.relests, block=c(2,3,2,3,3,3,2), dose=as.numeric(sub("^factor\\(d ose\\)","", 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 cur ve 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, col our=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.8 0), 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 Con trol") + 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 carrie d 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 co nfidence limits. ## The Bayes method has the advantage that we can c ompute 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 attribute s. ## ... 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","Mo del A2", "Pooled_chain3. txt"), file.path(".","Bayes","Mo del A2", "Pooled_index.t xt"), 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.tx t"), file.path(".","Bayes"," Model A3",

### 47

"s2_index.txt")

, thin=10,quiet=TRUE), read.coda(file.path(".","Bayes","Mo del A3","s2_chain2.txt"), file.path(".","Bayes","Mo del A3","s2_index.txt"), thin=10,quiet=TRUE), read.coda(file.path(".","Bayes","Mo del A3","s2_chain3.txt"), file.path(".","Bayes","Mo del 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 stru cture. 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 popula tion parameters by ## sampling repeatedly (so that the maximum se of t he 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,]),med ian(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 est imate) with CIs, ## along with block-specific dr curve, population d r 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(dos e), data=dta) ## Build a matrix of medians and confidence interva ls ## 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 varia ble t.relests2 <- as.data.frame(cbind(t.relests2, block=c(2,3,2,3,3,1,3,1,2), dose=as.numeric(sub("^factor\\(d ose\\)","", 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 valu e 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 es timates 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, r eplace=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, co lour=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 Con trol") ## Figure 2 png(file="RelChangefromControl.png", width=3.18*gol den, 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 p ercentiles 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()