The incorporation of animal and other nonhuman data into epidemiologic analyses is currently a topic of interest to many, including the International Agency for Research on Cancer (IARC) which classifies substantives with regard to their carcinogenicity. IARC convened a workshop on quantitative risk assessment in 2013; this workshop recommended formal consideration of mechanistic data and animal studies in quantitative risk assessment for humans, without outlining how this might be done.1 Bayesian analyses are increasingly used in a wide variety of fields to incorporate prior knowledge and beliefs, and to constrain parameter estimates to be within bounds believed to be reasonable, a priori. Here, we demonstrate the application of Bayesian methods to dose–response analysis using both human and animal studies and adjusting for exposure measurement error in the human studies, using silica exposure and lung cancer as an example. Silica is carcinogenic to humans (group 1) according to the International Agency for Research on Cancer. It is a common occupational exposure in many parts of the world, and the US Occupational Safety and Health Administration has recently lowered the occupational permissible limit in the United States, based on the relatively extensive human exposure–response data. In addition, there are several well-conducted animal studies of silica and lung cancer with multiple doses, providing dose–response information from experiments using controlled doses.
The first main feature of this article is the incorporation of both human and animal data in the same dose–response analysis. Information from experimental animal models provides prior data on the exposure–response between inhaled silica and subsequent development of lung carcinomas, but interspecies extrapolation of animal results to human are fraught with uncertainties due to differences in metabolism across species, typically greater doses to animals than are experienced by humans, different routes of exposure, and fewer concomitant exposures that could affect cancer risk. Despite these limitations, long-term cancer bioassays are predictive for identifying human carcinogens; the overwhelming majority of established human carcinogens have also been shown to cause cancer in animals.2,3 Moreover, regulatory agencies such as the US Environmental Protection Agency (EPA) and US Food and Drug Administration have long used animal cancer data to identify and regulate potential carcinogens and estimate plausible upper bounds on cancer risks for human exposures.4,5
Our analyses also include adjustments for measurement error in the absence of a gold standard error-free measurement, which might be used for calibration. Human exposure–response studies always involve measurement error. In the case of silica measurement, there are two obvious sources of possible error. The first is the assignment of silica concentrations to specific jobs/industries based on a sample of actual measurements of silica among workers in these same specific jobs/industries. The second is the conversion of older measurements of silica in units of millions of dust particles per cubic feet to silica gravimetric units of mg/m3. Conversion from dust particle counts to silica mass depends on the composition and particle size distribution of the dust, which vary across jobs; the conversion factor is based on the average value from a sample of dust measurements. Both of these sources of measurement error are likely to be errors primarily of the Berkson type, in which the true value is equal to the assigned value plus a statistically independent error term, rather than the classical type, in which the assigned value is equal to the true plus an independent error term.6 Berkson-type measurement error often arises in occupational and environmental settings, in which a group-level exposure value (e.g., based on average exposure for a particular job or neighborhood) is assigned to every individual in that group.7–10
In linear models, additive Berkson error with constant variance does not bias exposure–response coefficients but does increase their variance.6,7 However, in log-linear models when the error variance increases as exposure increases, and disease is rare, Berkson error can result in biased exposure–response coefficients, with either over- or under-estimation.8–12 Increasing error variance with increasing exposure is typical of many occupational settings, and there is evidence of this for occupational silica exposures.13 Furthermore, when measurement error increases with exposure, and when cases have higher exposure than noncases, error will also be a function of case status, that is, will be differential. Thus, cases will typically have more measurement error and higher true exposures than noncases in the same jobs.
Absent any gold standard measurement data on a subset of the cohort, Monte-Carlo sampling can be used to simulate the effect of nondifferential Berkson error.13,14 However, the additional complexities of accounting for the correlation of cases status with error are more easily confronted by building the error models as priors into a Bayesian analysis.15–17
We illustrate both these Bayesian applications in a large, multicentric, IARC-sponsored analysis of 10 silica-exposed cohorts.18
We used data from a previously published exposure–response analysis of lung cancer due and silica exposure in 10 retrospective silica-exposed cohorts with a total of 65,980 workers and 1,069 lung cancer deaths.18 The 10 cohorts included gold miners, granite workers, industrial sand workers, and other types of miners. Cohorts with appreciable exposures to other known lung carcinogens were excluded. Quantitative exposure estimates over time for all jobs held by workers in these cohorts were available19 via job-exposure matrices for each cohort. Application of the matrices to the work histories of each individual worker enabled the estimation of exposure levels across all work histories, and hence the quantitative estimation of cumulative exposure for all workers. For the pooled cohort data, the median cumulative lagged exposure was 2.42 mg/m3-years (IQR: 0.31, 7.78) among cases and 1.99 mg/m3-years (IQR: 0.19, 6.86) among matched controls.
We use nested case–control design via incidence density sampling and the conditional logistic regression likelihood L(β) to estimate the exposure–response coefficient (β) for silica exposure and lung cancer mortality:
where I is the number of cases, Ji is the number of individuals in the risk set for case i, xi1 is the silica exposure for case i, and xij is the silica exposure for individual j from the risk set for case i. The model was fit using both untransformed and log-transformed cumulative silica exposure (adding 1 mg/m3-day before converting to mg/m3-years and taking the log), both with 15-year lags. Up to 10 matched controls from the same cohort were randomly selected for each case, matching on sex, race, and age (within 5 years). We used 10 controls per case but note that the originally published analysis used 100 controls per case.18
We used laboratory animal data to characterize the prior on β, rather than using a naive diffuse prior. Our inclusion criteria were that each study had to have rodents exposed via inhalation and lung cancer as an endpoint, with observation over a 24-month period. Studies were identified by reviewing IARC’s Monographs 68 and 100C on silica.20,21 Three studies summarized in Monograph 68 fulfilled our criteria.22–24 Three other rodent inhalation studies were excluded due to lack of information regarding the exposure concentration, the length of exposure, or the length of follow-up.
Regarding the three included studies, Spiethoff et al.22 observed two groups of 82 female Wistar rats exposed to either 6 or 30 mg/m3 crystalline silica for 29 days, and after 24 months 8 and 13 rats developed lung cancer in the low- and high-exposure groups, respectively. Among 85 nonexposed rats none developed carcinoma of the lung. Muhle et al.23 observed one group of 50 males and 50 female Fischer rats for 25.5 months, during which time they were exposed to 1 mg/m3 crystalline silica for 96 weeks, and among which 13 lung carcinomas developed. In a comparison group of 50 male and 50 female nonexposed rats, there was one lung carcinoma. Holland et al.24 observed 60 female Fischer rats exposed to 12 mg/m3 crystalline silica for the rats’ lifespan (assumed to be approximately 24 months), among which 14 developed lung carcinomas, compared with 0 of 69 nonexposed rats. To establish a prior distribution, we ignored potential effects of different rat strains or genders, and pooled the data from the three studies.
We ran logistic models using the five dose groups from the combined studies (0, 25, 96, 124, and 996 mg/m3-weeks). Time-dependent data were not available to compute lagged exposure or rate ratios, so we used cumulative exposure without any lag, and assumed that the odds ratio from cumulative incidence data was approximately equal to the corresponding rate ratio, a reasonable assumption for rare outcomes.25 The predictor was cumulative exposure to silica (mg/m3-years) or log cumulative exposure, with a lifespan conversion assuming that one rat day was equivalent to 30 human days.26 In the log model, 1 mg/m3-year was added to cumulative dose to avoid taking the log of 0.
To develop a Bayesian prior for the human exposure–response coefficient from the animal exposure–response coefficient, we adopted a probabilistic model for a cross-species extrapolation factor as suggested by a recent National Research Council (NRC) panel which reviewed the EPA’s Integrated Risk Information System.27 Our model assumes that the hazard ratio for humans (HR) is equal to the hazard ratio for rats (HRrat) multiplied by a cross-species extrapolation factor U. It follows that, on the log scale, the human exposure–response coefficient (β) is equal to the rat exposure–response coefficient (βrat) plus the log of U, that is,
We chose a normal distribution for βrat, with the mean and standard deviation set to the point estimate and standard error from logistic regression of the rat data. Following the NRC recommendation, we chose a normally distributed prior on ln(U), centered on 0 and with a standard deviation of σU, that is,
where σU = 1.175, 0.561, or 0 for 10-, three-, or one-fold uncertainty, respectively.27 For example, with σU = 1.175 (10-fold uncertainty), the central 95% of the probability distribution for U ranges from 1/10 to 10, and the median value is 1. Cross-species extrapolation factors of 3 and 10 are typically used for noncancer endpoints27,28 but could also be reasonable for carcinogens given the good concordance between rat bioassays and human data.
As noted above, we focused on two sources of measurement error, both of which were errors of the Berkson type. The first was measurement error due to assignment of job-specific means to individual workers. In the silica cohort, the exposure level for specific subjects was not measured directly, but was assigned based on measurements for a sample of subjects in the same job, or based on estimated levels for specific jobs in the past when no measurements were available, via a job-exposure matrix.19 We assumed that the assigned job-specific mean was equal to the true mean exposure level for workers in that specific job. However, assignment of exposure level to individual workers using the mean level for specific jobs necessarily results in error in the assigned exposure level compared with the true exposure for each individual.
The second type of Berkson error was error due to incorrect conversion of dust to silica. Historically, measurements of dust were available to estimate the mean level of exposure per job. However, the silica concentration (mg/m3) is the metric of interest to regulators. Hence Mannetje et al.19 converted dust measurements to silica, based either on data provided in the literature or by the original authors of the six cohort studies. Often a single conversion factor was used (e.g., silica represented 30% of the respirable dust in the South African cohort), although sometimes the conversion factor varied by time period or facility. If a single conversion factor was used for a cohort, it may have failed to take into account differences in conversion factors by work area or job (e.g., if different areas in a mine had different kinds of dust). The use of a single conversion factor across jobs and work areas is likely to have introduced additional measurement error, albeit at the group (job) level rather than the individual level.
Error Model for Assignment of Mean Job Level to Each Worker
We assumed that the assigned exposure level for each worker for each work history (a job/area/time period combination specific for each worker) was in error. We assumed that the true individual worker exposure level varied log normally around the level μj assigned to all workers in each job/area/time period combination (hereafter referred to as “job”), with a constant coefficient of variation k:
where σj is the standard deviation of personal exposures for job j. The log mean (μlnj) and log standard deviation (σlnj) of a log-normal distribution have the following relationship29 to its arithmetic mean (μj) and standard deviation (σj):
which simplify to the following for σj = kμj:
We used k = 0.80, estimated from 37 measurements at one worksite of an industrial sand plant.30
Error Model for Conversion of Dust Concentrations to Silica Concentrations
A coefficient of variation of about 50% was reported for the dust to silica conversion factor based on 500 samples across 23 jobs in a South Africa gold mine31; we assumed that this degree of variation was applicable to all of our cohorts and normally distributed (with truncation to prevent unrealistic outliers). Other data are available for silica-exposed workers,32,33 which tend to support variation in this range for job-specific conversion factors.
Thus, the true exposure zij experienced by person i during job j is described with a series of probability models incorporating both types of error (group mean assignments and dust to silica conversion):
where xj is the measured group mean exposure for job j (shared by all workers with the same job/area/time period). Cumulative exposure for person i is computed as the product of zij and the time spent in job j, summed over all jobs held by that person up until 15 years before the case age. Notably, workers who share the same job/area/time period have the same conversion error, resulting in a shared-error structure that contrasts with independence assumptions underlying much of the measurement error literature.9,10
As noted earlier, measurement errors are differential by case status when error increases with level of exposure. However, the correlation coefficient between error and case status depends on the unknown exposure–response coefficient. One potential solution to this problem is a Bayesian analysis with a model that directly incorporates any correlations between response and error induced by the exposure–response relationship, without any need for the investigator to specify the correlation. Hence, we conducted a Bayesian analysis using the animal priors and Berkson measurement error structure described in the previous sections, with the two sources of exposure measurement error included as latent variables.
Samples from the posterior distribution of β were generated using Markov chain Monte Carlo techniques implemented with Just Another Gibbs Sampler software and the R programming environment.34,35 Because job was one of the two levels of measurement error and individuals had often worked at more than one job, cumulative exposures were repeatedly reconstructed from individual work histories with different draws from the latent measurement error variables at each Markov chain iteration. Models were run for 10,000 iterations with a 1,000 iteration burn-in period. Three separate chains were used to evaluate convergence using the Gelman-Rubin diagnostic, and means and standard deviations were computed for the posterior samples. Computer code for several of the analyses is provided in the eAppendix 1 (https://links.lww.com/EDE/B141).
The Figure shows the rodent carcinogenesis data, in terms of log odds, after converting a rat lifespan to a human lifespan, for both the observed data and the results of the cumulative exposure and log cumulative exposure models.
For the log cumulative dose model, the exposure–response coefficient βrat from logistic regression with the rat data was 0.51 per mg/m3-year (standard error 0.085), about nine times the corresponding hazard coefficient of 0.061 from frequentist analysis of the human data (Table 1). βrat for the cumulative dose model was 0.003 per mg/m3-year (standard error 0.001), about one-third of the corresponding human coefficient of 0.009. Although the animal data were sparse, which limited assessment of goodness of fit, the log cumulative dose model fit better than the untransformed cumulative dose model (Akaike information criterion 291.0 vs. 313.9).
Table 1 shows the results of the epidemiologic analyses, both for a frequentist analysis (no use of the animal data) and Bayesian analysis including a prior based on βrat and one-, three-, and 10-fold uncertainty in the cross-species extrapolation factor. Results are shown for pooled analysis of all epidemiologic studies and for two selected studies. For the pooled analyses using rat priors with three- and 10-fold uncertainty, the frequentist and Bayesian results were nearly identical. However, with one-fold uncertainty, the rat data were more influential; for each exposure metric the posterior mean is drawn toward the point estimate for βrat.
For individual epidemiologic studies, the impact of the animal prior varied. Results are shown for two example studies in Table 1; these studies were selected to highlight potential impacts of Bayesian incorporation of animal data. For example, for study 4 using log exposure, the frequentist coefficients were 0.160 per log mg/m3-year for humans and 0.51 per log mg/m3-year for rats. Each of the Bayesian posterior means for study 4 (0.174, 0.185, and 0.408 per log mg/m3-year) are weighted averages of the human and rat coefficients, with more weight given to the human coefficient as the uncertainty increases from one- to three-fold to 10-fold. Similar patterns are evident for each study and for both exposure metrics, although for many studies there is little difference in posterior means using three- or 10-fold uncertainty.
Table 2 shows the effect of a Bayesian adjustment for exposure measurement error, along with the rat prior using three-fold uncertainty. In the pooled analysis, there are only modest effects from measurement error adjustment for either of the two exposure metrics. For log exposure, the effect of measurement error adjustment is negligible, and for untransformed exposure there is only a 5% decrease in the posterior mean from 0.0092 to 0.0087 per mg/m3-year, that is, a slight move toward the rat coefficient of 0.003. Similarly, measurement error adjustment has modest impact in study 4, for which the posterior mean decreases by about 3% for each exposure metric. In contrast, in study 9, error adjustment results in an 18% decrease in the posterior mean for the effect of untransformed exposure. In each case, adjustment for measurement error in the human studies moves the posterior mean toward the rat coefficient (although in some cases the movement is negligible).
For some environmental toxicants, both human and animal studies are available to inform researchers and regulatory agencies about potential health effects of exposure. However, it is unclear how to weigh these different types of information in modeling exposure–response relationships. Animal studies conducted using randomized assignment to treatment, blinding, and other good laboratory practices have few threats to internal validity but questionable relevance for estimating exposure–response coefficients in humans. Human studies are conducted in more relevant populations but are typically conducted using observational designs in which exposures and confounders can be difficult to identify and measure accurately. Quantitative methods are needed for combining data from multiple studies in humans and animals in ways that respect their relative strengths and weaknesses, especially by regulatory agencies charged with developing a single exposure–response coefficient or “safe” exposure level for each toxicant.27
Although previous analyses have used Bayesian methods to combine toxicity data from human and animal studies,36,37 here we applied a new method of incorporating cross-species uncertainty using probabilistic interpretations of traditional extrapolation factors commonly used in risk assessment. Although we used software to generate posterior samples, for log-normal uncertainty models without measurement error adjustment (Table 1), the same results can be obtained without specialized software using inverse-variance weighted averages.27 These methods may be useful for a variety of chemicals for which both human and animal data are available to support quantitative dose–response modeling.
We have also accounted for the effects of Berkson exposure measurement error of the type and degree specified by our model, which can contribute bias and uncertainty in the exposure–response coefficient.8 Although we do not have any gold standard measurement to use as a basis for correcting for measurement error, we used Bayesian methods to go beyond a Monte-Carlo simulation of random error, via incorporation of the dependence of error on disease status (i.e., differential error).
Although measurement error adjustment and the animal priors with three- or 10-fold uncertainty did not have a major impact on the posterior exposure–response coefficient in our large sample pooled analyses, they did influence the results for some of the individual cohorts. This is consistent with general principles of Bayesian inference, in which the prior is given less weight as the sample size increases and the data are given less weight as measurement error increases. In situations for which the human studies are fewer or have smaller sample sizes, measurement error adjustment and the animal prior are expected to have larger impacts.
We did not restrict the parameter space for β. Some toxicologic dose–response models exclude any possibility of protective effects, which could be accomplished by truncating the prior at 0 or using a log-normal or gamma prior. For each of our pooled analyses, less than 0.02% of the posterior distribution falls below 0, and truncation would have negligible impact. However, it could have more of any impact on study 4, for which negative values constitute up to 20% of the posterior distribution.
Our model includes several prior parameters based on limited sampling, such as the 80% coefficient of variation for random error in the variation of individual worker exposures around the mean value for each job, based on measurements of 37 workers at one worksite. Although we did not account for it in our analyses, each of these parameters has some sampling uncertainty that could be addressed using an additional probability distribution and hyperparameters.17 Given the modest impact of our measurement error model on the results, we suspect that incorporating the additional sampling uncertainty would have minimal impact on our findings.
Another limitation of our analysis is the use of a 15-year lag to compute cumulative exposure in humans, but no such lag for rats. Incorporating a lag could decrease the cumulative exposure in one or more of the rat studies, which would increase the prior hazard coefficient (βrat). This is unlikely to affect the posterior for the pooled data, which appears largely unaffected by the prior when the cross-species extrapolation factor is greater than 1, but could have some impact on the results for the smaller studies.
We did not account for some potential threats to the validity of the exposure–response estimates, such as measurement error in the confounders, healthy worker effects, non-Berkson (i.e., classical) exposure measurement error, or unidentified coexposures. We chose to focus our efforts on exposure measurement error and cross-species extrapolation as the two most important sources of uncertainty for the studies on silica and lung cancer, but in other epidemiologic settings (particularly with nonoccupational exposures) uncontrolled confounding may be a more serious threat to validity. In theory, a probability model could be specified for any of these threats in each study, as suggested by Lash and Fink,13 and incorporated into the Bayesian analysis to combine the information across the human and animal studies.
The authors thank Dr. Kate Guyton of the International Agency for Research on Cancer for her valuable help in identifying and interpreting the key animal silica inhalation studies.
1. International Agency for Research on Cancer. Report of the IARC Advisory Group To Recommend On Quantitative Risk Characterization. 2014. Lyon, France: International Agency for Research on Cancer, (Internal Report 14/001). Available at: http://monographs.iarc.fr/ENG/Publications/internrep/14-001.pdf
. Accessed 14, November 2015.
2. Huff J. Long-term chemical carcinogenesis bioassays predict human cancer hazards. Issues, controversies, and uncertainties. Ann N Y Acad Sci. 1999;895:56–79.
3. Huff J, Jacobson MF, Davis DL. The limits of two-year bioassay exposure regimens for identifying chemical carcinogens. Environ Health Perspect. 2008;116:1439–1442.
4. US Environmental Protection Agency. Guidelines for Carcinogen Risk Assessment. 2005. Washington, DC: US Environmental Protection Agency; http://www2.epa.gov/sites/production/files/2013-09/documents/cancer_guidelines_final_3-25-05.pdf
. Accessed 14 November 2015.
5. US Food and Drug Administration, Center for Food Safety and Applied Nutrition. Toxicological Principles for the Safety Assessment of Food Ingredients. 2007. Washington, DC: US Food and Drug Administration; Available at: http://www.fda.gov/downloads/Food/GuidanceRegulation/UCM222779.pdf
. Accessed 14, November 2015.
6. Snedecor G, Cochran W. Statistical Methods. 1989.Ames, IA: Iowa St. University Press.
7. Seixas NS, Sheppard L. Maximizing accuracy and precision using individual and grouped exposure assessments. Scand J Work Environ Health. 1996;22:94–101.
8. Armstrong BG. Effect of measurement error on epidemiological studies of environmental and occupational exposures. Occup Environ Med. 1998;55:651–656.
9. Deddens J, Hornung R. Smith C, Christiani D, Kelsey K. Quantitative examples of continuous exposure measurement errors that bias risk estimates away from the null. Chemical Risk Assessment of Occupational Health. 1995:London: Auburn; 77–85.
10. Steenland K, Deddens JA, Zhao S. Biases in estimating the effect of cumulative exposure in log-linear models when estimated exposure levels are assigned. Scand J Work Environ Health. 2000;26:37–43.
11. Prentice R. Covariate measurement errors and parameter estimation in a failure time regression model. Biometrika. 1982;69:331–341.
12. Armstrong BG. The effects of measurement errors on relative risk regressions. Am J Epidemiol. 1990;132:1176–1184.
13. Lash TL, Fink AK. Semi-automated sensitivity analysis to assess systematic errors in observational data. Epidemiology. 2003;14:451–458.
14. Phillips CV. Quantifying and reporting uncertainty from systematic errors. Epidemiology. 2003;14:459–466.
15. Gustafson P. Measurement Error and Misclassification in Statistics and Epidemiology: Impacts and Bayesian Adjustments. 2004.Boca Raton, FL: Chapman & Hall/CRC.
16. Spiegelhalter D, Thomas A, Best N, Lunn D. WinBUGS Version 1.4.3. 2007. Available at: http://www.mrc-bsu.cam.ac.uk/bugs
. Accessed December 22, 2016.
17. Espino-Hernandez G, Gustafson P, Burstyn I. Bayesian adjustment for measurement error in continuous exposures in an individually matched case-control study. BMC Med Res Methodol. 2011;11:67.
18. Steenland K, Mannetje A, Boffetta P, et al.; International Agency for Research on Cancer. Pooled exposure-response analyses and risk assessment for lung cancer in 10 cohorts of silica-exposed workers: an IARC multicentre study. Cancer Causes Control. 2001;12:773–784.
19. Mannetje A’, Steenland K, Checkoway H, et al. Development of quantitative exposure data for a pooled exposure-response analysis of 10 silica cohorts. Am J Ind Med. 2002;42:73–86.
20. IARC (International Agency for Research on Cancer), Arsenic, Metals, Fibres and Dusts. 2012.Geneva, Monograph 100C, WHO.
21. IARC (International Agency for Research on Cancer), Silica, Some Silicates, Coal Dust and para-Aramid Fibrils. 1997.Geneva, Monograph 68, WHO.
22. Spiethoff A, Wesch H, Wegener K, Klimisch HJ. The effects of Thorotrast and quartz on the induction of lung tumors in rats. Health Phys. 1992;63:101–110.
23. Muhle H, Kittel B, Ernst H, Mohr U, Mermelstein R. Neoplastic lung lesions in rat after chronic exposure to crystalline silica. Scand J Work Environ Health. 1995;21(suppl 2):27–29.
24. Holland L, Wilson J, Tillery M, Smith D, Goldsmith D, Winn D, Shy C. Lung cancer in rats exposed to fibrogenic dusts. Silica, Silicosis, and Cancer: Controversy in Occupational Medicine. 1986.New York, NY: Praeger.
25. Rothman K, Greenland S, Lash T. Chapter 3: Measures of occurrence, modern epidemiology. 2008.Philadelphia, PA: Lippincott, Williams, and Wilkins.
26. Sengupta P. The laboratory rat: relating its age with human’s. Int J Prev Med. 2013;4:624–630, Table 2.
27. National Research Council, Chapter 7 in Review of EPAs Integrated Risk Information System process. 2014.Washington, DC: NRC.
28. Dourson ML, Stara JF. Regulatory history and experimental support of uncertainty (safety) factors. Regul Toxicol Pharmacol. 1983;3:224–238.
29. Atchinson J, Brown J. The Lognormal Distribution. 1963.Cambridge, England: Cambridge University Press.
30. Amandus H. A feasibility study of the adequacy of company records for a proposed NIOSH study of silicosis in industrial sand workers, NIOSH report, NIOSH DRDS. January 30, 1990.Morgantown, W Va.
31. Churchyard G, Pemba L, Magadia B, et al. Silicosis prevalence and exposure response relationships in older black miners in a S. African gold mine, SIMRAC report, March 7, 2003.
32. Kreiss K, Zhen B. Risk of silicosis in a Colorado mining community, Am J Ind Med. 1996;30:529.
33. Sanderson W, Steenland K, Deddens J, Historical respirable quartz exposures of industrial sand workers: 1946–1996. Am J Ind Med. 2000;38: 389–398.
34. Plummer M. JAGS Version 3.4.0 User Manual. 2013. Available at: http://www.stats.ox.ac.uk/~nicholls/MScMCMC14/jags_user_manual.pdf
. Accessed 14 November 2015.
35. Plummer M. rjags: Bayesian Graphical Models using MCMC. R package version 3–15. 2015. Available at: http://CRAN.R-project.org/package=rjags
. Accessed 14 August 2015.
36. Peters JL, Rushton L, Sutton AJ, Jones DR, Abrams KR, Mugglestone MA. Bayesian methods for the cross-design synthesis of epidemiological and toxicological evidence. J R Stat Soc C-Appl Stat. 2005;54:159–172.
37. DuMouchel W. Multivariate Bayesian logistic regression for analysis of clinical study safety issues. Stat Sci. 2012;27:319–339.