Stacks.cdc.gov



Supplementary Material to “QRAD (Quantal Risk Assessment Database): a database for exploring patterns in quantal dose-response data in risk assessment and its application to develop priors for Bayesian dose-response analysis”Matthew W. Wheeler1*, Walter W. Piegorsch2, and A. John Bailer3This supplementary document provides supporting material for QRAD (Quantal Risk Assessment Database): a database for exploring patterns in quantal dose-response data in risk assessment and its application to develop priors for Bayesian dose-response analysis by M.H. Wheeler, et al., (the “main document”). The various sections below address a variety of supplemental/supporting topics, and are not intended to flow naturally between each other. They are, however, presented in roughly the same order in which their counterpart topics appear in the main document. We assume the database has been imported as an R object into an R workspace, or as a csv file into SAS?.S.1.Basic summaries: simple frequency tables and data sourcesFor a further illustration of QRAD’s use, consider construction of a two-way frequency table summarizing features of the data sources, say, numbers of dose groups (as in main document Fig. 1) vs. data source (i.e., IRIS, CalEPA, etc.). We start by constructing a data frame that has a variable with the number of dose groups in each experiment (ID) and a variable describing the data source. The R dplyr package ADDIN EN.CITE <EndNote><Cite><Author>Wickham</Author><Year>2016</Year><RecNum>5396</RecNum><DisplayText>(Wickham and Francois, 2016)</DisplayText><record><rec-number>5396</rec-number><foreign-keys><key app="EN" db-id="tdvd0x2f0rt2d1e22rmxew0pwtdwaev0w5sx">5396</key></foreign-keys><ref-type name="Web Page">12</ref-type><contributors><authors><author>Wickham, Hadley</author><author>Francois, Romain</author></authors></contributors><titles><title><style face="italic" font="default" size="100%">dplyr</style><style face="normal" font="default" size="100%">: A grammar of data manipulation. R package version 0.5.0. computing</keyword><keyword>R package</keyword></keywords><dates><year>2016</year></dates><urls><related-urls><url>;(Wickham and Francois, 2016) is again our primary tool:Summary2 <- final.data %>% group_by(ID, data.source) %>% summarise(nDoseGrps = n())The resulting two-way frequency table is constructed via the table function, here,group_table <- with(Summary2, table(data.source, nDoseGrps))print( group_table )which results in the following cross-classified table (output edited): nDoseGrpsdata.source 2 3 4 5 6 7 8 10 13 16 IRIS 27 116 155 36 30 3 0 3 1 0 CalEPA 61 83 71 12 3 1 1 0 1 2 PPRTV 2 26 41 1 5 0 0 0 0 0 HEAST 5 6 10 0 4 0 0 0 0 0 ATSDR 0 1 2 0 0 0 0 0 0 0 OPP 0 1 14 9 0 0 0 0 0 0From this tabular output we see that only the IRIS and CalEPA data sources contributed experiments with more than 6 dose groups to the database. The fewest studies were contributions by ATSDR. The (sorted) row margins can be easily obtained as sort(margin.table(group_table,1))producing (output edited)data.sourceATSDR OPP HEAST PPRTV CalEPA IRIS 3 24 25 75 235 371 which displays the relative contribution of the different data sources to the database in increasing order. Namely, ATSDR only contributed 3 studies to this database while IRIS contributed 371. As an aside, the frequency output / numerator of the relative frequencies in main document §2.3 on number of IDs can be produced from this contingency table via the commandmargin.table(group.table,2).S.2.Subsetting data to only have more than two dose groupsReturning to the Summary1 object from the main document’s §2.3, we find that 95 of the 733 of the quantal data sets (approximately 13%) in our database contain only two doses:table(Summary1$nDoseGrps)?2 ??3 ??4 ??5 ??6 ??7 ??8 ?10 ?13 ?1695 233 293 ?58 ?42 ??4 ??1 ??3 ??2 ??2In effect, these two-dose data sets are precluded from use in studying a dose-response pattern, although we feel they are still valuable as indicators of a toxicological effect in the experiments from which they were drawn. For investigators wishing to study dose responses in data sets with only three or more doses, those data sets with two doses can be ignored or removed. For instance, the R function filter (part of the dplyr package) can be used to select data sets with more than two doses. Now, we select these data sets with a row for each unique dose group for a particular experiment using the R code below, Summary3 <- final.data %>% group_by(ID) %>% mutate(nDoseGrps = n()) %>% select(ID, chemical, nDoseGrps, dose, r.dose, x, n) %>% filter(nDoseGrps > 2)One can easily confirm that chemicals with 3 or more dose groups are selected, e.g., viasort(unique(Summary3$nDoseGrps))we find the following numbers of dose groups in the subset database (output edited): 3 4 5 6 7 8 10 13 16Further operations on this new data frame can be conducted in similar fashion to those seen above.S.3.Applying a function to all experiments in the databaseFor a more-involved example of the type of toxicological calculations we might consider for these data, we fit the simple logistic regression model corresponding to main document Equation (6) to each data set in QRAD and then calculate a median effective dose, ED50, for the fit PEVuZE5vdGU+PENpdGU+PEF1dGhvcj5QaWVnb3JzY2g8L0F1dGhvcj48WWVhcj4yMDA1PC9ZZWFy

PjxSZWNOdW0+MTA3ODU8L1JlY051bT48U3VmZml4PmAsIMKnNC4xLjE8L1N1ZmZpeD48RGlzcGxh

eVRleHQ+KFBpZWdvcnNjaCBhbmQgQmFpbGVyLCAyMDA1LCDCpzQuMS4xKTwvRGlzcGxheVRleHQ+

PHJlY29yZD48cmVjLW51bWJlcj4xMDc4NTwvcmVjLW51bWJlcj48Zm9yZWlnbi1rZXlzPjxrZXkg

YXBwPSJFTiIgZGItaWQ9InRkdmQweDJmMHJ0MmQxZTIycm14ZXcwcHd0ZHdhZXYwdzVzeCI+MTA3

ODU8L2tleT48L2ZvcmVpZ24ta2V5cz48cmVmLXR5cGUgbmFtZT0iQm9vayI+NjwvcmVmLXR5cGU+

PGNvbnRyaWJ1dG9ycz48YXV0aG9ycz48YXV0aG9yPlBpZWdvcnNjaCwgV2FsdGVyIFc8L2F1dGhv

cj48YXV0aG9yPkJhaWxlciwgQSBKPC9hdXRob3I+PC9hdXRob3JzPjwvY29udHJpYnV0b3JzPjx0

aXRsZXM+PHRpdGxlPkFuYWx5emluZyBFbnZpcm9ubWVudGFsIERhdGE8L3RpdGxlPjwvdGl0bGVz

PjxrZXl3b3Jkcz48a2V5d29yZD5lbnZpcm9ubWV0cmljczwva2V5d29yZD48a2V5d29yZD5tdWx0

aXBsZSBsaW5lYXIgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5BTk9WQTwva2V5d29yZD48

a2V5d29yZD5BTkNPVkE8L2tleXdvcmQ+PGtleXdvcmQ+cmFuZG9tLWVmZmVjdHMgbW9kZWxzPC9r

ZXl3b3JkPjxrZXl3b3JkPnBvbHlub21pYWwgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5u

b25saW5lYXIgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5waWVjZXdpc2UgcmVncmVzc2lv

biBtb2RlbHM8L2tleXdvcmQ+PGtleXdvcmQ+Z3Jvd3RoIGN1cnZlczwva2V5d29yZD48a2V5d29y

ZD5nb21wZXJ0eiBtb2RlbDwva2V5d29yZD48a2V5d29yZD5sb2dpc3RpYyBncm93dGggY3VydmVz

PC9rZXl3b3JkPjxrZXl3b3JkPndlaWJ1bGwgZ3Jvd3RoIGN1cnZlczwva2V5d29yZD48a2V5d29y

ZD5yYXRpb25hbCBwb2x5bm9taWFsczwva2V5d29yZD48a2V5d29yZD5NaWNoYWVsaXMtTWVudGVu

IG1vZGVsPC9rZXl3b3JkPjxrZXl3b3JkPk1vcmdhbi1NZXJjZXItRmxvZGluIG1vZGVsPC9rZXl3

b3JkPjxrZXl3b3JkPmdlbmVyYWxpemVkIGxpbmVhciBtb2RlbHM8L2tleXdvcmQ+PGtleXdvcmQ+

ZGV2aWFuY2UgZnVuY3Rpb248L2tleXdvcmQ+PGtleXdvcmQ+bWF4aW11bSBxdWFzaS1saWtlbGlo

b29kPC9rZXl3b3JkPjxrZXl3b3JkPk1RTDwva2V5d29yZD48a2V5d29yZD5HRUU8L2tleXdvcmQ+

PGtleXdvcmQ+bG9naXN0aWMgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5vdmVyZGlzcGVy

c2lvbjwva2V5d29yZD48a2V5d29yZD5leHRyYS1iaW5vbWlhbCB2YXJpYWJpbGl0eTwva2V5d29y

ZD48a2V5d29yZD5jb3VudCBkYXRhPC9rZXl3b3JkPjxrZXl3b3JkPmV4dHJhLXBvaXNzb24gdmFy

aWFiaWxpdHk8L2tleXdvcmQ+PGtleXdvcmQ+Z2FtbWEgR0xpTXM8L2tleXdvcmQ+PGtleXdvcmQ+

cXVhbnRpdGF0aXZlIHJpc2sgYXNzZXNzbWVudDwva2V5d29yZD48a2V5d29yZD5wb3RlbmN5IGVz

dGltYXRpb248L2tleXdvcmQ+PGtleXdvcmQ+ZG9zZS1yZXNwb25zZSBkYXRhPC9rZXl3b3JkPjxr

ZXl3b3JkPm1lZGlhbiBlZmZlY3RpdmUgZG9zZTwva2V5d29yZD48a2V5d29yZD5hZGRpdGlvbmFs

IHJpc2s8L2tleXdvcmQ+PGtleXdvcmQ+ZXhjZXNzIHJpc2s8L2tleXdvcmQ+PGtleXdvcmQ+cmlz

ayBhbmFseXNpczwva2V5d29yZD48a2V5d29yZD5RUkE8L2tleXdvcmQ+PGtleXdvcmQ+YmVuY2ht

YXJrIGFuYWx5c2lzPC9rZXl3b3JkPjxrZXl3b3JkPmJlbmNobWFyayBkb3NlPC9rZXl3b3JkPjxr

ZXl3b3JkPkJNREw8L2tleXdvcmQ+PGtleXdvcmQ+dW5jZXJ0YWludHkgYW5hbHlzaXM8L2tleXdv

cmQ+PGtleXdvcmQ+dW5jZXJ0YWludHkgZmFjdG9yczwva2V5d29yZD48a2V5d29yZD5zZW5zaXRp

dml0eSBhbmFseXNpczwva2V5d29yZD48a2V5d29yZD5jb3JyZWxhdGlvbiByYXRpb3M8L2tleXdv

cmQ+PGtleXdvcmQ+dGVtcG9yYWwgZGF0YTwva2V5d29yZD48a2V5d29yZD5hdXRvcmVncmVzc2l2

ZSBtb2RlbGluZzwva2V5d29yZD48a2V5d29yZD50aW1lIHNlcmllczwva2V5d29yZD48a2V5d29y

ZD5oYXJtb25pYyByZWdyZXNzaW9uPC9rZXl3b3JkPjxrZXl3b3JkPkZvdXJpZXIgYW5hbHlzaXM8

L2tleXdvcmQ+PGtleXdvcmQ+YXV0b2NvcnJlbGF0aW9uIGZ1bmN0aW9uPC9rZXl3b3JkPjxrZXl3

b3JkPkFDRjwva2V5d29yZD48a2V5d29yZD5BUiBtb2RlbDwva2V5d29yZD48a2V5d29yZD5pbnRl

cnZlbnRpb24gYW5hbHlzaXM8L2tleXdvcmQ+PGtleXdvcmQ+c2Vhc29uYWxpdHk8L2tleXdvcmQ+

PGtleXdvcmQ+bG9uZ2l0dWRpbmFsIGdyb3d0aCBkYXRhPC9rZXl3b3JkPjxrZXl3b3JkPm1peGVk

IG1vZGVsczwva2V5d29yZD48a2V5d29yZD5zcGF0aWFsbHkgY29ycmVsYXRlZCBkYXRhPC9rZXl3

b3JkPjxrZXl3b3JkPnNwYXRpYWwgcG9pbnQgcGF0dGVybnM8L2tleXdvcmQ+PGtleXdvcmQ+Y29t

cGxldGUgc3BhdGlhbCByYW5kb21uZXNzPC9rZXl3b3JkPjxrZXl3b3JkPlJpcGxleeKAmXMgSyBm

dW5jdGlvbjwva2V5d29yZD48a2V5d29yZD5zcGF0aWFsIGF1dG9jb3JyZWxhdGlvbjwva2V5d29y

ZD48a2V5d29yZD5Nb3JhbiZhcG9zO3MgSSBjb2VmZmljaWVudDwva2V5d29yZD48a2V5d29yZD5H

ZWFyeSZhcG9zO3MgYyBjb2VmZmljaWVudDwva2V5d29yZD48a2V5d29yZD5zZW1pdmFyaW9ncmFt

PC9rZXl3b3JkPjxrZXl3b3JkPmtyaWdpbmc8L2tleXdvcmQ+PGtleXdvcmQ+c3BhdGlhbCBwcmVk

aWN0aW9uPC9rZXl3b3JkPjxrZXl3b3JkPmNvbWJpbmluZyBlbnZpcm9ubWVudGFsIGluZm9ybWF0

aW9uPC9rZXl3b3JkPjxrZXl3b3JkPmNvbWJpbmluZyBwLXZhbHVlczwva2V5d29yZD48a2V5d29y

ZD5lZmZlY3Qgc2l6ZSBlc3RpbWF0aW9uPC9rZXl3b3JkPjxrZXl3b3JkPm1ldGEtYW5hbHlzaXM8

L2tleXdvcmQ+PGtleXdvcmQ+cHVibGljYXRpb24gYmlhczwva2V5d29yZD48a2V5d29yZD5oaXN0

b3JpY2FsIGNvbnRyb2xzPC9rZXl3b3JkPjxrZXl3b3JkPmVudmlyb25tZW50YWwgc2FtcGxpbmcg

c3VydmV5IHNhbXBsaW5nPC9rZXl3b3JkPjxrZXl3b3JkPnN5c3RlbWF0aWMgc2FtcGxpbmc8L2tl

eXdvcmQ+PGtleXdvcmQ+c3RyYXRpZmllZCByYW5kb20gc2FtcGxpbmc8L2tleXdvcmQ+PGtleXdv

cmQ+Y2x1c3RlciBzYW1wbGluZzwva2V5d29yZD48a2V5d29yZD5jYXB0dXJlLXJlY2FwdHVyZSBz

YW1wbGluZzwva2V5d29yZD48a2V5d29yZD5xdWFkcmF0IHNhbXBsaW5nPC9rZXl3b3JkPjxrZXl3

b3JkPmxpbmUtaW50ZXJjZXB0IHNhbXBsaW5nPC9rZXl3b3JkPjxrZXl3b3JkPnJhbmtlZCBzZXQg

c2FtcGxpbmc8L2tleXdvcmQ+PGtleXdvcmQ+Y29tcG9zaXRlIHNhbXBsaW5nPC9rZXl3b3JkPjxr

ZXl3b3JkPnRhZy1yZWNhcHR1cmUgc2FtcGxpbmc8L2tleXdvcmQ+PC9rZXl3b3Jkcz48ZGF0ZXM+

PHllYXI+MjAwNTwveWVhcj48L2RhdGVzPjxwdWItbG9jYXRpb24+Q2hpY2hlc3RlcjwvcHViLWxv

Y2F0aW9uPjxwdWJsaXNoZXI+Sm9obiBXaWxleSAmYW1wOyBTb25zPC9wdWJsaXNoZXI+PGlzYm4+

MC00NzAtODQ4MzYtNyAoYWNpZC1mcmVlIHBhcGVyKSYjeEQ7OTc4LTAtNDcwLTg0ODM2LTQgKGFj

aWQtZnJlZSBwYXBlcikmI3hEOzA0NzAwMTIyMzQgKGVsZWN0cm9uaWMgYmsuKSYjeEQ7OTc4LTAt

NDcwLTAxMjIzLTkgKGVsZWN0cm9uaWMgYmsuKTwvaXNibj48Y2FsbC1udW0+R0U0NS5TMjUgUDU0

PC9jYWxsLW51bT48bGFiZWw+QjUmI3hEOyAgcGRmPC9sYWJlbD48dXJscz48cmVsYXRlZC11cmxz

Pjx1cmw+aHR0cDovL3d3dy53aWxleS5jb20vV2lsZXlDREEvV2lsZXlUaXRsZS9wcm9kdWN0Q2Qt

MDQ3MDg0ODM2Ny5odG1sPC91cmw+PHVybD5odHRwOi8vd3d3My5pbnRlcnNjaWVuY2Uud2lsZXku

Y29tL2NnaS1iaW4vYm9va2hvbWUvMTEyMTM0NjM2PC91cmw+PHVybD5odHRwOi8vd3d3My5pbnRl

cnNjaWVuY2Uud2lsZXkuY29tLmV6cHJveHkubGlicmFyeS5hcml6b25hLmVkdS9jZ2ktYmluL2Jv

b2tob21lLzExMjEzNDYzNjwvdXJsPjwvcmVsYXRlZC11cmxzPjwvdXJscz48L3JlY29yZD48L0Np

dGU+PC9FbmROb3RlPgB=

ADDIN EN.CITE PEVuZE5vdGU+PENpdGU+PEF1dGhvcj5QaWVnb3JzY2g8L0F1dGhvcj48WWVhcj4yMDA1PC9ZZWFy

PjxSZWNOdW0+MTA3ODU8L1JlY051bT48U3VmZml4PmAsIMKnNC4xLjE8L1N1ZmZpeD48RGlzcGxh

eVRleHQ+KFBpZWdvcnNjaCBhbmQgQmFpbGVyLCAyMDA1LCDCpzQuMS4xKTwvRGlzcGxheVRleHQ+

PHJlY29yZD48cmVjLW51bWJlcj4xMDc4NTwvcmVjLW51bWJlcj48Zm9yZWlnbi1rZXlzPjxrZXkg

YXBwPSJFTiIgZGItaWQ9InRkdmQweDJmMHJ0MmQxZTIycm14ZXcwcHd0ZHdhZXYwdzVzeCI+MTA3

ODU8L2tleT48L2ZvcmVpZ24ta2V5cz48cmVmLXR5cGUgbmFtZT0iQm9vayI+NjwvcmVmLXR5cGU+

PGNvbnRyaWJ1dG9ycz48YXV0aG9ycz48YXV0aG9yPlBpZWdvcnNjaCwgV2FsdGVyIFc8L2F1dGhv

cj48YXV0aG9yPkJhaWxlciwgQSBKPC9hdXRob3I+PC9hdXRob3JzPjwvY29udHJpYnV0b3JzPjx0

aXRsZXM+PHRpdGxlPkFuYWx5emluZyBFbnZpcm9ubWVudGFsIERhdGE8L3RpdGxlPjwvdGl0bGVz

PjxrZXl3b3Jkcz48a2V5d29yZD5lbnZpcm9ubWV0cmljczwva2V5d29yZD48a2V5d29yZD5tdWx0

aXBsZSBsaW5lYXIgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5BTk9WQTwva2V5d29yZD48

a2V5d29yZD5BTkNPVkE8L2tleXdvcmQ+PGtleXdvcmQ+cmFuZG9tLWVmZmVjdHMgbW9kZWxzPC9r

ZXl3b3JkPjxrZXl3b3JkPnBvbHlub21pYWwgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5u

b25saW5lYXIgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5waWVjZXdpc2UgcmVncmVzc2lv

biBtb2RlbHM8L2tleXdvcmQ+PGtleXdvcmQ+Z3Jvd3RoIGN1cnZlczwva2V5d29yZD48a2V5d29y

ZD5nb21wZXJ0eiBtb2RlbDwva2V5d29yZD48a2V5d29yZD5sb2dpc3RpYyBncm93dGggY3VydmVz

PC9rZXl3b3JkPjxrZXl3b3JkPndlaWJ1bGwgZ3Jvd3RoIGN1cnZlczwva2V5d29yZD48a2V5d29y

ZD5yYXRpb25hbCBwb2x5bm9taWFsczwva2V5d29yZD48a2V5d29yZD5NaWNoYWVsaXMtTWVudGVu

IG1vZGVsPC9rZXl3b3JkPjxrZXl3b3JkPk1vcmdhbi1NZXJjZXItRmxvZGluIG1vZGVsPC9rZXl3

b3JkPjxrZXl3b3JkPmdlbmVyYWxpemVkIGxpbmVhciBtb2RlbHM8L2tleXdvcmQ+PGtleXdvcmQ+

ZGV2aWFuY2UgZnVuY3Rpb248L2tleXdvcmQ+PGtleXdvcmQ+bWF4aW11bSBxdWFzaS1saWtlbGlo

b29kPC9rZXl3b3JkPjxrZXl3b3JkPk1RTDwva2V5d29yZD48a2V5d29yZD5HRUU8L2tleXdvcmQ+

PGtleXdvcmQ+bG9naXN0aWMgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5vdmVyZGlzcGVy

c2lvbjwva2V5d29yZD48a2V5d29yZD5leHRyYS1iaW5vbWlhbCB2YXJpYWJpbGl0eTwva2V5d29y

ZD48a2V5d29yZD5jb3VudCBkYXRhPC9rZXl3b3JkPjxrZXl3b3JkPmV4dHJhLXBvaXNzb24gdmFy

aWFiaWxpdHk8L2tleXdvcmQ+PGtleXdvcmQ+Z2FtbWEgR0xpTXM8L2tleXdvcmQ+PGtleXdvcmQ+

cXVhbnRpdGF0aXZlIHJpc2sgYXNzZXNzbWVudDwva2V5d29yZD48a2V5d29yZD5wb3RlbmN5IGVz

dGltYXRpb248L2tleXdvcmQ+PGtleXdvcmQ+ZG9zZS1yZXNwb25zZSBkYXRhPC9rZXl3b3JkPjxr

ZXl3b3JkPm1lZGlhbiBlZmZlY3RpdmUgZG9zZTwva2V5d29yZD48a2V5d29yZD5hZGRpdGlvbmFs

IHJpc2s8L2tleXdvcmQ+PGtleXdvcmQ+ZXhjZXNzIHJpc2s8L2tleXdvcmQ+PGtleXdvcmQ+cmlz

ayBhbmFseXNpczwva2V5d29yZD48a2V5d29yZD5RUkE8L2tleXdvcmQ+PGtleXdvcmQ+YmVuY2ht

YXJrIGFuYWx5c2lzPC9rZXl3b3JkPjxrZXl3b3JkPmJlbmNobWFyayBkb3NlPC9rZXl3b3JkPjxr

ZXl3b3JkPkJNREw8L2tleXdvcmQ+PGtleXdvcmQ+dW5jZXJ0YWludHkgYW5hbHlzaXM8L2tleXdv

cmQ+PGtleXdvcmQ+dW5jZXJ0YWludHkgZmFjdG9yczwva2V5d29yZD48a2V5d29yZD5zZW5zaXRp

dml0eSBhbmFseXNpczwva2V5d29yZD48a2V5d29yZD5jb3JyZWxhdGlvbiByYXRpb3M8L2tleXdv

cmQ+PGtleXdvcmQ+dGVtcG9yYWwgZGF0YTwva2V5d29yZD48a2V5d29yZD5hdXRvcmVncmVzc2l2

ZSBtb2RlbGluZzwva2V5d29yZD48a2V5d29yZD50aW1lIHNlcmllczwva2V5d29yZD48a2V5d29y

ZD5oYXJtb25pYyByZWdyZXNzaW9uPC9rZXl3b3JkPjxrZXl3b3JkPkZvdXJpZXIgYW5hbHlzaXM8

L2tleXdvcmQ+PGtleXdvcmQ+YXV0b2NvcnJlbGF0aW9uIGZ1bmN0aW9uPC9rZXl3b3JkPjxrZXl3

b3JkPkFDRjwva2V5d29yZD48a2V5d29yZD5BUiBtb2RlbDwva2V5d29yZD48a2V5d29yZD5pbnRl

cnZlbnRpb24gYW5hbHlzaXM8L2tleXdvcmQ+PGtleXdvcmQ+c2Vhc29uYWxpdHk8L2tleXdvcmQ+

PGtleXdvcmQ+bG9uZ2l0dWRpbmFsIGdyb3d0aCBkYXRhPC9rZXl3b3JkPjxrZXl3b3JkPm1peGVk

IG1vZGVsczwva2V5d29yZD48a2V5d29yZD5zcGF0aWFsbHkgY29ycmVsYXRlZCBkYXRhPC9rZXl3

b3JkPjxrZXl3b3JkPnNwYXRpYWwgcG9pbnQgcGF0dGVybnM8L2tleXdvcmQ+PGtleXdvcmQ+Y29t

cGxldGUgc3BhdGlhbCByYW5kb21uZXNzPC9rZXl3b3JkPjxrZXl3b3JkPlJpcGxleeKAmXMgSyBm

dW5jdGlvbjwva2V5d29yZD48a2V5d29yZD5zcGF0aWFsIGF1dG9jb3JyZWxhdGlvbjwva2V5d29y

ZD48a2V5d29yZD5Nb3JhbiZhcG9zO3MgSSBjb2VmZmljaWVudDwva2V5d29yZD48a2V5d29yZD5H

ZWFyeSZhcG9zO3MgYyBjb2VmZmljaWVudDwva2V5d29yZD48a2V5d29yZD5zZW1pdmFyaW9ncmFt

PC9rZXl3b3JkPjxrZXl3b3JkPmtyaWdpbmc8L2tleXdvcmQ+PGtleXdvcmQ+c3BhdGlhbCBwcmVk

aWN0aW9uPC9rZXl3b3JkPjxrZXl3b3JkPmNvbWJpbmluZyBlbnZpcm9ubWVudGFsIGluZm9ybWF0

aW9uPC9rZXl3b3JkPjxrZXl3b3JkPmNvbWJpbmluZyBwLXZhbHVlczwva2V5d29yZD48a2V5d29y

ZD5lZmZlY3Qgc2l6ZSBlc3RpbWF0aW9uPC9rZXl3b3JkPjxrZXl3b3JkPm1ldGEtYW5hbHlzaXM8

L2tleXdvcmQ+PGtleXdvcmQ+cHVibGljYXRpb24gYmlhczwva2V5d29yZD48a2V5d29yZD5oaXN0

b3JpY2FsIGNvbnRyb2xzPC9rZXl3b3JkPjxrZXl3b3JkPmVudmlyb25tZW50YWwgc2FtcGxpbmcg

c3VydmV5IHNhbXBsaW5nPC9rZXl3b3JkPjxrZXl3b3JkPnN5c3RlbWF0aWMgc2FtcGxpbmc8L2tl

eXdvcmQ+PGtleXdvcmQ+c3RyYXRpZmllZCByYW5kb20gc2FtcGxpbmc8L2tleXdvcmQ+PGtleXdv

cmQ+Y2x1c3RlciBzYW1wbGluZzwva2V5d29yZD48a2V5d29yZD5jYXB0dXJlLXJlY2FwdHVyZSBz

YW1wbGluZzwva2V5d29yZD48a2V5d29yZD5xdWFkcmF0IHNhbXBsaW5nPC9rZXl3b3JkPjxrZXl3

b3JkPmxpbmUtaW50ZXJjZXB0IHNhbXBsaW5nPC9rZXl3b3JkPjxrZXl3b3JkPnJhbmtlZCBzZXQg

c2FtcGxpbmc8L2tleXdvcmQ+PGtleXdvcmQ+Y29tcG9zaXRlIHNhbXBsaW5nPC9rZXl3b3JkPjxr

ZXl3b3JkPnRhZy1yZWNhcHR1cmUgc2FtcGxpbmc8L2tleXdvcmQ+PC9rZXl3b3Jkcz48ZGF0ZXM+

PHllYXI+MjAwNTwveWVhcj48L2RhdGVzPjxwdWItbG9jYXRpb24+Q2hpY2hlc3RlcjwvcHViLWxv

Y2F0aW9uPjxwdWJsaXNoZXI+Sm9obiBXaWxleSAmYW1wOyBTb25zPC9wdWJsaXNoZXI+PGlzYm4+

MC00NzAtODQ4MzYtNyAoYWNpZC1mcmVlIHBhcGVyKSYjeEQ7OTc4LTAtNDcwLTg0ODM2LTQgKGFj

aWQtZnJlZSBwYXBlcikmI3hEOzA0NzAwMTIyMzQgKGVsZWN0cm9uaWMgYmsuKSYjeEQ7OTc4LTAt

NDcwLTAxMjIzLTkgKGVsZWN0cm9uaWMgYmsuKTwvaXNibj48Y2FsbC1udW0+R0U0NS5TMjUgUDU0

PC9jYWxsLW51bT48bGFiZWw+QjUmI3hEOyAgcGRmPC9sYWJlbD48dXJscz48cmVsYXRlZC11cmxz

Pjx1cmw+aHR0cDovL3d3dy53aWxleS5jb20vV2lsZXlDREEvV2lsZXlUaXRsZS9wcm9kdWN0Q2Qt

MDQ3MDg0ODM2Ny5odG1sPC91cmw+PHVybD5odHRwOi8vd3d3My5pbnRlcnNjaWVuY2Uud2lsZXku

Y29tL2NnaS1iaW4vYm9va2hvbWUvMTEyMTM0NjM2PC91cmw+PHVybD5odHRwOi8vd3d3My5pbnRl

cnNjaWVuY2Uud2lsZXkuY29tLmV6cHJveHkubGlicmFyeS5hcml6b25hLmVkdS9jZ2ktYmluL2Jv

b2tob21lLzExMjEzNDYzNjwvdXJsPjwvcmVsYXRlZC11cmxzPjwvdXJscz48L3JlY29yZD48L0Np

dGU+PC9FbmROb3RlPgB=

ADDIN EN.CITE.DATA (Piegorsch and Bailer, 2005, §4.1.1). That is, given proportions, Yi/ni, and a series of doses, di (i = 1, ..., m), rescale the doses to ui = di/max{di}—the variable r.dose in the database—and assume a binomial model for the observed counts: Yi ~ indep. Bin.(ni,πi). Here the response probability πi is modeled as a function of the rescaled dose: πi = π(ui). The simple logistic dose-response model is π(u) = 1/(1 + exp{–?0 – ?1u}), with ED50 given by –?0/?1. The unknown ?-parameters and the ED50 are estimated from the data; maximum likelihood is a favored approach PEVuZE5vdGU+PENpdGU+PEF1dGhvcj5QaWVnb3JzY2g8L0F1dGhvcj48WWVhcj4yMDA1PC9ZZWFy

PjxSZWNOdW0+MTA3ODU8L1JlY051bT48U3VmZml4PmAsIMKnQS40LjM8L1N1ZmZpeD48RGlzcGxh

eVRleHQ+KFBpZWdvcnNjaCBhbmQgQmFpbGVyLCAyMDA1LCDCp0EuNC4zKTwvRGlzcGxheVRleHQ+

PHJlY29yZD48cmVjLW51bWJlcj4xMDc4NTwvcmVjLW51bWJlcj48Zm9yZWlnbi1rZXlzPjxrZXkg

YXBwPSJFTiIgZGItaWQ9InRkdmQweDJmMHJ0MmQxZTIycm14ZXcwcHd0ZHdhZXYwdzVzeCI+MTA3

ODU8L2tleT48L2ZvcmVpZ24ta2V5cz48cmVmLXR5cGUgbmFtZT0iQm9vayI+NjwvcmVmLXR5cGU+

PGNvbnRyaWJ1dG9ycz48YXV0aG9ycz48YXV0aG9yPlBpZWdvcnNjaCwgV2FsdGVyIFc8L2F1dGhv

cj48YXV0aG9yPkJhaWxlciwgQSBKPC9hdXRob3I+PC9hdXRob3JzPjwvY29udHJpYnV0b3JzPjx0

aXRsZXM+PHRpdGxlPkFuYWx5emluZyBFbnZpcm9ubWVudGFsIERhdGE8L3RpdGxlPjwvdGl0bGVz

PjxrZXl3b3Jkcz48a2V5d29yZD5lbnZpcm9ubWV0cmljczwva2V5d29yZD48a2V5d29yZD5tdWx0

aXBsZSBsaW5lYXIgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5BTk9WQTwva2V5d29yZD48

a2V5d29yZD5BTkNPVkE8L2tleXdvcmQ+PGtleXdvcmQ+cmFuZG9tLWVmZmVjdHMgbW9kZWxzPC9r

ZXl3b3JkPjxrZXl3b3JkPnBvbHlub21pYWwgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5u

b25saW5lYXIgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5waWVjZXdpc2UgcmVncmVzc2lv

biBtb2RlbHM8L2tleXdvcmQ+PGtleXdvcmQ+Z3Jvd3RoIGN1cnZlczwva2V5d29yZD48a2V5d29y

ZD5nb21wZXJ0eiBtb2RlbDwva2V5d29yZD48a2V5d29yZD5sb2dpc3RpYyBncm93dGggY3VydmVz

PC9rZXl3b3JkPjxrZXl3b3JkPndlaWJ1bGwgZ3Jvd3RoIGN1cnZlczwva2V5d29yZD48a2V5d29y

ZD5yYXRpb25hbCBwb2x5bm9taWFsczwva2V5d29yZD48a2V5d29yZD5NaWNoYWVsaXMtTWVudGVu

IG1vZGVsPC9rZXl3b3JkPjxrZXl3b3JkPk1vcmdhbi1NZXJjZXItRmxvZGluIG1vZGVsPC9rZXl3

b3JkPjxrZXl3b3JkPmdlbmVyYWxpemVkIGxpbmVhciBtb2RlbHM8L2tleXdvcmQ+PGtleXdvcmQ+

ZGV2aWFuY2UgZnVuY3Rpb248L2tleXdvcmQ+PGtleXdvcmQ+bWF4aW11bSBxdWFzaS1saWtlbGlo

b29kPC9rZXl3b3JkPjxrZXl3b3JkPk1RTDwva2V5d29yZD48a2V5d29yZD5HRUU8L2tleXdvcmQ+

PGtleXdvcmQ+bG9naXN0aWMgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5vdmVyZGlzcGVy

c2lvbjwva2V5d29yZD48a2V5d29yZD5leHRyYS1iaW5vbWlhbCB2YXJpYWJpbGl0eTwva2V5d29y

ZD48a2V5d29yZD5jb3VudCBkYXRhPC9rZXl3b3JkPjxrZXl3b3JkPmV4dHJhLXBvaXNzb24gdmFy

aWFiaWxpdHk8L2tleXdvcmQ+PGtleXdvcmQ+Z2FtbWEgR0xpTXM8L2tleXdvcmQ+PGtleXdvcmQ+

cXVhbnRpdGF0aXZlIHJpc2sgYXNzZXNzbWVudDwva2V5d29yZD48a2V5d29yZD5wb3RlbmN5IGVz

dGltYXRpb248L2tleXdvcmQ+PGtleXdvcmQ+ZG9zZS1yZXNwb25zZSBkYXRhPC9rZXl3b3JkPjxr

ZXl3b3JkPm1lZGlhbiBlZmZlY3RpdmUgZG9zZTwva2V5d29yZD48a2V5d29yZD5hZGRpdGlvbmFs

IHJpc2s8L2tleXdvcmQ+PGtleXdvcmQ+ZXhjZXNzIHJpc2s8L2tleXdvcmQ+PGtleXdvcmQ+cmlz

ayBhbmFseXNpczwva2V5d29yZD48a2V5d29yZD5RUkE8L2tleXdvcmQ+PGtleXdvcmQ+YmVuY2ht

YXJrIGFuYWx5c2lzPC9rZXl3b3JkPjxrZXl3b3JkPmJlbmNobWFyayBkb3NlPC9rZXl3b3JkPjxr

ZXl3b3JkPkJNREw8L2tleXdvcmQ+PGtleXdvcmQ+dW5jZXJ0YWludHkgYW5hbHlzaXM8L2tleXdv

cmQ+PGtleXdvcmQ+dW5jZXJ0YWludHkgZmFjdG9yczwva2V5d29yZD48a2V5d29yZD5zZW5zaXRp

dml0eSBhbmFseXNpczwva2V5d29yZD48a2V5d29yZD5jb3JyZWxhdGlvbiByYXRpb3M8L2tleXdv

cmQ+PGtleXdvcmQ+dGVtcG9yYWwgZGF0YTwva2V5d29yZD48a2V5d29yZD5hdXRvcmVncmVzc2l2

ZSBtb2RlbGluZzwva2V5d29yZD48a2V5d29yZD50aW1lIHNlcmllczwva2V5d29yZD48a2V5d29y

ZD5oYXJtb25pYyByZWdyZXNzaW9uPC9rZXl3b3JkPjxrZXl3b3JkPkZvdXJpZXIgYW5hbHlzaXM8

L2tleXdvcmQ+PGtleXdvcmQ+YXV0b2NvcnJlbGF0aW9uIGZ1bmN0aW9uPC9rZXl3b3JkPjxrZXl3

b3JkPkFDRjwva2V5d29yZD48a2V5d29yZD5BUiBtb2RlbDwva2V5d29yZD48a2V5d29yZD5pbnRl

cnZlbnRpb24gYW5hbHlzaXM8L2tleXdvcmQ+PGtleXdvcmQ+c2Vhc29uYWxpdHk8L2tleXdvcmQ+

PGtleXdvcmQ+bG9uZ2l0dWRpbmFsIGdyb3d0aCBkYXRhPC9rZXl3b3JkPjxrZXl3b3JkPm1peGVk

IG1vZGVsczwva2V5d29yZD48a2V5d29yZD5zcGF0aWFsbHkgY29ycmVsYXRlZCBkYXRhPC9rZXl3

b3JkPjxrZXl3b3JkPnNwYXRpYWwgcG9pbnQgcGF0dGVybnM8L2tleXdvcmQ+PGtleXdvcmQ+Y29t

cGxldGUgc3BhdGlhbCByYW5kb21uZXNzPC9rZXl3b3JkPjxrZXl3b3JkPlJpcGxleeKAmXMgSyBm

dW5jdGlvbjwva2V5d29yZD48a2V5d29yZD5zcGF0aWFsIGF1dG9jb3JyZWxhdGlvbjwva2V5d29y

ZD48a2V5d29yZD5Nb3JhbiZhcG9zO3MgSSBjb2VmZmljaWVudDwva2V5d29yZD48a2V5d29yZD5H

ZWFyeSZhcG9zO3MgYyBjb2VmZmljaWVudDwva2V5d29yZD48a2V5d29yZD5zZW1pdmFyaW9ncmFt

PC9rZXl3b3JkPjxrZXl3b3JkPmtyaWdpbmc8L2tleXdvcmQ+PGtleXdvcmQ+c3BhdGlhbCBwcmVk

aWN0aW9uPC9rZXl3b3JkPjxrZXl3b3JkPmNvbWJpbmluZyBlbnZpcm9ubWVudGFsIGluZm9ybWF0

aW9uPC9rZXl3b3JkPjxrZXl3b3JkPmNvbWJpbmluZyBwLXZhbHVlczwva2V5d29yZD48a2V5d29y

ZD5lZmZlY3Qgc2l6ZSBlc3RpbWF0aW9uPC9rZXl3b3JkPjxrZXl3b3JkPm1ldGEtYW5hbHlzaXM8

L2tleXdvcmQ+PGtleXdvcmQ+cHVibGljYXRpb24gYmlhczwva2V5d29yZD48a2V5d29yZD5oaXN0

b3JpY2FsIGNvbnRyb2xzPC9rZXl3b3JkPjxrZXl3b3JkPmVudmlyb25tZW50YWwgc2FtcGxpbmcg

c3VydmV5IHNhbXBsaW5nPC9rZXl3b3JkPjxrZXl3b3JkPnN5c3RlbWF0aWMgc2FtcGxpbmc8L2tl

eXdvcmQ+PGtleXdvcmQ+c3RyYXRpZmllZCByYW5kb20gc2FtcGxpbmc8L2tleXdvcmQ+PGtleXdv

cmQ+Y2x1c3RlciBzYW1wbGluZzwva2V5d29yZD48a2V5d29yZD5jYXB0dXJlLXJlY2FwdHVyZSBz

YW1wbGluZzwva2V5d29yZD48a2V5d29yZD5xdWFkcmF0IHNhbXBsaW5nPC9rZXl3b3JkPjxrZXl3

b3JkPmxpbmUtaW50ZXJjZXB0IHNhbXBsaW5nPC9rZXl3b3JkPjxrZXl3b3JkPnJhbmtlZCBzZXQg

c2FtcGxpbmc8L2tleXdvcmQ+PGtleXdvcmQ+Y29tcG9zaXRlIHNhbXBsaW5nPC9rZXl3b3JkPjxr

ZXl3b3JkPnRhZy1yZWNhcHR1cmUgc2FtcGxpbmc8L2tleXdvcmQ+PC9rZXl3b3Jkcz48ZGF0ZXM+

PHllYXI+MjAwNTwveWVhcj48L2RhdGVzPjxwdWItbG9jYXRpb24+Q2hpY2hlc3RlcjwvcHViLWxv

Y2F0aW9uPjxwdWJsaXNoZXI+Sm9obiBXaWxleSAmYW1wOyBTb25zPC9wdWJsaXNoZXI+PGlzYm4+

MC00NzAtODQ4MzYtNyAoYWNpZC1mcmVlIHBhcGVyKSYjeEQ7OTc4LTAtNDcwLTg0ODM2LTQgKGFj

aWQtZnJlZSBwYXBlcikmI3hEOzA0NzAwMTIyMzQgKGVsZWN0cm9uaWMgYmsuKSYjeEQ7OTc4LTAt

NDcwLTAxMjIzLTkgKGVsZWN0cm9uaWMgYmsuKTwvaXNibj48Y2FsbC1udW0+R0U0NS5TMjUgUDU0

PC9jYWxsLW51bT48bGFiZWw+QjUmI3hEOyAgcGRmPC9sYWJlbD48dXJscz48cmVsYXRlZC11cmxz

Pjx1cmw+aHR0cDovL3d3dy53aWxleS5jb20vV2lsZXlDREEvV2lsZXlUaXRsZS9wcm9kdWN0Q2Qt

MDQ3MDg0ODM2Ny5odG1sPC91cmw+PHVybD5odHRwOi8vd3d3My5pbnRlcnNjaWVuY2Uud2lsZXku

Y29tL2NnaS1iaW4vYm9va2hvbWUvMTEyMTM0NjM2PC91cmw+PHVybD5odHRwOi8vd3d3My5pbnRl

cnNjaWVuY2Uud2lsZXkuY29tLmV6cHJveHkubGlicmFyeS5hcml6b25hLmVkdS9jZ2ktYmluL2Jv

b2tob21lLzExMjEzNDYzNjwvdXJsPjwvcmVsYXRlZC11cmxzPjwvdXJscz48L3JlY29yZD48L0Np

dGU+PC9FbmROb3RlPgB=

ADDIN EN.CITE PEVuZE5vdGU+PENpdGU+PEF1dGhvcj5QaWVnb3JzY2g8L0F1dGhvcj48WWVhcj4yMDA1PC9ZZWFy

PjxSZWNOdW0+MTA3ODU8L1JlY051bT48U3VmZml4PmAsIMKnQS40LjM8L1N1ZmZpeD48RGlzcGxh

eVRleHQ+KFBpZWdvcnNjaCBhbmQgQmFpbGVyLCAyMDA1LCDCp0EuNC4zKTwvRGlzcGxheVRleHQ+

PHJlY29yZD48cmVjLW51bWJlcj4xMDc4NTwvcmVjLW51bWJlcj48Zm9yZWlnbi1rZXlzPjxrZXkg

YXBwPSJFTiIgZGItaWQ9InRkdmQweDJmMHJ0MmQxZTIycm14ZXcwcHd0ZHdhZXYwdzVzeCI+MTA3

ODU8L2tleT48L2ZvcmVpZ24ta2V5cz48cmVmLXR5cGUgbmFtZT0iQm9vayI+NjwvcmVmLXR5cGU+

PGNvbnRyaWJ1dG9ycz48YXV0aG9ycz48YXV0aG9yPlBpZWdvcnNjaCwgV2FsdGVyIFc8L2F1dGhv

cj48YXV0aG9yPkJhaWxlciwgQSBKPC9hdXRob3I+PC9hdXRob3JzPjwvY29udHJpYnV0b3JzPjx0

aXRsZXM+PHRpdGxlPkFuYWx5emluZyBFbnZpcm9ubWVudGFsIERhdGE8L3RpdGxlPjwvdGl0bGVz

PjxrZXl3b3Jkcz48a2V5d29yZD5lbnZpcm9ubWV0cmljczwva2V5d29yZD48a2V5d29yZD5tdWx0

aXBsZSBsaW5lYXIgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5BTk9WQTwva2V5d29yZD48

a2V5d29yZD5BTkNPVkE8L2tleXdvcmQ+PGtleXdvcmQ+cmFuZG9tLWVmZmVjdHMgbW9kZWxzPC9r

ZXl3b3JkPjxrZXl3b3JkPnBvbHlub21pYWwgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5u

b25saW5lYXIgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5waWVjZXdpc2UgcmVncmVzc2lv

biBtb2RlbHM8L2tleXdvcmQ+PGtleXdvcmQ+Z3Jvd3RoIGN1cnZlczwva2V5d29yZD48a2V5d29y

ZD5nb21wZXJ0eiBtb2RlbDwva2V5d29yZD48a2V5d29yZD5sb2dpc3RpYyBncm93dGggY3VydmVz

PC9rZXl3b3JkPjxrZXl3b3JkPndlaWJ1bGwgZ3Jvd3RoIGN1cnZlczwva2V5d29yZD48a2V5d29y

ZD5yYXRpb25hbCBwb2x5bm9taWFsczwva2V5d29yZD48a2V5d29yZD5NaWNoYWVsaXMtTWVudGVu

IG1vZGVsPC9rZXl3b3JkPjxrZXl3b3JkPk1vcmdhbi1NZXJjZXItRmxvZGluIG1vZGVsPC9rZXl3

b3JkPjxrZXl3b3JkPmdlbmVyYWxpemVkIGxpbmVhciBtb2RlbHM8L2tleXdvcmQ+PGtleXdvcmQ+

ZGV2aWFuY2UgZnVuY3Rpb248L2tleXdvcmQ+PGtleXdvcmQ+bWF4aW11bSBxdWFzaS1saWtlbGlo

b29kPC9rZXl3b3JkPjxrZXl3b3JkPk1RTDwva2V5d29yZD48a2V5d29yZD5HRUU8L2tleXdvcmQ+

PGtleXdvcmQ+bG9naXN0aWMgcmVncmVzc2lvbjwva2V5d29yZD48a2V5d29yZD5vdmVyZGlzcGVy

c2lvbjwva2V5d29yZD48a2V5d29yZD5leHRyYS1iaW5vbWlhbCB2YXJpYWJpbGl0eTwva2V5d29y

ZD48a2V5d29yZD5jb3VudCBkYXRhPC9rZXl3b3JkPjxrZXl3b3JkPmV4dHJhLXBvaXNzb24gdmFy

aWFiaWxpdHk8L2tleXdvcmQ+PGtleXdvcmQ+Z2FtbWEgR0xpTXM8L2tleXdvcmQ+PGtleXdvcmQ+

cXVhbnRpdGF0aXZlIHJpc2sgYXNzZXNzbWVudDwva2V5d29yZD48a2V5d29yZD5wb3RlbmN5IGVz

dGltYXRpb248L2tleXdvcmQ+PGtleXdvcmQ+ZG9zZS1yZXNwb25zZSBkYXRhPC9rZXl3b3JkPjxr

ZXl3b3JkPm1lZGlhbiBlZmZlY3RpdmUgZG9zZTwva2V5d29yZD48a2V5d29yZD5hZGRpdGlvbmFs

IHJpc2s8L2tleXdvcmQ+PGtleXdvcmQ+ZXhjZXNzIHJpc2s8L2tleXdvcmQ+PGtleXdvcmQ+cmlz

ayBhbmFseXNpczwva2V5d29yZD48a2V5d29yZD5RUkE8L2tleXdvcmQ+PGtleXdvcmQ+YmVuY2ht

YXJrIGFuYWx5c2lzPC9rZXl3b3JkPjxrZXl3b3JkPmJlbmNobWFyayBkb3NlPC9rZXl3b3JkPjxr

ZXl3b3JkPkJNREw8L2tleXdvcmQ+PGtleXdvcmQ+dW5jZXJ0YWludHkgYW5hbHlzaXM8L2tleXdv

cmQ+PGtleXdvcmQ+dW5jZXJ0YWludHkgZmFjdG9yczwva2V5d29yZD48a2V5d29yZD5zZW5zaXRp

dml0eSBhbmFseXNpczwva2V5d29yZD48a2V5d29yZD5jb3JyZWxhdGlvbiByYXRpb3M8L2tleXdv

cmQ+PGtleXdvcmQ+dGVtcG9yYWwgZGF0YTwva2V5d29yZD48a2V5d29yZD5hdXRvcmVncmVzc2l2

ZSBtb2RlbGluZzwva2V5d29yZD48a2V5d29yZD50aW1lIHNlcmllczwva2V5d29yZD48a2V5d29y

ZD5oYXJtb25pYyByZWdyZXNzaW9uPC9rZXl3b3JkPjxrZXl3b3JkPkZvdXJpZXIgYW5hbHlzaXM8

L2tleXdvcmQ+PGtleXdvcmQ+YXV0b2NvcnJlbGF0aW9uIGZ1bmN0aW9uPC9rZXl3b3JkPjxrZXl3

b3JkPkFDRjwva2V5d29yZD48a2V5d29yZD5BUiBtb2RlbDwva2V5d29yZD48a2V5d29yZD5pbnRl

cnZlbnRpb24gYW5hbHlzaXM8L2tleXdvcmQ+PGtleXdvcmQ+c2Vhc29uYWxpdHk8L2tleXdvcmQ+

PGtleXdvcmQ+bG9uZ2l0dWRpbmFsIGdyb3d0aCBkYXRhPC9rZXl3b3JkPjxrZXl3b3JkPm1peGVk

IG1vZGVsczwva2V5d29yZD48a2V5d29yZD5zcGF0aWFsbHkgY29ycmVsYXRlZCBkYXRhPC9rZXl3

b3JkPjxrZXl3b3JkPnNwYXRpYWwgcG9pbnQgcGF0dGVybnM8L2tleXdvcmQ+PGtleXdvcmQ+Y29t

cGxldGUgc3BhdGlhbCByYW5kb21uZXNzPC9rZXl3b3JkPjxrZXl3b3JkPlJpcGxleeKAmXMgSyBm

dW5jdGlvbjwva2V5d29yZD48a2V5d29yZD5zcGF0aWFsIGF1dG9jb3JyZWxhdGlvbjwva2V5d29y

ZD48a2V5d29yZD5Nb3JhbiZhcG9zO3MgSSBjb2VmZmljaWVudDwva2V5d29yZD48a2V5d29yZD5H

ZWFyeSZhcG9zO3MgYyBjb2VmZmljaWVudDwva2V5d29yZD48a2V5d29yZD5zZW1pdmFyaW9ncmFt

PC9rZXl3b3JkPjxrZXl3b3JkPmtyaWdpbmc8L2tleXdvcmQ+PGtleXdvcmQ+c3BhdGlhbCBwcmVk

aWN0aW9uPC9rZXl3b3JkPjxrZXl3b3JkPmNvbWJpbmluZyBlbnZpcm9ubWVudGFsIGluZm9ybWF0

aW9uPC9rZXl3b3JkPjxrZXl3b3JkPmNvbWJpbmluZyBwLXZhbHVlczwva2V5d29yZD48a2V5d29y

ZD5lZmZlY3Qgc2l6ZSBlc3RpbWF0aW9uPC9rZXl3b3JkPjxrZXl3b3JkPm1ldGEtYW5hbHlzaXM8

L2tleXdvcmQ+PGtleXdvcmQ+cHVibGljYXRpb24gYmlhczwva2V5d29yZD48a2V5d29yZD5oaXN0

b3JpY2FsIGNvbnRyb2xzPC9rZXl3b3JkPjxrZXl3b3JkPmVudmlyb25tZW50YWwgc2FtcGxpbmcg

c3VydmV5IHNhbXBsaW5nPC9rZXl3b3JkPjxrZXl3b3JkPnN5c3RlbWF0aWMgc2FtcGxpbmc8L2tl

eXdvcmQ+PGtleXdvcmQ+c3RyYXRpZmllZCByYW5kb20gc2FtcGxpbmc8L2tleXdvcmQ+PGtleXdv

cmQ+Y2x1c3RlciBzYW1wbGluZzwva2V5d29yZD48a2V5d29yZD5jYXB0dXJlLXJlY2FwdHVyZSBz

YW1wbGluZzwva2V5d29yZD48a2V5d29yZD5xdWFkcmF0IHNhbXBsaW5nPC9rZXl3b3JkPjxrZXl3

b3JkPmxpbmUtaW50ZXJjZXB0IHNhbXBsaW5nPC9rZXl3b3JkPjxrZXl3b3JkPnJhbmtlZCBzZXQg

c2FtcGxpbmc8L2tleXdvcmQ+PGtleXdvcmQ+Y29tcG9zaXRlIHNhbXBsaW5nPC9rZXl3b3JkPjxr

ZXl3b3JkPnRhZy1yZWNhcHR1cmUgc2FtcGxpbmc8L2tleXdvcmQ+PC9rZXl3b3Jkcz48ZGF0ZXM+

PHllYXI+MjAwNTwveWVhcj48L2RhdGVzPjxwdWItbG9jYXRpb24+Q2hpY2hlc3RlcjwvcHViLWxv

Y2F0aW9uPjxwdWJsaXNoZXI+Sm9obiBXaWxleSAmYW1wOyBTb25zPC9wdWJsaXNoZXI+PGlzYm4+

MC00NzAtODQ4MzYtNyAoYWNpZC1mcmVlIHBhcGVyKSYjeEQ7OTc4LTAtNDcwLTg0ODM2LTQgKGFj

aWQtZnJlZSBwYXBlcikmI3hEOzA0NzAwMTIyMzQgKGVsZWN0cm9uaWMgYmsuKSYjeEQ7OTc4LTAt

NDcwLTAxMjIzLTkgKGVsZWN0cm9uaWMgYmsuKTwvaXNibj48Y2FsbC1udW0+R0U0NS5TMjUgUDU0

PC9jYWxsLW51bT48bGFiZWw+QjUmI3hEOyAgcGRmPC9sYWJlbD48dXJscz48cmVsYXRlZC11cmxz

Pjx1cmw+aHR0cDovL3d3dy53aWxleS5jb20vV2lsZXlDREEvV2lsZXlUaXRsZS9wcm9kdWN0Q2Qt

MDQ3MDg0ODM2Ny5odG1sPC91cmw+PHVybD5odHRwOi8vd3d3My5pbnRlcnNjaWVuY2Uud2lsZXku

Y29tL2NnaS1iaW4vYm9va2hvbWUvMTEyMTM0NjM2PC91cmw+PHVybD5odHRwOi8vd3d3My5pbnRl

cnNjaWVuY2Uud2lsZXkuY29tLmV6cHJveHkubGlicmFyeS5hcml6b25hLmVkdS9jZ2ktYmluL2Jv

b2tob21lLzExMjEzNDYzNjwvdXJsPjwvcmVsYXRlZC11cmxzPjwvdXJscz48L3JlY29yZD48L0Np

dGU+PC9FbmROb3RlPgB=

ADDIN EN.CITE.DATA (Piegorsch and Bailer, 2005, §A.4.3). In R, we can fit the logistic model to each experiment via the glm function with its family= binomial option, using the rescaled doses r.dose. The estimated ED50 based on this fit is then the ratio –?0/?1 from the R output.We begin by creating a new data frame that includes two new variables: (i) nobs: a count of the number of observations in a dose group with the response and (ii) nNotObs: a count of the number of observations in a dose group without the response:model_final<- final.data %>%? mutate(nobs=obs, nNotObs=n-obs)Next, we split the data frame into lists by the original study’s ID and fit the logistic model separately to each study using R’s glm function:split_model_final <- split(model_final, model_final$ID) ?model_logistic_reg_fits <- lapply( split_model_final, ??????????????????????? function(x) glm(cbind(nobs, nNotObs)~r.dose, ????????????????????????????????data=x, family=binomial) )Note that this will generate a number of R errors and warnings (e.g., algorithm did not converge, or fitted probabilities numerically 0 or 1 occurred). Some of these data sets are not fit well by a linear logistic model—see below.We first extract the coefficients from each fit to estimate the ED50:model_coef <- lapply(model_logistic_reg_fits, coef) After the model coefficients are extracted from the fits, the ED50s are calculated and stored as a vector viaED50.est <- lapply(model_coef, function(x) -x[1]/x[2]) ?ED50_vec <- unlist(ED50.est)This produces 27 estimated ED50s less than zero and five ED50s greater than 5.0. Both reflect instability in the estimates since, as constructed here, (i)?a dose cannot drop below zero and (ii)?the doses have been scaled to have a maximum value of 1.0. For this analysis, we exclude any of these 32 possible values, which results in a new data frame containing 701 ED50 estimates. Viz.,ED50_tmp <- data.frame(ED50=ED50_vec)ED50_df <- ED50_tmp %>% filter(ED50>=0, ED50<=5)After excluding the missing ED50s in a new data frame noMissED50.df with complete information, we construct a histogram of the final, estimated ED50s in Fig. S.1 using R’s stock hist function. As might be expected, the sample histogram displays a substantial right skew. A similar graphic can be constructed via the ggplot2 function geom_histogram with binwidth specified. Sample R code is# distribution of ED50ggplot(ED50_df,aes(x=ED50)) + geom_histogram(color='darkgrey',fill='white', binwidth=0.25) + theme_minimal()Fig. S.1. Histogram of estimated median lethal doses (ED50) from QRAD.It is also illustrative to examine the various empirical dose-response patterns from the various studies, stratified by the different data sources. This can provide insights into problems encountered when attempting formal model fits such as a logistic regression. The dose-response profiles for each experiment were plotted as a line graph with the display faceted by the data sources. Sample ggplot commands to produce this display areggplot(final.data, aes(x=r.dose,y=obs/n, group=ID)) + ?geom_line(alpha=.2) + ?theme_bw() +?xlab('Scaled Dose') +?ylab('Proportion Responding') + facet_wrap(~data.source)The results appear in Fig. S.2. While many of the dose-response patterns therein appear monotonic and non-decreasing, there remain a large number of experiments with a drop or reduction in the response at higher dose levels, colorfully known as an ‘umbrella pattern’ ADDIN EN.CITE <EndNote><Cite><Author>Mack</Author><Year>1981</Year><RecNum>190</RecNum><DisplayText>(Mack and Wolfe, 1981)</DisplayText><record><rec-number>190</rec-number><foreign-keys><key app="EN" db-id="tdvd0x2f0rt2d1e22rmxew0pwtdwaev0w5sx">190</key></foreign-keys><ref-type name="Journal Article">17</ref-type><contributors><authors><author>Mack, G A</author><author>Wolfe, D A</author></authors></contributors><titles><title><style face="italic" font="default" size="100%">K</style><style face="normal" font="default" size="100%">-sample rank tests for umbrella alternatives</style></title><secondary-title>J. Amer. Statist. Assoc.</secondary-title><alt-title>Journal of the American Statistical Association</alt-title></titles><periodical><full-title>J. Amer. Statist. Assoc.</full-title><abbr-1>Journal of the American Statistical Association</abbr-1></periodical><alt-periodical><full-title>J. Amer. Statist. Assoc.</full-title><abbr-1>Journal of the American Statistical Association</abbr-1></alt-periodical><pages>175-181</pages><volume>76</volume><number>373</number><keywords><keyword>dose downturn</keyword><keyword>nonparametric methods</keyword><keyword>distribution-free tests</keyword><keyword>power</keyword><keyword>one-factor ANOVA</keyword></keywords><dates><year>1981</year></dates><label> pdf</label><urls><related-urls><url>;(Mack and Wolfe, 1981). This is often seen when dose-related toxicity or other competing risks lead to mortality in the experimental subject before a non-lethal outcome such as cancer is observed. Fig. S.2. QRAD dose-response patterns, stratified by data sourceWe can modify Fig. S.2 to highlight the patterns associated with a non-significant test of trend. In the last example of this section, we first fit a separate logistic regression to each experiment (ID) in the database. We extract the p-value associated with a test on the significance of the slope from the logistic regression and then highlight in color those dose-response patterns based on this quantity. To start, we use lapply to fit all of the logistic regressions and extract the coefficients and p-value (the 2nd row and 4th column of the matrix of results obtained from the coefficient extractor function):model_final<- final.data %>% mutate(nobs=obs, nNotObs=n-obs)split_model_final <- split(model_final, model_final$ID) model_logistic_reg_fits <- lapply( split_model_final, function(x) glm(cbind(nobs, nNotObs)~r.dose, data=x, family=binomial) )model_summ <- lapply(model_logistic_reg_fits, summary)model_coef <- lapply(model_summ, coefficients)model_Trend_Pvalue <- lapply(model_coef, function(x) x[2,4])These quantities are placed into a data frame that is then joined with the original database in preparation of plotting. We restrict attention in this plot to experiments with 3, 4, 5, or 6 dose groups.Data_3to6grps <- final.data %>% group_by(ID, data.source) %>% mutate(nDoseGrps = n()) %>% filter(nDoseGrps > 2 & nDoseGrps<7)# Add p-value from logistic regression fitData_3to6grps_plus <- left_join(Data_3to6grps,trend_df)Data_3to6grps_plus %>% select(ID, dose, r.dose, obs, n, nDoseGrps)For instance, we plot dose-response patterns with p-values from the logistic slope test larger than 0.10 for the CalEPA experiments using the following code:Data_3to6grps_plus %>% mutate(NSig = (Pvalue > 0.10 )) %>% filter(data.source=="CalEPA") %>% ggplot(aes(x=r.dose,y=obs/n, group=ID)) + geom_line(aes(color=NSig, width=2), alpha=.5) + theme_bw() + xlab('Scaled Dose') + ylab('Proportion Responding') + facet_wrap(~nDoseGrps) + scale_color_manual(values=c("light grey", "red"), guide=FALSE)which produces Fig. S.3.Fig. S.3. QRAD dose-response patterns, stratified by number of dose groups for CalEPA studies. Patterns with a p-value > 0.10 from a logistic regression test of the slope parameter are colored red in this plot.While some highlighted patterns in this plot are flat and consistent with what we would expect from a case where the response does not change over dose (e.g. 0/48, 1/49, 3/49), others are surprising (e.g., 0/19, 0/46, 19/45). These latter cases illustrate enigmatic issues with the logistic regression test when zero or flat response patterns are observed at all but the highest dose(s).S.4.Using EPA BMDS software to fit multistage modelsAs our database is based upon quantal-data studies employed in the regulatory process, interest may exist in comparing potential differences in Benchmark Dose (BMD) estimates ADDIN EN.CITE <EndNote><Cite><Author>Crump</Author><Year>1984</Year><RecNum>5371</RecNum><DisplayText>(Crump, 1984)</DisplayText><record><rec-number>5371</rec-number><foreign-keys><key app="EN" db-id="tdvd0x2f0rt2d1e22rmxew0pwtdwaev0w5sx">5371</key></foreign-keys><ref-type name="Journal Article">17</ref-type><contributors><authors><author>Crump, Kenny S</author></authors></contributors><titles><title>A new method for determining allowable daily intake</title><secondary-title>Fundam. Appl. Toxicol.</secondary-title><alt-title>Fundamental and Applied Toxicology</alt-title></titles><alt-periodical><full-title>Fund. Appl. Toxicol.</full-title><abbr-1>Fundamental and Applied Toxicology</abbr-1></alt-periodical><pages>854-871</pages><volume>4</volume><number>5</number><keywords><keyword>history of statistics</keyword><keyword>introduction of benchmark dose</keyword><keyword>BMD</keyword><keyword>multiple BMRs</keyword><keyword>BMD0.01</keyword><keyword>BMD01</keyword><keyword>BMD05</keyword><keyword>BMD10</keyword><keyword>quantal response data (Table 1)</keyword><keyword>benchmark analysis</keyword><keyword>QRA</keyword><keyword>quantitative risk assessment</keyword><keyword>low-dose extrapolation</keyword><keyword>MTDADI</keyword><keyword>RfC</keyword><keyword>RfD</keyword></keywords><dates><year>1984</year></dates><label> pdf</label><urls><related-urls><url>(84)90107-6</url></related-urls></urls><custom1>PMID6510615</custom1></record></Cite></EndNote>(Crump, 1984) and in their corresponding 100(1–α)% lower bounds (BMDLs), which are often starting points in quantitative risk assessments PEVuZE5vdGU+PENpdGU+PEF1dGhvcj5JemFkaTwvQXV0aG9yPjxZZWFyPjIwMTI8L1llYXI+PFJl

Y051bT4zNjE3PC9SZWNOdW0+PERpc3BsYXlUZXh0PihJemFkaTxzdHlsZSBmYWNlPSJpdGFsaWMi

PiBldCBhbC48L3N0eWxlPiwgMjAxMjsgVS5LLiBDb21taXR0ZWUgb24gQ2FyY2lub2dlbmljaXR5

IChDT0MpLCAyMDE0KTwvRGlzcGxheVRleHQ+PHJlY29yZD48cmVjLW51bWJlcj4zNjE3PC9yZWMt

bnVtYmVyPjxmb3JlaWduLWtleXM+PGtleSBhcHA9IkVOIiBkYi1pZD0idGR2ZDB4MmYwcnQyZDFl

MjJybXhldzBwd3Rkd2FldjB3NXN4Ij4zNjE3PC9rZXk+PC9mb3JlaWduLWtleXM+PHJlZi10eXBl

IG5hbWU9IkpvdXJuYWwgQXJ0aWNsZSI+MTc8L3JlZi10eXBlPjxjb250cmlidXRvcnM+PGF1dGhv

cnM+PGF1dGhvcj5JemFkaSwgSG9kYTwvYXV0aG9yPjxhdXRob3I+R3J1bmR5LCBKZWFuIEU8L2F1

dGhvcj48YXV0aG9yPkJvc2UsIFJhbmphbjwvYXV0aG9yPjwvYXV0aG9ycz48L2NvbnRyaWJ1dG9y

cz48dGl0bGVzPjx0aXRsZT5FdmFsdWF0aW9uIG9mIHRoZSBiZW5jaG1hcmsgZG9zZSBmb3IgcG9p

bnQgb2YgZGVwYXJ0dXJlIGRldGVybWluYXRpb24gZm9yIGEgdmFyaWV0eSBvZiBjaGVtaWNhbCBj

bGFzc2VzIGluIGFwcGxpZWQgcmVndWxhdG9yeSBzZXR0aW5nczwvdGl0bGU+PHNlY29uZGFyeS10

aXRsZT5SaXNrIEFuYWx5Ljwvc2Vjb25kYXJ5LXRpdGxlPjxhbHQtdGl0bGU+UmlzayBBbmFseXNp

czwvYWx0LXRpdGxlPjwvdGl0bGVzPjxhbHQtcGVyaW9kaWNhbD48ZnVsbC10aXRsZT5SaXNrIEFu

YWwuPC9mdWxsLXRpdGxlPjxhYmJyLTE+UmlzayBBbmFseXNpczwvYWJici0xPjwvYWx0LXBlcmlv

ZGljYWw+PHBhZ2VzPjgzMC04MzU8L3BhZ2VzPjx2b2x1bWU+MzI8L3ZvbHVtZT48bnVtYmVyPjU8

L251bWJlcj48a2V5d29yZHM+PGtleXdvcmQ+YmVuY2htYXJrIGFuYWx5c2lzPC9rZXl3b3JkPjxr

ZXl3b3JkPkJNRDwva2V5d29yZD48a2V5d29yZD5CTUM8L2tleXdvcmQ+PGtleXdvcmQ+cG9pbnQg

b2YgZGVwYXJ0dXJlPC9rZXl3b3JkPjxrZXl3b3JkPlBPRDwva2V5d29yZD48a2V5d29yZD5yZWd1

bGF0b3J5IGFwcGxpY2F0aW9uczwva2V5d29yZD48L2tleXdvcmRzPjxkYXRlcz48eWVhcj4yMDEy

PC95ZWFyPjwvZGF0ZXM+PGxhYmVsPiAgcGRmPC9sYWJlbD48dXJscz48cmVsYXRlZC11cmxzPjx1

cmw+aHR0cDovL2R4LmRvaS5vcmcvMTAuMTExMS9qLjE1MzktNjkyNC4yMDExLjAxNzMyLng8L3Vy

bD48L3JlbGF0ZWQtdXJscz48L3VybHM+PC9yZWNvcmQ+PC9DaXRlPjxDaXRlPjxBdXRob3I+VS5L

LiBDb21taXR0ZWUgb24gQ2FyY2lub2dlbmljaXR5IChDT0MpPC9BdXRob3I+PFllYXI+MjAxNDwv

WWVhcj48UmVjTnVtPjEyMDQ2PC9SZWNOdW0+PHJlY29yZD48cmVjLW51bWJlcj4xMjA0NjwvcmVj

LW51bWJlcj48Zm9yZWlnbi1rZXlzPjxrZXkgYXBwPSJFTiIgZGItaWQ9InRkdmQweDJmMHJ0MmQx

ZTIycm14ZXcwcHd0ZHdhZXYwdzVzeCI+MTIwNDY8L2tleT48L2ZvcmVpZ24ta2V5cz48cmVmLXR5

cGUgbmFtZT0iUmVwb3J0Ij4yNzwvcmVmLXR5cGU+PGNvbnRyaWJ1dG9ycz48YXV0aG9ycz48YXV0

aG9yPlUuSy4gQ29tbWl0dGVlIG9uIENhcmNpbm9nZW5pY2l0eSAoQ09DKSw8L2F1dGhvcj48L2F1

dGhvcnM+PC9jb250cmlidXRvcnM+PHRpdGxlcz48dGl0bGU+RGVmaW5pbmcgYSBQb2ludCBvZiBE

ZXBhcnR1cmUgYW5kIFBvdGVuY3kgRXN0aW1hdGVzIGluIENhcmNpbm9nZW5pYyBEb3NlIFJlc3Bv

bnNlPC90aXRsZT48L3RpdGxlcz48bnVtYmVyPkNPQy9HIDA1PC9udW1iZXI+PGtleXdvcmRzPjxr

ZXl3b3JkPnJpc2sgYW5hbHlzaXM8L2tleXdvcmQ+PGtleXdvcmQ+UVJBPC9rZXl3b3JkPjxrZXl3

b3JkPkJNRDwva2V5d29yZD48a2V5d29yZD5CTURMPC9rZXl3b3JkPjxrZXl3b3JkPm9uZS1zaWRl

ZCBjb25maWRlbmNlIGxpbWl0PC9rZXl3b3JkPjxrZXl3b3JkPmRvc2UgcmVzcG9uc2U8L2tleXdv

cmQ+PGtleXdvcmQ+YmVuY2htYXJrIGFuYWx5c2lzPC9rZXl3b3JkPjxrZXl3b3JkPm5vbmNhbmNl

ciBlbmRwb2ludHM8L2tleXdvcmQ+PGtleXdvcmQ+bm9ucXVhbnRhbCBkYXRhPC9rZXl3b3JkPjxr

ZXl3b3JkPmNhcmNpbm9nZW5lc2lzPC9rZXl3b3JkPjxrZXl3b3JkPnF1YW50YWwgZGF0YTwva2V5

d29yZD48a2V5d29yZD5ub25xdWFudGFsIGRhdGE8L2tleXdvcmQ+PGtleXdvcmQ+UE9EPC9rZXl3

b3JkPjxrZXl3b3JkPlZTRDwva2V5d29yZD48L2tleXdvcmRzPjxkYXRlcz48eWVhcj4yMDE0PC95

ZWFyPjwvZGF0ZXM+PHB1Yi1sb2NhdGlvbj5Mb25kb248L3B1Yi1sb2NhdGlvbj48cHVibGlzaGVy

PlB1YmxpYyBIZWFsdGggRW5nbGFuZCwgQ29tbWl0dGVlIG9uIENhcmNpbm9nZW5pY2l0eSBvZiBD

aGVtaWNhbHMgaW4gRm9vZCwgQ29uc3VtZXIgUHJvZHVjdHMgYW5kIHRoZSBFbnZpcm9ubWVudCA8

L3B1Ymxpc2hlcj48aXNibj5DT0MvRyAwNTwvaXNibj48bGFiZWw+ICBwZGY8L2xhYmVsPjx3b3Jr

LXR5cGU+VGVjaG5pY2FsIFJlcG9ydDwvd29yay10eXBlPjx1cmxzPjxyZWxhdGVkLXVybHM+PHVy

bD5odHRwczovL3d3dy5nb3YudWsvZ292ZXJubWVudC9wdWJsaWNhdGlvbnMvY2FyY2lub2dlbmlj

LWRvc2UtcmVzcG9uc2UtZGVmaW5pbmctYS1wb2ludC1vZi1kZXBhcnR1cmUtYW5kLXBvdGVuY3kt

ZXN0aW1hdGVzPC91cmw+PC9yZWxhdGVkLXVybHM+PC91cmxzPjwvcmVjb3JkPjwvQ2l0ZT48L0Vu

ZE5vdGU+

ADDIN EN.CITE PEVuZE5vdGU+PENpdGU+PEF1dGhvcj5JemFkaTwvQXV0aG9yPjxZZWFyPjIwMTI8L1llYXI+PFJl

Y051bT4zNjE3PC9SZWNOdW0+PERpc3BsYXlUZXh0PihJemFkaTxzdHlsZSBmYWNlPSJpdGFsaWMi

PiBldCBhbC48L3N0eWxlPiwgMjAxMjsgVS5LLiBDb21taXR0ZWUgb24gQ2FyY2lub2dlbmljaXR5

IChDT0MpLCAyMDE0KTwvRGlzcGxheVRleHQ+PHJlY29yZD48cmVjLW51bWJlcj4zNjE3PC9yZWMt

bnVtYmVyPjxmb3JlaWduLWtleXM+PGtleSBhcHA9IkVOIiBkYi1pZD0idGR2ZDB4MmYwcnQyZDFl

MjJybXhldzBwd3Rkd2FldjB3NXN4Ij4zNjE3PC9rZXk+PC9mb3JlaWduLWtleXM+PHJlZi10eXBl

IG5hbWU9IkpvdXJuYWwgQXJ0aWNsZSI+MTc8L3JlZi10eXBlPjxjb250cmlidXRvcnM+PGF1dGhv

cnM+PGF1dGhvcj5JemFkaSwgSG9kYTwvYXV0aG9yPjxhdXRob3I+R3J1bmR5LCBKZWFuIEU8L2F1

dGhvcj48YXV0aG9yPkJvc2UsIFJhbmphbjwvYXV0aG9yPjwvYXV0aG9ycz48L2NvbnRyaWJ1dG9y

cz48dGl0bGVzPjx0aXRsZT5FdmFsdWF0aW9uIG9mIHRoZSBiZW5jaG1hcmsgZG9zZSBmb3IgcG9p

bnQgb2YgZGVwYXJ0dXJlIGRldGVybWluYXRpb24gZm9yIGEgdmFyaWV0eSBvZiBjaGVtaWNhbCBj

bGFzc2VzIGluIGFwcGxpZWQgcmVndWxhdG9yeSBzZXR0aW5nczwvdGl0bGU+PHNlY29uZGFyeS10

aXRsZT5SaXNrIEFuYWx5Ljwvc2Vjb25kYXJ5LXRpdGxlPjxhbHQtdGl0bGU+UmlzayBBbmFseXNp

czwvYWx0LXRpdGxlPjwvdGl0bGVzPjxhbHQtcGVyaW9kaWNhbD48ZnVsbC10aXRsZT5SaXNrIEFu

YWwuPC9mdWxsLXRpdGxlPjxhYmJyLTE+UmlzayBBbmFseXNpczwvYWJici0xPjwvYWx0LXBlcmlv

ZGljYWw+PHBhZ2VzPjgzMC04MzU8L3BhZ2VzPjx2b2x1bWU+MzI8L3ZvbHVtZT48bnVtYmVyPjU8

L251bWJlcj48a2V5d29yZHM+PGtleXdvcmQ+YmVuY2htYXJrIGFuYWx5c2lzPC9rZXl3b3JkPjxr

ZXl3b3JkPkJNRDwva2V5d29yZD48a2V5d29yZD5CTUM8L2tleXdvcmQ+PGtleXdvcmQ+cG9pbnQg

b2YgZGVwYXJ0dXJlPC9rZXl3b3JkPjxrZXl3b3JkPlBPRDwva2V5d29yZD48a2V5d29yZD5yZWd1

bGF0b3J5IGFwcGxpY2F0aW9uczwva2V5d29yZD48L2tleXdvcmRzPjxkYXRlcz48eWVhcj4yMDEy

PC95ZWFyPjwvZGF0ZXM+PGxhYmVsPiAgcGRmPC9sYWJlbD48dXJscz48cmVsYXRlZC11cmxzPjx1

cmw+aHR0cDovL2R4LmRvaS5vcmcvMTAuMTExMS9qLjE1MzktNjkyNC4yMDExLjAxNzMyLng8L3Vy

bD48L3JlbGF0ZWQtdXJscz48L3VybHM+PC9yZWNvcmQ+PC9DaXRlPjxDaXRlPjxBdXRob3I+VS5L

LiBDb21taXR0ZWUgb24gQ2FyY2lub2dlbmljaXR5IChDT0MpPC9BdXRob3I+PFllYXI+MjAxNDwv

WWVhcj48UmVjTnVtPjEyMDQ2PC9SZWNOdW0+PHJlY29yZD48cmVjLW51bWJlcj4xMjA0NjwvcmVj

LW51bWJlcj48Zm9yZWlnbi1rZXlzPjxrZXkgYXBwPSJFTiIgZGItaWQ9InRkdmQweDJmMHJ0MmQx

ZTIycm14ZXcwcHd0ZHdhZXYwdzVzeCI+MTIwNDY8L2tleT48L2ZvcmVpZ24ta2V5cz48cmVmLXR5

cGUgbmFtZT0iUmVwb3J0Ij4yNzwvcmVmLXR5cGU+PGNvbnRyaWJ1dG9ycz48YXV0aG9ycz48YXV0

aG9yPlUuSy4gQ29tbWl0dGVlIG9uIENhcmNpbm9nZW5pY2l0eSAoQ09DKSw8L2F1dGhvcj48L2F1

dGhvcnM+PC9jb250cmlidXRvcnM+PHRpdGxlcz48dGl0bGU+RGVmaW5pbmcgYSBQb2ludCBvZiBE

ZXBhcnR1cmUgYW5kIFBvdGVuY3kgRXN0aW1hdGVzIGluIENhcmNpbm9nZW5pYyBEb3NlIFJlc3Bv

bnNlPC90aXRsZT48L3RpdGxlcz48bnVtYmVyPkNPQy9HIDA1PC9udW1iZXI+PGtleXdvcmRzPjxr

ZXl3b3JkPnJpc2sgYW5hbHlzaXM8L2tleXdvcmQ+PGtleXdvcmQ+UVJBPC9rZXl3b3JkPjxrZXl3

b3JkPkJNRDwva2V5d29yZD48a2V5d29yZD5CTURMPC9rZXl3b3JkPjxrZXl3b3JkPm9uZS1zaWRl

ZCBjb25maWRlbmNlIGxpbWl0PC9rZXl3b3JkPjxrZXl3b3JkPmRvc2UgcmVzcG9uc2U8L2tleXdv

cmQ+PGtleXdvcmQ+YmVuY2htYXJrIGFuYWx5c2lzPC9rZXl3b3JkPjxrZXl3b3JkPm5vbmNhbmNl

ciBlbmRwb2ludHM8L2tleXdvcmQ+PGtleXdvcmQ+bm9ucXVhbnRhbCBkYXRhPC9rZXl3b3JkPjxr

ZXl3b3JkPmNhcmNpbm9nZW5lc2lzPC9rZXl3b3JkPjxrZXl3b3JkPnF1YW50YWwgZGF0YTwva2V5

d29yZD48a2V5d29yZD5ub25xdWFudGFsIGRhdGE8L2tleXdvcmQ+PGtleXdvcmQ+UE9EPC9rZXl3

b3JkPjxrZXl3b3JkPlZTRDwva2V5d29yZD48L2tleXdvcmRzPjxkYXRlcz48eWVhcj4yMDE0PC95

ZWFyPjwvZGF0ZXM+PHB1Yi1sb2NhdGlvbj5Mb25kb248L3B1Yi1sb2NhdGlvbj48cHVibGlzaGVy

PlB1YmxpYyBIZWFsdGggRW5nbGFuZCwgQ29tbWl0dGVlIG9uIENhcmNpbm9nZW5pY2l0eSBvZiBD

aGVtaWNhbHMgaW4gRm9vZCwgQ29uc3VtZXIgUHJvZHVjdHMgYW5kIHRoZSBFbnZpcm9ubWVudCA8

L3B1Ymxpc2hlcj48aXNibj5DT0MvRyAwNTwvaXNibj48bGFiZWw+ICBwZGY8L2xhYmVsPjx3b3Jr

LXR5cGU+VGVjaG5pY2FsIFJlcG9ydDwvd29yay10eXBlPjx1cmxzPjxyZWxhdGVkLXVybHM+PHVy

bD5odHRwczovL3d3dy5nb3YudWsvZ292ZXJubWVudC9wdWJsaWNhdGlvbnMvY2FyY2lub2dlbmlj

LWRvc2UtcmVzcG9uc2UtZGVmaW5pbmctYS1wb2ludC1vZi1kZXBhcnR1cmUtYW5kLXBvdGVuY3kt

ZXN0aW1hdGVzPC91cmw+PC9yZWxhdGVkLXVybHM+PC91cmxzPjwvcmVjb3JkPjwvQ2l0ZT48L0Vu

ZE5vdGU+

ADDIN EN.CITE.DATA (Izadi et al., 2012; U.K. Committee on Carcinogenicity (COC), 2014). We do this using the multistage model as implemented in the U.S. EPA Benchmark Dose Software (BMDS) suite ADDIN EN.CITE <EndNote><Cite><Author>U.S. EPA</Author><Year>2016</Year><RecNum>10402</RecNum><DisplayText>(U.S. EPA, 2016)</DisplayText><record><rec-number>10402</rec-number><foreign-keys><key app="EN" db-id="tdvd0x2f0rt2d1e22rmxew0pwtdwaev0w5sx">10402</key></foreign-keys><ref-type name="Report">27</ref-type><contributors><authors><author>U.S. EPA,</author></authors></contributors><titles><title>Benchmark Dose Software (BMDS) User Manual Version 2.6.0.1</title></titles><keywords><keyword>risk analysis</keyword><keyword>QRA</keyword><keyword>BMD</keyword><keyword>BMDL</keyword><keyword>BMDS software</keyword><keyword>dose response</keyword><keyword>benchmark analysis</keyword><keyword>quantal data</keyword><keyword>Weibull model</keyword><keyword>point of departure</keyword><keyword>gamma model</keyword><keyword>nonquantal data</keyword><keyword>exponential model</keyword><keyword>Hill model</keyword><keyword>polynomial model</keyword><keyword>power model</keyword></keywords><dates><year>2016</year></dates><pub-location>Research Triangle Park, NC</pub-location><publisher>National Center for Environmental Assessment, U.S. Environmental Protection Agency</publisher><label> pdf</label><work-type>Technical Report</work-type><urls><related-urls><url>;(U.S. EPA, 2016). Expanding the multistage response function from main document Equation (2), we study multistage models for the response probability π(u) containing a first-degree, second-degree, or third-degree polynomial. That is, π(u) = ? + (1 – ?)(1 – exp{–?1u – ??u2 – ?3u3}),where γ is the background response rate and the coefficients ?1, ??, and ?3 represent terms for the first-, second-, and third-degree components of the polynomial predictor, respectively. For this analysis, the first-degree multistage model sets ?? = ?? = 0, while the second-degree model sets only ??= 0. For all model fits the coefficients are constrained to be greater than or equal to zero, which for u ≥ 0 imposes the necessary restriction that 0 ≤ π(u) ≤ 1. It also forces π(u) to be monotone increasing.As the parameters are constrained, it is possible for higher order models to degenerate into other models and one may question the necessity of using higher order terms for regulatory purposes. To investigate this, an R script is given below to compare all three models. The following analysis and corresponding R code assumes that the files multistage.exe, bmd.dll, libgcc_s_dw2-1.dll, libgfortran-3.dll and libquadmath-0.dll have been copied from the BMDS install directory and placed in the directory from which this R script is executed. The user-defined function BMDCompute, which is employed below, calls the multistage model and is given in the next supplement Section. The BMDs are computed by setting the extra risk under the multistage model, RE(u) = 1 – exp{–?1u – ?2u2 – ?3u3}, equal to 10%.# fit multistage models with cubic, quadratic or linear componentsBMD.MULTISTAGE3 <-final.data %>% group_by(ID) ?%>%???????????????? select(chemical,r.dose,n,obs) %>%???????????????? do(BMDCompute(.,3))#The unique function is used to remove an extra BMD value that may be #returned by the procedureBMD.MULTISTAGE3<-unique(BMD.MULTISTAGE3)BMD.MULTISTAGE2 <-final.data %>% group_by(ID) ?%>%??????????????????????? select(chemical,r.dose,n,obs) %>%??????????????????????? do(BMDCompute(.,2))BMD.MULTISTAGE2<-unique(BMD.MULTISTAGE2)BMD.MULTISTAGE1 <-final.data %>% group_by(ID) ?%>%???????????????????? ???? select(chemical,r.dose,n,obs) %>%???????????????????? do(BMDCompute(.,1)) ?BMD.MULTISTAGE3<-unique(BMD.MULTISTAGE3)# collect the BMD estimates from the 3 multistage models and # ??calculate their correlationsV = cov(cbind(BMD.MULTISTAGE1[,3],???????????? BMD.MULTISTAGE2[,3],???????????? BMD.MULTISTAGE3[,3]),use="plete.obs")cov2cor(V)This yields the estimated BMD correlation matrix BMD BMD BMDBMD 1.0000000 0.8534608 0.8132934BMD 0.8534608 1.0000000 0.9514441BMD 0.8132934 0.9514441 1.0000000Next, the R code # collect the BMDLs from the 3 multistage models and # ??calculate their correlationsV = cov(cbind(BMD.MULTISTAGE1[,4],???????????? BMD.MULTISTAGE2[,4],???????????? BMD.MULTISTAGE3[,4]),use="plete.obs")cov2cor(V)yields the BMDL correlation matrix BMDL BMDL BMDLBMDL 1.0000000 0.9116890 0.8690537BMDL 0.9116890 1.0000000 0.9607213BMDL 0.8690537 0.9607213 1.0000000We see that the BMD estimates from the database exhibit rather high pairwise correlation coefficients: e.g., the correlation between the first-order model’s BMD estimates and those from the second-order model is found to be 0.85. Similarly, the correlation between the BMDLs from those two models is even higher, at 0.91. Further exploration would show that for these data sets, 40% of the first-order and second-order BMDLs are identical, where the quadratic coefficient is being estimated as zero for the second-order multistage model. Roughly similar values are evidenced for the other pairwise correlations. For example, with the third-order model there is even less divergence from the lower, second-order model: both the BMD estimates and the BMDLs show increases in their correlations, to 0.95 and 0.96, respectively.This analysis shows that there is a complex interrelationship between the three model forms, and that a multistage model higher than second degree may not be necessary. This corresponds to precious indications that at most a second-degree multistage model is adequate for dichotomous dose-response relationships ADDIN EN.CITE <EndNote><Cite><Author>Nitcheva</Author><Year>2007</Year><RecNum>11301</RecNum><DisplayText>(Nitcheva<style face="italic"> et al.</style>, 2007)</DisplayText><record><rec-number>11301</rec-number><foreign-keys><key app="EN" db-id="tdvd0x2f0rt2d1e22rmxew0pwtdwaev0w5sx">11301</key></foreign-keys><ref-type name="Journal Article">17</ref-type><contributors><authors><author>Nitcheva, Danila K</author><author>Piegorsch, Walter W</author><author>West, R Webser</author></authors></contributors><titles><title>On use of the multistage dose-response model for assessing laboratory animal carcinogenicity</title><secondary-title>Regul. Toxicol. Pharmacol.</secondary-title><alt-title>Regulatory Toxicology and Pharmacology</alt-title></titles><periodical><full-title>Regul. Toxicol. Pharmacol.</full-title><abbr-1>Regulatory Toxicology and Pharmacology</abbr-1></periodical><alt-periodical><full-title>Regul. Toxicol. Pharmacol.</full-title><abbr-1>Regulatory Toxicology and Pharmacology</abbr-1></alt-periodical><pages>135-147</pages><volume>48</volume><number>2</number><keywords><keyword>bootstrap hypothesis testing</keyword><keyword>cancer</keyword><keyword>dose-response modeling</keyword><keyword>multistage model</keyword><keyword>quantal response</keyword></keywords><dates><year>2007</year></dates><label>J97&#xD; pdf</label><urls><related-urls><url>;(Nitcheva et al., 2007). It requires the file `multistage.exe,’ which is available with BMDS 2.7. ADDIN S.5.R BMDCompute FUNCTION# the following function performs the following operations...# Call the EPA BMDS software to fit multistage models.# Write the command file that is submitted to BMDS# Input:# data := data frame with the following structure# deg := degree of the multistage model that was fit# Output:# dataframe containing two columns: BMD, BMDLlibrary (dplyr)BMDCompute<- function(data,deg){ print(data) # Check if there are more than two data points # also check if there is a non 0% or 100% response. if ((nrow(data) > 2) * (sum((data[,5]/data[,4] < 1) * ( data[,5]/data[,4] > 0 )) > 0) ){ X = c( 'Multistage', 'BMDS_Model_Run', 'NoDataFile','.\\run.out', sprintf('%d %d',nrow(data),deg), '500 1.00E-08 1.00E-08 0 1 1 0 0', '0.1000 0 0.95', paste(rep('-9999 ',deg+1),collapse=""), '0', paste(rep('-9999 ',deg+1),collapse=""), 'Dose Effect NEGATIVE_RESPONSE') for (i in 1:nrow(data)){ X = c(X,sprintf('%f %d %d',data[i,3],as.numeric(data[i,5]),as.numeric(data[i,4]-data[i,5]))) } write(X,"out.(d)",sep="\n") system("multistage out.(d)") infile <- file("out.out", "rt") file.lines=readLines(infile) if (length(grep("BMD computation failed",file.lines)) > 0){ close(infile) return(data.frame(chemname = data[1,2],BMD = NA, BMDL=NA)) }else{ BMD = as.numeric(trimws(sub("BMD =*","",file.lines[grep("BMD =*",file.lines)]), "both")) BMDL = as.numeric(trimws(sub("BMDL =*","",file.lines[grep("BMDL =*",file.lines)]),"both")) close(infile) return(data.frame(chemname = data[1,2],BMD=BMD,BMDL=BMDL)) } } else { return(data.frame(chemname = data[1,2],BMD = NA, BMDL=NA)) }} #end of function ADDIN S.6.Empirical distributions and recommended priors for all models in table 1To complement the graphical summaries in main document Figs. 2–4 of the empirical distributions seen in QRAD, we present here similar graphic displays for the remaining model/parameter combinations in main document Table 1. For easy referral, these are as follows:(1)Quantal-linear ??→ see Fig. S.4Quantal-linear ?1?→ see Fig. S.5(2)Multi-stage ??→ see Fig. S.6Multi-stage ?1?→ see Fig. S.7Multi-stage ?2?→ see Fig. S.8(3)Weibull ??→ see Fig. S.9Weibull ??→ see main document Fig. 4Weibull ??→ see Fig. S.10(4)Gamma ??→ see Fig. S.11Gamma ??→ see Fig. S.12Gamma ??→ see Fig. S.13(5)Logistic ?0?→ see Fig. S.14Logistic ?1?→ see Fig. S.15(6)Log-logistic ??→ see Fig. S.16Log-logistic ?0?→ see main document Fig. 3Log-logistic ?1?→ see Fig. S.17(7)Probit ?0?→ see Fig. S.18Probit ?1?→ see main document Fig. 2(8)Log-probit ??→ see Fig. S.19Log-probit ?0?→ see Fig. S.20Log-probit ?1?→ see Fig. S.21Fig. S.4: Empirical distribution (histogram) relating to background response parameter ? for the Quantal-linear dose-response model (1) and corresponding recommended prior (dark curve), ? ~ Beta(0.313, 2.545), from Table 1. Fig. S.5: Empirical distribution (histogram) relating to ?1 parameter for the Quantal-linear dose-response model (1) and corresponding recommended prior (dark curve), ?1 ~ LN(–0.470, 1.688), from Table 1. Fig. S.6: Empirical distribution (histogram) relating to background response parameter ? for the Multi-stage dose-response model (2) and corresponding recommended prior (dark curve), ? ~ Beta(0.326, 2.582), from Table 1. Fig. S.7: Empirical distribution (histogram) relating to ?1 parameter for the Multi-stage dose-response model (2) and corresponding recommended prior (dark curve), ?1 ~ LN(–1.006, 4.949), from Table 1. Fig. S.8: Empirical distribution (histogram) relating to ?2 parameter for the Multi-stage dose-response model (2) and corresponding recommended prior (dark curve), ?2 ~ LN(–1.016, 5.013), from Table 1. Fig. S.9: Empirical distribution (histogram) relating to background response parameter ? for the Weibull dose-response model (3) and corresponding recommended prior (dark curve), ? ~ Beta(0.271, 2.583), from Table 1. Fig. S.10: Empirical distribution (histogram) relating to ? parameter for the Weibull dose-response model (3) and corresponding recommended prior (dark curve), ? ~ LN(–0.464, 1.535), from Table 1. Fig. S.11: Empirical distribution (histogram) of background response parameter ? for the gamma dose-response model (4) and corresponding recommended prior (dark curve), ? ~ Beta(0.276, 2.572), from Table 1. Fig. S.12: Empirical distribution (histogram) of shape parameter α for the gamma dose-response model (4) and corresponding recommended prior (dark curve), based on ? = ?–0.2 ~ LN(–0.153, 2.637), from Table 1. Fig. S.13: Empirical distribution (histogram) relating to ? parameter for the gamma dose-response model (4) and corresponding recommended prior (dark curve), ? ~ LN(–0.587, 5.659), from Table 1. Fig. S.14. Empirical distribution (histogram) relating to parameter β0 for the logistic dose-response model (5) and corresponding recommended prior (dark curve), β0 ~ N(–2.526, 3.733), from Table 1.Fig. S.15. Empirical distribution (histogram) relating to parameter β1 for the logistic dose-response model (5) and corresponding recommended prior (dark curve), β1 ~ LN(1.018, 0.603), from Table 1.Fig. S.16: Empirical distribution (histogram) of background response parameter ? for the log-logistic dose-response model (6) and corresponding recommended prior (dark curve), ? ~ Beta(0.275, 2.571), from Table 1. Fig. S.17. Empirical distribution (histogram) relating to parameter β1 for the log-logistic dose-response model (6) and corresponding recommended prior (dark curve), β1 ~ LN(0.274, 0.960), from Table 1.Fig. S.18. Empirical distribution (histogram) relating to parameter β0 for the probit dose-response model (7) and corresponding recommended prior (dark curve), β0 ~ N(–1.66, 3.264), from Table 1.Fig. S.19. Empirical distribution (histogram) of background response parameter ? for the log-probit dose-response model (8) and corresponding recommended prior (dark curve), ? ~ Beta(0.333, 2.877), from Table 1.Fig. S.20. Empirical distribution (histogram) relating to parameter β0 for the log-probit dose-response model (8) and corresponding recommended prior (dark curve), β0 ~ N(–0.514,10.337), from Table 1.Fig. S.21. Empirical distribution (histogram) relating to parameter β1 for the log-probit dose-response model (8) and corresponding recommended prior (dark curve), β1 ~ LN(–0.186, 0.579), from Table 1. ADDIN S.7.R Code to for Developing PriorsThe following code uses the database to develop the priors. library(MABMDDichotomous)library(nloptr)library(dplyr)load(url(""))#Weibull Dose-Responseweibull <- function(BETAS,n,o,d){ gamma <- BETAS[1] alpha <- BETAS[3] beta <- BETAS[2] p <- gamma + (1-gamma)*(1-exp(-beta*d^alpha)) p[d==0] = gamma return(p)}#Binomial likelihoodweibbinLike <- function(par,n,o,d){ p <- weibull(par,n,o,d) log.like = o*log(p) + (n-o)*log(1-p) log.like[p==0] = log(1e-8) log.like[p==1] = log(1e-8) return(-sum(log.like))}#Weibull Dose-Responsegamma.D <- function(BETAS,n,o,d){ gamma <- BETAS[1] alpha <- BETAS[3] beta <- BETAS[2] p <- gamma + (1-gamma)*pgamma(beta*d,alpha) p[d==0] = gamma return(p)}#Binomial likelihoodgambinLike <- function(par,n,o,d){ p <- gamma.D(par,n,o,d) log.like = o*log(p) + (n-o)*log(1-p) log.like[p==0] = log(1e-8) log.like[p==1] = log(1e-8) return(-sum(log.like))}BMDComputeWeibull<- function(data){ rval3<-nloptr(c(0.1,1,1),eval_f=weibbinLike, lb=c(0,0,0.2), ub=c(1,18,18), n=data$n,o=data$obs,d=data$r.dose, opts = list("algorithm"="NLOPT_LN_BOBYQA",xtol_rel = 1e-6,maxeval=3000)) return(data.frame(ID = data$ID[1],gamma=rval3$solution[1],beta=rval3$solution[2],alpha=rval3$solution[3],AIC = 2*rval3$objective+2*3)) }BMDComputeGamma<- function(data){ rval3<-nloptr(c(0.1,1,1),eval_f=gambinLike, lb=c(0,0,0.2), ub=c(1,18,18), n=data$n,o=data$obs,d=data$r.dose, opts = list("algorithm"="NLOPT_LN_BOBYQA",xtol_rel = 1e-6,maxeval=3000)) return(data.frame(ID = data$ID[1],gamma=rval3$solution[1],beta=rval3$solution[2],alpha=rval3$solution[3],AIC = 2*rval3$objective+2*3)) }####################################################################333########################################################################log-logisticlog.logistic <- function(BETAS,n,o,d){ gamma <- BETAS[1] beta0 <- BETAS[2] beta1 <- BETAS[3] p <- gamma + (1-gamma)*(1/(1+exp(-beta0-beta1*log(d)))) p[d==0] = gamma return(p)}loglogitlike <- function(par,n,o,d){ p <- log.logistic(par,n,o,d) log.like = o*log(p) + (n-o)*log(1-p) log.like[p==0] = log(1e-8) log.like[p==1] = log(1e-8) return(-sum(log.like))}BMDComputeLlogistic <- function(data){ rval3<-nloptr(c(0.1,1,1),eval_f=loglogitlike, lb=c(0,-18,0.2), ub=c(1,18,18), n=data$n,o=data$obs,d=data$r.dose, opts = list("algorithm"="NLOPT_LN_BOBYQA",xtol_rel = 1e-6,maxeval=3000)) return(data.frame(ID = data$ID[1],gamma=rval3$solution[1],beta0=rval3$solution[2],beta1=rval3$solution[3],AIC = 2*rval3$objective+2*3)) }#####################################################################################################################log-probitlog.probit <- function(BETAS,n,o,d){ gamma <- BETAS[1] beta0 <- BETAS[2] beta1 <- BETAS[3] p <- gamma + (1-gamma)*pnorm(beta0+beta1*log(d)) p[d==0] = gamma return(p)}logprobitlike <- function(par,n,o,d){ p <- log.probit(par,n,o,d) log.like = o*log(p) + (n-o)*log(1-p) log.like[p==0] = log(1e-8) log.like[p==1] = log(1e-8) return(-sum(log.like))}BMDComputeLprobit <- function(data){ rval3<-nloptr(c(0.1,0,1),eval_f=logprobitlike, lb=c(0,-18,0.2), ub=c(1,18,18), n=data$n,o=data$obs,d=data$r.dose, opts = list("algorithm"="NLOPT_LN_BOBYQA",xtol_rel = 1e-6,maxeval=3000)) return(data.frame(ID = data$ID[1],gamma=rval3$solution[1],beta0=rval3$solution[2],beta1=rval3$solution[3],AIC = 2*rval3$objective+2*3)) }###########################################################################Quantal Linear and Multistage###########################################################################Binomial likelihood Logistic#log-probitmultistage <- function(BETAS,n,o,d){ gamma <- BETAS[1] beta1 <- BETAS[2] beta2 <- BETAS[3] p <- gamma + (1-gamma)*(1-exp(-beta1*d-beta2*d^2)) return(p)}qlinear <- function(BETAS,n,o,d){ gamma <- BETAS[1] beta1 <- BETAS[2] p <- gamma + (1-gamma)*(1-exp(-beta1*d)) return(p)}multistagelike <- function(par,n,o,d){ p <- multistage(par,n,o,d) log.like = o*log(p) + (n-o)*log(1-p) log.like[p==0] = log(1e-8) log.like[p==1] = log(1e-8) return(-sum(log.like))}qlinearlike <- function(par,n,o,d){ p <- qlinear(par,n,o,d) log.like = o*log(p) + (n-o)*log(1-p) log.like[p==0] = log(1e-8) log.like[p==1] = log(1e-8) return(-sum(log.like))}BMDComputeMstage <- function(data){ rval3<-nloptr(c(0.1,1,1),eval_f=multistagelike, lb=c(0,0,0), ub=c(1,18,18), n=data$n,o=data$obs,d=data$r.dose, opts = list("algorithm"="NLOPT_LN_BOBYQA",xtol_rel = 1e-6,maxeval=3000)) return(data.frame(ID = data$ID[1],gamma=rval3$solution[1],beta1=rval3$solution[2],beta1=rval3$solution[3],AIC = 2*rval3$objective+2*3)) }BMDComputeQlinear<- function(data){ rval3<-nloptr(c(0.1,1),eval_f=qlinearlike, lb=c(0,0), ub=c(1,18), n=data$n,o=data$obs,d=data$r.dose, opts = list("algorithm"="NLOPT_LN_BOBYQA",xtol_rel = 1e-6,maxeval=3000)) return(data.frame(ID = data$ID[1],gamma=rval3$solution[1],beta1=rval3$solution[2],AIC = 2*rval3$objective+2*2)) }###########################################################################Quantal Linear and Multistage###########################################################################Binomial likelihood Logistic#log-probitprobit <- function(BETAS,n,o,d){ beta0 <- BETAS[1] beta1 <- BETAS[2] p <- pnorm(beta0+beta1*d) return(p)}logistic <- function(BETAS,n,o,d){ beta0 <- BETAS[1] beta1 <- BETAS[2] p <- 1/(1+ exp(-beta0 -beta1*d)) return(p)}probitlike <- function(par,n,o,d){ p <- probit(par,n,o,d) log.like = o*log(p) + (n-o)*log(1-p) log.like[p==0] = log(1e-8) log.like[p==1] = log(1e-8) return(-sum(log.like))}logisticlike <- function(par,n,o,d){ p <- logistic(par,n,o,d) log.like = o*log(p) + (n-o)*log(1-p) log.like[p==0] = log(1e-8) log.like[p==1] = log(1e-8) return(-sum(log.like))}BMDComputelogistic <- function(data){ rval3<-nloptr(c(0,0),eval_f=logisticlike, lb=c(-18,0), ub=c(18,18), n=data$n,o=data$obs,d=data$r.dose, opts = list("algorithm"="NLOPT_LN_BOBYQA",xtol_rel = 1e-6,maxeval=3000)) return(data.frame(ID = data$ID[1],beta0=rval3$solution[1],beta1=rval3$solution[2],AIC = 2*rval3$objective+2*2))}BMDComputeprobit<- function(data){ rval3<-nloptr(c(-1.2,0),eval_f=probitlike, lb=c(-18,0), ub=c(18,18), n=data$n,o=data$obs,d=data$r.dose, opts = list("algorithm"="NLOPT_LN_BOBYQA",xtol_rel = 1e-6,maxeval=3000)) return(data.frame(ID = data$ID[1],beta0=rval3$solution[1],beta1=rval3$solution[2],AIC = 2*rval3$objective+2*2)) }results.probit <- final.data %>% group_by(ID) %>% select(ID,r.dose,n,obs) %>% filter(length(r.dose) >= 4) %>% filter(sum(n) > 100) %>% do(BMDComputeprobit(.))results.logistic <- final.data %>% group_by(ID) %>% select(ID,r.dose,n,obs) %>% filter(length(r.dose) >= 4) %>% filter(sum(n) > 100) %>% do(BMDComputelogistic(.))results.llogistic <- final.data %>% group_by(ID) %>% select(ID,r.dose,n,obs) %>% filter(length(r.dose) >= 4) %>% filter(sum(n) > 100) %>% do(BMDComputeLlogistic(.))results.multistage <- final.data %>% group_by(ID) %>% select(ID,r.dose,n,obs) %>% filter(length(r.dose) >= 4) %>% filter(sum(n) > 100) %>% do(BMDComputeMstage(.))results.Qlinear <- final.data %>% group_by(ID) %>% select(ID,r.dose,n,obs) %>% filter(length(r.dose) >= 4) %>% filter(sum(n) > 100) %>% do(BMDComputeQlinear(.))results.lprobit <- final.data %>% group_by(ID) %>% select(ID,r.dose,n,obs) %>% filter(length(r.dose) >= 4) %>% filter(sum(n) > 100) %>% do(BMDComputeLprobit(.))results.weibull <- final.data %>% group_by(ID) %>% select(ID,r.dose,n,obs) %>% filter(length(r.dose) >= 4) %>% filter(sum(n) > 100) %>% do(BMDComputeWeibull(.))results.gamma <- final.data %>% group_by(ID) %>% select(ID,r.dose,n,obs) %>% filter(length(r.dose) >= 4) %>% filter(sum(n) > 100) %>% do(BMDComputeGamma(.))AIC.MAT<-cbind(results.weibull$AIC,results.gamma$AIC,results.llogistic$AIC,results.logistic$AIC, results.probit$AIC, results.lprobit$AIC,results.multistage$AIC,results.Qlinear$AIC)mins <- apply(AIC.MAT,1,which.min)model <- matrix(1:8,nrow=nrow(AIC.MAT),ncol=8,byrow=TRUE)priorWeights3 <- colMeans(model==mins)names(priorWeights3) <- c("Weibull","Gamma","Log-Logistic","Logistic","Probit","Log-Probit","Multistage","QLinear")###########################################################################################################################3#Create the emperical prior distributions###################################################results.llogisticpdf(file="prior-histograms.pdf")beta0 <- ggplot(results.llogistic, aes(x = beta0)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.4) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + ggtitle("Log-Logistic") + xlab(expression(paste(beta[0]))) + stat_function(fun = dnorm, colour = "black", args = list(mean = mean(results.llogistic$beta0, na.rm = TRUE), sd = sd(results.llogistic$beta0, na.rm = TRUE)),size=1.25)beta0 mean(results.llogistic$beta0)sd(results.llogistic$beta0)beta0 <- ggplot(results.llogistic, aes(x = beta1)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.4) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + xlab(expression(paste(beta[1]))) + ggtitle("Log-Logistic") + scale_x_continuous(limits = c(0, 10)) + stat_function(fun = dlnorm, colour = "black", args = list(mean = mean(log(results.llogistic$beta1), na.rm = TRUE), sd = sd(log(results.llogistic$beta1), na.rm = TRUE)),size=1.25)beta0mean(log(results.llogistic$beta1),na.rm = T)sd(log(results.llogistic$beta1),na.rm = T)####################################################### Log Probit#############################################33beta0 <- ggplot(results.lprobit, aes(x = beta0)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.4) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + ggtitle("Log-Probit") + xlab(expression(paste(beta[0]))) + stat_function(fun = dnorm, colour = "black", args = list(mean = mean(results.lprobit$beta0, na.rm = TRUE), sd = sd(results.lprobit$beta0, na.rm = TRUE)),size=1.25)beta0 mean(results.lprobit$beta0)sd(results.lprobit$beta0)beta0 <- ggplot(results.lprobit, aes(x = beta1)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.4) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + xlab(expression(paste(beta[1]))) + ggtitle("Log-Probit") + scale_x_continuous(limits = c(0, 10)) + stat_function(fun = dlnorm, colour = "black", args = list(mean = mean(log(results.lprobit$beta1), na.rm = TRUE), sd = sd(log(results.lprobit$beta1), na.rm = TRUE)),size=1.25)beta0mean(log(results.lprobit$beta1),na.rm = T)sd(log(results.lprobit$beta1),na.rm = T)#################################################################################33##################################################3#Create the emperical prior distributions###################################################beta0 <- ggplot(results.logistic[(results.logistic$beta0 > -10)*(results.logistic$beta0 < 10)==1,], aes(x = beta0)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.4) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + ggtitle("Logistic") + xlab(expression(paste(beta[0]))) + stat_function(fun = dnorm, colour = "black", args = list(mean = mean(results.logistic$beta0[(results.logistic$beta0 > -10)*(results.logistic$beta0 < 10)==1], na.rm = TRUE), sd = sd(results.logistic$beta0[(results.logistic$beta0 > -10)*(results.logistic$beta0 < 10)==1], na.rm = TRUE)),size=1.25)beta0 mean(results.logistic$beta0[(results.logistic$beta0 > -10)*(results.logistic$beta0 < 10)==1])sd(results.logistic$beta0[(results.logistic$beta0 > -10)*(results.logistic$beta0 < 10)==1])beta0 <- ggplot(results.logistic, aes(x = beta1)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.4) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + xlab(expression(paste(beta[1]))) + ggtitle("Logistic") + scale_x_continuous(limits = c(0, 10)) + stat_function(fun = dlnorm, colour = "black", args = list(mean = mean(log(results.logistic$beta1+0.1), na.rm = TRUE), sd = sd(log(results.logistic$beta1+0.1), na.rm = TRUE)),size=1.25)beta0mean(log(results.logistic$beta1+0.1),na.rm = T)sd(log(results.logistic$beta1+0.1),na.rm = T)####################################################### Log Probit#############################################33beta0 <- ggplot(results.probit[(results.probit$beta0 > -5)*(results.probit$beta0 < 5)==1,], aes(x = beta0)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.4) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + ggtitle("Probit") + xlab(expression(paste(beta[0]))) + stat_function(fun = dnorm, colour = "black", args = list(mean = mean(results.probit$beta0[(results.probit$beta0 > -5)*(results.probit$beta0 < 5)==1], na.rm = TRUE), sd = sd(results.probit$beta0[(results.probit$beta0 > -5)*(results.probit$beta0 < 5)==1], na.rm = TRUE)),size=1.25)beta0 mean(results.probit$beta0[(results.probit$beta0 > -5)*(results.probit$beta0 < 5)==1])sd(results.probit$beta0[(results.probit$beta0 > -5)*(results.probit$beta0 < 5)==1])beta0 <- ggplot(results.probit[results.probit$beta1 >0,], aes(x = beta1)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.4) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + xlab(expression(paste(beta[1]))) + ggtitle("Probit") + scale_x_continuous(limits = c(0, 10)) + stat_function(fun = dlnorm, colour = "black", args = list(mean = mean(log(results.probit$beta1[results.probit$beta1 >0]), na.rm = TRUE), sd = sd(log(results.probit$beta1[results.probit$beta1 >0]), na.rm = TRUE)),size=1.25)beta0mean(log(results.probit$beta1),na.rm = T)sd(log(results.probit$beta1),na.rm = T)######################################################################################################################################## Log Quantal Linear######################################################beta0 <- ggplot(results.Qlinear[results.Qlinear$beta1 > 0,], aes(x = beta1)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.25) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + ggtitle("Quantal-Linear") + xlab(expression(paste(beta[1]))) + stat_function(fun = dlnorm, colour = "black", args = list(mean = mean(log(results.Qlinear$beta1[results.Qlinear$beta1 > 0]), na.rm = TRUE), sd = sd(log(results.Qlinear$beta1[results.Qlinear$beta1 > 0]), na.rm = TRUE)),size=1.25)beta0 mean(log(results.Qlinear$beta1[results.Qlinear$beta1 > 0]),na.rm = T)sd(log(results.Qlinear$beta1[results.Qlinear$beta1 > 0]),na.rm = T)######################################################################################################################################## Multistage######################################################beta0 <- ggplot(results.multistage[results.multistage$beta1 > 0,], aes(x = beta1)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.25) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + ggtitle("Multistage 2-Degree") + xlab(expression(paste(beta[1]))) + stat_function(fun = dlnorm, colour = "black", args = list(mean = mean(log(results.multistage$beta1[results.multistage$beta1 > 0]), na.rm = TRUE), sd = sd(log(results.multistage$beta1[results.multistage$beta1 > 0]), na.rm = TRUE)),size=1.25)beta0 mean(log(results.multistage$beta1[results.multistage$beta1 > 0]),na.rm = T)sd(log(results.multistage$beta1[results.multistage$beta1 > 0]),na.rm = T)beta0 <- ggplot(results.multistage[results.multistage$beta1.1 > 0,], aes(x = beta1.1)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.25) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + ggtitle("Multistage 2-Degree") + xlab(expression(paste(beta[2]))) + stat_function(fun = dlnorm, colour = "black", args = list(mean = mean(log(results.multistage$beta1.1[results.multistage$beta1.1 > 0]), na.rm = TRUE), sd = sd(log(results.multistage$beta1.1[results.multistage$beta1.1 > 0]), na.rm = TRUE)),size=1.25)beta0 mean(log(results.multistage$beta1.1[results.multistage$beta1.1 > 0]),na.rm = T)sd(log(results.multistage$beta1.1[results.multistage$beta1.1 > 0]),na.rm = T)####################################################### Gamma######################################################beta0 <- ggplot(results.gamma[results.gamma$beta > 0,], aes(x = beta)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.20) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + ggtitle("Gamma") + xlab(expression(paste(beta))) + stat_function(fun = dlnorm, colour = "black", args = list(mean = mean(log(results.gamma$beta[results.gamma$beta > 0]), na.rm = TRUE), sd = sd(log(results.gamma$beta[results.gamma$beta > 0]), na.rm = TRUE)),size=1.25)beta0 mean(log(results.gamma$beta[results.gamma$beta > 0]),na.rm = T)sd(log(results.gamma$beta[results.gamma$beta > 0]),na.rm = T)beta0 <- ggplot(results.gamma[results.gamma$alpha > 0,], aes(x = alpha)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.30) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + ggtitle("Gamma") + xlab(expression(paste(alpha))) + stat_function(fun = dlnorm, colour = "black", args = list(mean = mean(log(results.gamma$alpha[results.gamma$alpha > 0]), na.rm = TRUE), sd = sd(log(results.gamma$alpha[results.gamma$alpha > 0]), na.rm = TRUE)),size=1.25)beta0 mean(log(results.gamma$alpha[results.gamma$alpha > 0]),na.rm = T)sd(log(results.gamma$alpha[results.gamma$alpha > 0]),na.rm = T)####################################################### weibull######################################################beta0 <- ggplot(results.weibull[results.weibull$beta > 0,], aes(x = beta)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.20) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + ggtitle("Weibull") + xlab(expression(paste(beta))) + stat_function(fun = dlnorm, colour = "black", args = list(mean = mean(log(results.weibull$beta[results.weibull$beta > 0]), na.rm = TRUE), sd = sd(log(results.weibull$beta[results.weibull$beta > 0]), na.rm = TRUE)),size=1.25)beta0 mean(log(results.weibull$beta[results.weibull$beta > 0]),na.rm = T)sd(log(results.weibull$beta[results.weibull$beta > 0]),na.rm = T)beta0 <- ggplot(results.weibull[results.weibull$alpha > 0,], aes(x = alpha)) + geom_histogram(aes(y = ..density..),center = 0,colour="grey20",fill="grey90",binwidth = 0.30) + theme_bw() + theme(panel.border = element_blank(), axis.line = element_line(colour = "black")) + ggtitle("Weibull") + xlab(expression(paste(alpha))) + stat_function(fun = dlnorm, colour = "black", args = list(mean = mean(log(results.weibull$alpha[results.weibull$alpha > 0]), na.rm = TRUE), sd = sd(log(results.weibull$alpha[results.weibull$alpha > 0]), na.rm = TRUE)),size=1.25)beta0 mean(log(results.weibull$alpha[results.weibull$alpha > 0]),na.rm = T)sd(log(results.weibull$alpha[results.weibull$alpha > 0]),na.rm = T)dev.off() ADDIN S.7.R Code for Weibull Bayesian Analysis for 3-monochloropropane-1,2-diolThe following code is the R/RSTAN code for the Bayesian analysis. library(MABMDDichotomous)library(dplyr)library(rstan)N <- c(50,50,50,50)D <- c(0,1.97,8.27,29.5)/29.5O <- c(1,11,21,36)Adata <- list(len=length(N),y=O,n=N,d=D)wei_model.p <- "data { int<lower=0> len; //number of dose groups int<lower=0> y[len]; //observed number of cases int<lower=0> n[len]; //number of subjects real<lower=0> d[len]; //dose levels}parameters { real <lower=0,upper=1> gamma; real <lower = 0> beta; real <lower = 0> alpha;}model { alpha ~ lognormal(-0.243,1.61); //note: STD devations are specified here beta ~ lognormal(-0.587,1.57);// variances are in the manuscript gamma ~ beta(0.271, 2.583); //likelihood for (i in 1:len){ y[i] ~ binomial(n[i],gamma+(1-gamma)*(1-exp(-beta*(d[i])^(0.2+alpha)))); }}"mcmc_iter.wei = 25000mcmc_warmup.wei = 1000fit_wei = stan(model_code=wei_model.p,data=Adata,iter=mcmc_iter.wei,warmup=mcmc_warmup.wei,chains=1,seed=8675309)list_of_draws <- extract(fit_wei)BMD.P <- (-log(0.9)/list_of_draws$beta)^(1/(list_of_draws$alpha+0.2))*29.5wei_model <- "data {int<lower=0> len; //number of dose groupsint<lower=0> y[len]; //observed number of casesint<lower=0> n[len]; //number of subjectsreal<lower=0> d[len]; //dose levels}parameters {real <lower=0,upper=1> gamma;real <lower = 0> beta; real <lower = 0> alpha;}model {alpha ~ uniform(0.2,18); //note: STD devations are specified herebeta ~ uniform(0,18);// variances are in the manuscriptgamma ~ uniform(0,18);//likelihoodfor (i in 1:len){y[i] ~ binomial(n[i],gamma+(1-gamma)*(1-exp(-beta*(d[i])^(alpha))));}}"mcmc_iter.wei = 25000mcmc_warmup.wei = 1000fit_wei2 = stan(model_code=wei_model,data=Adata,iter=mcmc_iter.wei,warmup=mcmc_warmup.wei,chains=1,seed=8675309)list_of_draws2 <- extract(fit_wei2)BMD <- (-log(0.9)/list_of_draws2$beta)^(1/(list_of_draws2$alpha))*29.5########################################################################################################gamma###################################################N <- c(50,50,50,50)D <- c(0,1.97,8.27,29.5)/29.5O <- c(1,11,21,36)Adata <- list(len=length(N),y=O,n=N,d=D)gamma_model.p <- "data {int<lower=0> len; //number of dose groupsint<lower=0> y[len]; //observed number of casesint<lower=0> n[len]; //number of subjectsreal<lower=0> d[len]; //dose levels} parameters { real <lower=0,upper=1> gamma; real <lower = 0> beta; real <lower = 0> alpha; }model { alpha ~ lognormal(-0.243,1.60); //note: STD devations are specified here beta ~ lognormal(-0.587,1.62);// variances are in the manuscript gamma ~ beta(0.276, 2.572);//likelihoodfor (i in 1:len){if (d[i] > 0){ y[i] ~ binomial(n[i],gamma+(1-gamma)*gamma_cdf(beta*d[i],alpha+0.2,1));}else{ y[i] ~ binomial(n[i],gamma);}}}"mcmc_iter.wei = 25000mcmc_warmup.wei = 1000fit_gamma = stan(model_code=gamma_model.p,data=Adata,iter=mcmc_iter.wei,warmup=mcmc_warmup.wei,chains=1,seed=8675309)list_of_draws <- extract(fit_gamma)GAMMA.BMD.P <- (qgamma(0.1,list_of_draws$alpha+0.2,1)/list_of_draws$beta)*29.5gam_model <- "data {int<lower=0> len; //number of dose groupsint<lower=0> y[len]; //observed number of casesint<lower=0> n[len]; //number of subjectsreal<lower=0> d[len]; //dose levels}parameters {real <lower=0,upper=1> gamma;real <lower = 0> beta; real <lower = 0> alpha;}model {alpha ~ uniform(0.2,18); //note: STD devations are specified herebeta ~ uniform(0,18);// variances are in the manuscriptgamma ~ uniform(0,18);//likelihood for (i in 1:len){ if (d[i] > 0){ y[i] ~ binomial(n[i],gamma+(1-gamma)*gamma_cdf(beta*d[i],alpha,1)); }else{ y[i] ~ binomial(n[i],gamma); } }}"mcmc_iter.wei = 25000mcmc_warmup.wei = 1000fit_gamma2 = stan(model_code=gam_model,data=Adata,iter=mcmc_iter.wei,warmup=mcmc_warmup.wei,chains=1,seed=8675309)list_of_draws2 <- extract(fit_gamma2)GAMMA.BMD <- (qgamma(0.1,list_of_draws2$alpha,1)/list_of_draws2$beta)*29.5References ADDIN EN.REFLIST Crump, K. S. (1984). A new method for determining allowable daily intake. Fundamental and Applied Toxicology 4, 854-871. PMCID/PMID: PMID6510615.Izadi, H., Grundy, J. E., and Bose, R. (2012). Evaluation of the benchmark dose for point of departure determination for a variety of chemical classes in applied regulatory settings. Risk Analysis 32, 830-835.Mack, G. A., and Wolfe, D. A. (1981). K-sample rank tests for umbrella alternatives. Journal of the American Statistical Association 76, 175-181.Nitcheva, D. K., Piegorsch, W. W., and West, R. W. (2007). On use of the multistage dose-response model for assessing laboratory animal carcinogenicity. Regulatory Toxicology and Pharmacology 48, 135-147. PMCID/PMID: PMC2040324.Piegorsch, W. W., and Bailer, A. J. (2005). Analyzing Environmental Data. Chichester: John Wiley & Sons.U.K. Committee on Carcinogenicity (COC) (2014). Defining a Point of Departure and Potency Estimates in Carcinogenic Dose Response. Technical Report number COC/G 05. Public Health England, Committee on Carcinogenicity of Chemicals in Food, Consumer Products and the Environment London.U.S. EPA (2016). Benchmark Dose Software (BMDS) User Manual Version 2.6.0.1. Technical Report. National Center for Environmental Assessment, U.S. Environmental Protection Agency, Research Triangle Park, NC.Wickham, H., and Francois, R. (2016). dplyr: A grammar of data manipulation. R package version 0.5.0. . . ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download