Systematic Dose-Response of Environmental Epidemiologic Studies: Dose and Response Pre-Analysis

Meta-analysis approaches can be used to assess the human risks due to exposure to environmental chemicals when there are numerous high-quality epidemiologic studies of priority outcomes in a database. However, methodological issues related to how different studies report effect measures and incorporate exposure into their analyses arise that complicate the pooled analysis of multiple studies. As such, there are “pre-analysis” steps that are often necessary to prepare summary data reported in epidemiologic studies for dose-response analysis. This paper uses epidemiologic studies of arsenic-induced health effects as a case example and addresses the issues surrounding the estimation of mean doses from censored dose- or exposure-intervals reported in the literature (e.g., estimation of mean doses from high exposures that are only reported as an open-ended interval), calculation of a common dose metric for use in a dose-response meta-analysis (one that takes into consideration inter-individual variability), and calculation of response “effective counts” that inherently account for confounders. The methods herein may be generalizable to 1) the analysis of other environmental contaminants with a suitable database of epidemiologic studies, and 2) any meta-analytic approach used to pool information across studies. A second companion paper detailing the use of “pre-analyzed” data in a hierarchical Bayesian dose-response model and techniques for extrapolating risks to target populations follows.


Introduction
While dose-response analyses for the purpose of human health risk assessment have typically been conducted by applying dose-response modeling techniques to a single, perhaps "best" animal toxicological data set, data from epidemiologic studies have been used less frequently. As of 2007, approximately 8% (44) of the 545 chemical assessments in the IRIS database were identified as using human data in the derivation of toxicity values (Persad and Cooper 2008) and several more recent assessments (ethylene oxide, hydrogen cyanide, Libby amphibole asbestos, 2,3,7,8-tetrachlodibenzo-p-dioxin, tetrachloroethylene and trichloroethylene) make use of epidemiologic data for quantitative purposes (U.S. EPA 2019). The use of epidemiologic data is preferable over the use of animal data because it avoids issues related to interspecies extrapolation to account for differences in toxicokinetics and toxicodynamics (U.S. EPA 2005). However, when using epidemiologic data, there is often a need to integrate information from a variety of sources. Meta-analysis explicitly involves simultaneous modeling of studies, yielding pooled estimates, in order to leverage information from multiple studies, possibly further decreasing uncertainties surrounding the use of individual studies (i.e., human variability, etc.).
However, it is often the case that risk assessors are limited to information they can obtain from the peer-reviewed literature and overwhelmingly this results in only summary data being available. As such, the following methods should prove valuable for increasing the role of epidemiologic studies in quantitative dose-response analyses in human health risk assessment. This is the first in a series of two papers that will address the use of epidemiologic studies for dose-response analyses for the purpose of identifying target population risks in a human health risk assessment. The approaches outlined in these two papers expand upon previous methods applied by EPA in IRIS assessments for the evaluation of epidemiologic data (U.S. EPA 2017) in several respects. An important feature of this approach is that the studies suitable for combination in our dose-response meta-analysis (DRMA) can be either cohort studies or case-control studies. In particular, we want to emphasize and illustrate what we consider to be important, and potentially novel, components of such an analysis through a case study of inorganic arsenic drinking water exposure and bladder cancer morbidity. Figure 1 displays the overall flow of the DRMA analysis. In this first manuscript, we address issues that typically arise when dealing with summarized (e.g., published) analyses, shown in blue in Figure 1. We categorize these as "pre-analysis" steps, in the sense that they prepare the data so that they are amenable to dose-response modeling and characterize uncertainties associated with the data presentations available to the modeler. These steps are preliminary steps that allow for a systematic dose-response analysis that combines data from multiple studies.
The subsequent paper will focus on the choice of models, the statistical rationale for their use, and how they can be suitable for studies of different types (cohort and case-control). Details of a Bayesian DRMA will be provided, including important derivations that allow or facilitate analysis (orange section of Figure 1). That paper will also demonstrate how the results of the dose-response modeling will be extrapolated to a target population of interest using lifetable analysis approaches (shown in green in Figure 1).
All of the steps illustrated on Figure 1 will be performed in the context of epidemiologic studies relating arsenic exposure (through drinking water) to bladder cancer incidence. The specific studies included in this analysis are discussed briefly in the following section. Subsequent sections of this manuscript will address dose estimation for interval-censored data, including the important consideration of uncertainty in dose metric values, and the use of reported numbers of cases to calculate effective counts that reflect adjustments for potential confounders. This manuscript ends with a discussion of pertinent issues and of the linkage to the DRMA implementation detailed in the second paper in this series.

Case-Study: Epidemiologic Studies of Arsenic and Bladder Cancer
The case study underlying the methodologies discussed here is based on published reports that link arsenic exposure to bladder cancer morbidity. Multiple health agencies have concluded that the relationship between arsenic exposure and bladder cancer is causal (NRC, 2013;ATSDR, 2016;ATSDR 2007;NTP, 2016;IARC 2012;WHO 2011a,b); we take that as a given for our analyses.
Multiple bladder cancer studies considered to have sufficient information to support a doseresponse analysis will be used in the actual dose-response portion of this analysis (see second paper, Allen et al. 2020). It should be noted that these studies do not necessarily represent the final set of epidemiologic studies that might be chosen for a formal human health risk assessment. For the purposes of this paper (i.e., demonstration of the pre-analysis steps), we have focused on only two such studies, a cohort study (Chen et al. 2010) and a case-control study (Meliker et al. 2010), as the pre-analysis steps differ slightly given different study designs. The second paper (Allen et al., 2020) describing the Bayesian DRMA methods will provide more details on other studies included in the analysis.
In this case, the two studies selected here to demonstrate the dose and response conversions each had three or more exposure groups defined in the published reports, including the "referent" group (i.e., the group to which other groups are compared in order to calculate the relative risk or odds ratio). It is worth noting, however, that in a DRMA context, there is no need to exclude studies with only two exposure groups, nor would inclusion of two-group studies necessarily limit the "complexity" of the models applied. The fact that several studies are modeled together provides data at more than two exposure levels and allows for more inclusiveness with respect to study selection and dose-response model forms.

Pre-analyses for the Estimation of a Common Dose
One of the novel aspects that distinguishes the current approach from other meta-analyses of epidemiologic data, such as the recently published arsenic meta-analyses for bladder and lung cancer by Lynch et al. (2017a;b) and for cardiovascular diseases by Moon et al. (2017), is the conversion of route-specific (e.g., oral) exposure metrics commonly used in epidemiologic risk assessments (e.g., drinking water concentrations, measures of cumulative exposure) to a common dose metric (e.g., total daily oral intake) that is presumed to be more closely related to the responses for the health effect in question. Conversion of disparate exposure metrics to a single unified dose metric removes the restriction of choosing a "best" exposure metric and thus allows for a larger set of studies to be used in the DRMA. When drinking water is the primary exposure route, Baris et al. (2016) emphasized the importance of incorporating estimates of water intake for the reliable discernment of dose-response relationships, particularly for low to moderately exposed populations. As examples, the cohort study (Chen et al. 2010) and the case-control study (Meliker et al. 2010) vary with respect to the units in which exposure or dose were presented: ug/L-year and ug/day, respectively. As such, a primary requirement for the pre-analysis is the specification of a common unit of dose and conversion of all study doses to that scale. The units of interest for this analysis were average daily μg inorganic arsenic (iAs) intake per kg bodyweight (average daily μg/kg).
In practice, we combined the requirement for a common dose metric with the need to estimate a mean dose for each group that is merely represented by a range of dose values in the actual studies. In most published reports, groups are commonly defined by the range of exposures experienced by the individuals (or person-years) in the group; they are "interval censored" with respect to exposure. Thus, we first computed means (in units as given in the published papers) and then applied the conversions to those means. This process is summarized in Figure 2, which shows that accounting for uncertainties in the estimation of the means or of the conversions is an integral part of this exercise. For the mean estimation, the process is summarized as follows.

Group Means and Uncertainty
Suppose that a study has G exposure groups defined as (c 0 , c 1 ), (c 1 , c 2 ), … , (c G−1 , c G ); c 0 may equal 0 and c G may equal ∞. Let τ, with elements τ g , be the vector of observed proportions of individuals in groups g = 1, …, G. As discussed below, for a case-control study, only the control counts were used to define the τ values; otherwise all individuals were used to define τ.
We assume that individual exposures follow a lognormal distribution with parameters μ and σ (respectively, the mean and standard deviation of natural logarithms of exposure); the same assumption was made by Moon et al. (2017) in their analysis of arsenic risks. We estimate these parameters by maximum likelihood (ML), maximizing the following loglikelihood expression appropriate for the censored data: where N is the total number of observations under consideration and φ(x | μ, σ) is the cumulative distribution function for the Normal distribution with mean μ and standard deviation σ. Let μ' and σ' be the ML estimates (MLEs) given τ and the values of c g (g = 0, …, G). A simple likelihood maximization routine, implemented with an Excel spreadsheet (see Supplemental Material), was used to estimate μ' and σ'. Given μ' and σ', the mean within a given exposure interval (c g , c g+1 ) is given by mean(g) = e μ′ + σ′ 2 2 × θ(U 1 (g) − σ′) − θ(U 0 (g) − σ′) θ(U 1 (g)) − θ(U 0 (g)) , Eq. 1 where U 1 (g) = (ln(c g + 1 ) − μ′) σ′ , U 0 (g) = (ln(c g ) − μ′) σ′ , and θ() is the cumulative distribution function for the standard normal distribution (Söderlind 2013).
In addition to the maximum likelihood estimate of an exposure for each dose group, we have computed what are referred to as "low-end" and "high-end" exposure estimates. Those estimates were based on profile likelihood bounds for the mean exposure for individuals in the high-exposure group. First, we computed the maximized log-likelihood (MLL) (the loglikelihood associated with MLEs of the parameters μ and σ). Then the mean value for the high-group was maximized by modifying μ and σ subject to the constraint that where LL is the log-likelihood associated with the modified μ and σ. This is a chi-squaredbased (1 degree of freedom) 95% confidence upper bound on the high-group mean. The same procedure was followed to determine a 95% confidence lower bound on the highgroup mean by minimizing the mean value for the high group under the same constraint (Eq. 2). Supplemental Files S3 and S4 demonstrate fully these calculations.
Note that we are only specifically maximizing (or minimizing) the mean for the highest exposure group. In some instances, this may result in a change in the estimated mean for other groups (primarily the lowest exposure group) when the modified μ and σ are input into Eq. 1. We have opted for this procedure for two primary reasons. First, the highest exposure group is almost always the one that is most uncertain, often presented as an open-ended interval. The mean for that group is not easily nor consistently estimated by techniques that can be applied to the other exposure groups (e.g., by taking the midpoint of the interval). Second, the dose and response observations for the highest exposure group often have a large impact on dose-response model estimation. Especially for the slope parameters in a model, the positioning of the highest dose-response pair can be very influential. Thus, our procedure focuses on that influential observation and derives bounds specifically tailored to determine the full range of possible values for exposure consistent with the data (and our assumption of lognormality of exposures over the studied population) for that observation. Additional discussion about the choice of the lognormal distribution is provided in the Discussion section.

Dose Conversions and Uncertainty
As discussed above in Section 2.2, estimation of a common dose metric (as opposed to exposure metric) for all studies is imperative in order to account for the effect of dose on the estimated response. However, epidemiologic studies often use different, but related, exposure-or dose-metrics such as exposure concentration, cumulative exposure, and even biomarkers of exposure such as internal tissue concentrations or urinary concentrations of the chemical of interest. In this specific case-report, if the data presented in all the published studies had been in the units of interest (daily average μg/kg), no conversion would be necessary and DRMA could be performed with the study-reported data. However, that was not the case for any of the studies under consideration here, nor would it be in general practice.
First, all of the factors required to convert to μg/kg were identified. As an example, consider the Chen et al. (2010) study of bladder cancer. The tabulation of results reported in that study reported cumulative drinking water exposures, in (μg iAs/L drinking water) · years. The conversion from cumulative exposures to average daily dose (μg/kg) was carried out as follows: where the terms in that expression are: DI = dietary intake (average daily μg/kg); f = fraction of time (over lifetime up through the study) spent consuming well water (unitless); WCR = water consumption rate (L/kg); WE = well water concentration (μg/L); and LE = low-end water concentration (μg/L). The variable f was calculated as the ratio of the assumed average duration of well exposure (ADWE), generally the reported duration of drinking well water (RDWE; yr), to the average age at diagnosis (AAD; yr). It was assumed that when drinking non-well water, the subjects consumed water with the low-end water concentration. The parameter WE was derived separately for each group by dividing the reported cumulative exposure (CE; μg/L · yr) for that group by the RDWE: WE = CE RDWE . The values used in Eq.
3 above ideally would come from study-specific data reported in the study of interest, but could also be drawn from other suitable sources (see Table 1  Importantly, the factors affecting conversion were not treated as single values. Rather, a distribution was assumed over the individuals in the study and the specific type of distribution was allowed to differ for different variables. For instance, lognormal distributions were assumed for DI, WCR, RDWE, and LE, whereas a beta distribution was used for the f variable. Those distributions were expressed in terms of inter-individual variability and were based on knowledge of the population from which the study participants were drawn (e.g., exposure factors handbook values or study-specific data). Thus, when converting the group-specific mean cumulative values to group-specific mean μg/kg values, we sampled from the relevant distributions (in the example above, for DI, f, WCR, LE, and RDWE) and computed the resulting conversion for each individual in the exposure group in question (in cases where the n per dose group was greater than 1000, the average was computed using 1000 samples to ease computational burden). Then we averaged over all the individuals in the group, to come up with a group mean conversion (see Table 2). A Monte Carlo (MC) simulation was conducted in which this process was repeated 1,000 times to derive a distribution of averaged dose values. The median as well as the 2.5 th and 97.5 th percentiles were pulled from the MC results to characterize the best, low-end and high-end dose values (Figure 2; Tables 3 and 4). As expected, given the methods used to estimate dose group means (see Section 2.3), the resulting differences in the high exposure group mean were substantial when using either the "MLE", "high", or "low" dose scenarios whereas the low exposure group mean distribution were not (see Table 4, Figure 3).
Other items to note include the fact that we have included dietary exposure in all the conversions. Thus, our dose estimates represent total inorganic arsenic intake from all oral routes, not necessarily just exposure from drinking water. Our approach explicitly considers the sample sizes in each exposure group; by averaging over the number of individuals in each group, we have automatically considered that groups with more observations will be less uncertain (about the mean group-specific conversions) than groups with fewer observations.
Finally, note that each study, possibly reporting different exposure summaries, is handled differently depending on the reported units. For example, the Meliker et al. (2010) study used daily arsenic exposure (in units of μg/day) as the dose metric rather than cumulative exposure in units of (μg iAs/L drinking water) · years as did Chen et al. (2010). Thus, the conversion to average daily μg/kg (dose) was carried out for the Meliker et al. (2010) study as follows: where the terms in that expression are as above with the addition of BW = body weight (kg).
The variable f was estimated as described above, but in the case of the Meliker et al. (2010) study, the parameter WE was derived separately for each group by dividing the reported daily exposure (μg/day) for that group by a BW value selected from a normal distribution centered at the average bodyweight from the EPA Exposure Factors handbook (2011).
All dose conversions were computed via Excel, using the MC simulation add-in Yasai (www.yasai.rutgers.edu). While other software programs or languages (i.e., R or Python) are more powerful in some regards and could be used to implement the dose-conversions, Excel was chosen for its ubiquity of use and because Excel workbooks are useful for organizing the analyses for presentation. The conversion implementations are shown in full detail for all studies used in the Bayesian DRMA analysis in the Supplemental Material (Files S1-S9). Supplemental files S3 and S4 provide full details for the Chen et al. (2010) and Meliker et al. (2010) studies.

Response Pre-Analysis
For both cohort and case-control studies, published manuscripts almost always report adjusted relative risks (RR's) or odds ratios (OR's), respectively. The adjusted results attempt to factor out the effects of other, possibly confounding, variables, in order to estimate the effect specifically associated with the exposure of interest, in this case arsenic exposure.
The Bayesian approach that we have adopted for the dose-response analysis (see second paper) is based on likelihoods of observing a particular number of cases. For example, the number of observed cases in a cohort study is commonly modeled as coming from a Poisson distribution.
To deal with this requirement, we have computed adjusted counts of cases and controls (or cases and expected numbers) -we refer to them throughout this manuscript as "effective count(s)" to avoid confusion with other adjustments that may be part of the analysis. The derivation of effective counts had one very specific goal: to construct a set of counts that reflects only the effect of arsenic. It attempts to construct a data set that would be as if all groups under consideration differed only with respect to arsenic dose but were uniform with respect to the other variables for which RR's or OR's were adjusted. As will be shown below, we are in fact adjusting so as to mimic data that might have been collected had all covariates (other than dose) in all groups been the same as those observed in the referent group. In the context of a case-control study, for example, the effective counts could be viewed as a single 2xG table (G = number of dose groups) of the effective counts of cases and controls which could be considered to represent what data one would have gotten from sets of individuals who were homogeneous with respect to other covariates. Therefore, "effective counts" are the data that would have resulted in the adjusted OR or RR values had confounding not occurred in the study population.
Calculation of effective counts has been the focus of multiple papers: Greenland and Longnecker (1992) reported on a method that retained original sample sizes but obtained the adjusted relative risks, Hamling et al. (2008) used a method that allowed sample sizes to change but resulted in the adjusted RRs and the standard error of the logRR, and Orsini et al.
(2012) provided corrected equations for the variance of the logRR and concluded that either adjustment improved estimation (compared to using unadjusted counts). Both Greenland and Longnecker (1992) and Orsini et al. (2012) make clear that "relative risk" includes the metrics of OR, hazard ratio, etc., and so the methods can be applied to case-control studies, incidence rate studies, or cumulative incidence studies (see Rothman and Greenland (1998), for definitions of these study types, which are discussed further below). The definition of the variance terms for logRR varies across these study types, but they are all amenable to effective count calculations.
The term "effective count" refers to the fact that when making the adjustments to OR and RR, the impact of the other variables is removed. But the estimation of the associations between confounders and dose, and between confounders and the endpoint, "uses up" some of the degrees of freedom associated with the initial sample size. In essence, the approximation of the "otherwise homogeneous set of individuals" mentioned above results in a smaller sample size for evaluating just the arsenic effect. The magnitude of that effect depends (among other things) on how strongly the set of confounders for which RR's or OR's have been adjusted are associated with the arsenic dose and with the endpoint of interest (Rothman and Greenland 1998).
In practice, one typically has two sources of information from published literature from which effective counts can be estimated. The first source consists of the values of the adjusted RR's or OR's themselves. Whatever counts we derive should result in the values of the adjusted ratios reported, when one computes a "simple" ratio from the effective counts. The second source of information is obtained from estimates of the standard errors (or confidence limits) reported for the RR's or OR's (or for log(RR) or log(OR)). Those input values are shown in Tables 5 and 6 for our case study. Examples will illustrate the procedure for effective count computation.

Incidence Rate Cohort
Studies-This approach was implemented for cohort studies where observed numbers of cases and expected numbers of cases were presented (or derivable) and the study used an internal referent group for defining the relative risks.
Consider the data shown in Table 5, for the Chen et al. (2010) dataset. The first group is the internal referent group (adjusted RR = 1 by definition). Orsini et al. (2012) give the standard deviation of the logRR(i) values from such a study as: Here Cases(i) refers to the number of cases in group i, with the referent group being group 0. We will fix the number of cases in the referent group; reasons for that decision are discussed at the end of this section. Then it is easy to solve for Cases(i) in Eq. 5 and then compute the expected 1 number as Cases(i) adjusted RR 2 . Derivation of a value for SE(logRR(i)) is based on the reported confidence limits for the logRR(i) values, as follows. The standard procedure for estimating a 95% upper confidence bound for RR's is this (Rothman and Greenland, 1998 from which Eq.5 gives the effective count for Cases(1), which equals 2.84. In practice (and for results shown in Table 5, for example), we have averaged the results of using the upper bound and lower bound for the confidence interval, to reduce the effect of round-off error in reported values. This is equivalent to equating the width of the confidence interval to 2 × 1.96 × SE (on the log scale). We computed the SE values for the two sides separately (and then averaged) as that facilitates the identification of errors or typographical mistakes in the reported values. The two SE estimates should be essentially the same, to the number of 1 "Expected" here refers to what the expected number in group i would be if it had the same exposure as the reference group, but with its own specific confounder profile. 2 Adjusted RR here refers to the adjusted relative risk value reported in the included studies. since, even with effective counts, RR = number of cases divided by expected number. The same procedure is followed for all the other dose groups. The results obtained for this example are displayed in Table 5. Computation of RR using the effective cases and expected effective counts reported in Table 5 results in RR's and confidence bounds that match those reported as "adjusted" values in Chen et al. (2010), as was expected.
In this particular case, the effective counts for the cases are very similar to the raw counts. Adjustment for the other covariates had little effect on effective counts in this instance. Table 6, obtained from a published report of case-control study of arsenic and bladder cancer (Meliker et al. 2010). Here, we want to replace the case numbers with effective counts (call them a 0 to a 2 , for the three groups with 0 corresponding to the referent group as usual). Moreover, we want to replace the control numbers with effective counts (call them b 0 to b 2 , for the three groups).

Case-Control Studies-Consider the data in
By the definition of odds ratios, and of standard errors for log(OR), we know that we want, in this example, Furthermore, the basis for the confidence limits (Rothman and Greenland,1998) is in the following equations: SE(log(OR 1 )) = 1 a 0 + 1 b 0 + 1 a 1 + 1 b 1 = 0.14 SE(log(OR 2 )) = 1 a 0 SEs are calculated from the confidence limits bounds, with the calculation based on the upper bound averaged with the corresponding calculation based on the lower bound, as was described for the cohort studies (Eq. 5 but with "OR" replacing "RR"). There are four equations for six unknowns (a0,…,a2, b0,…,b2), so we must specify two additional constraints. The first constraint is that we will set a 0 equal to the observed number of cases in the referent group (189 in this example). The second constraint relates to the number of

EPA Author Manuscript
EPA Author Manuscript controls in the referent group, as a proportion of the total number of controls in the 3 groups. Specifically, we specify b 0 ,…,b 2 to satisfy the equation b 0 (b 0 + b 1 + b 2 ) = r 0 .
Eq. 8 The resulting system of equations requires a relatively simple numerical solution. The spreadsheet provided in the Supplemental Material sets up the equations and uses Excel's Solver capability to illustrate how the effective counts were derived. The rationale for this constraint is presented in the Discussion section below.
In this example a 0 = 189 r 0 = 252 (252 + 234 + 48) = 0.472 Table 6 also shows the effective counts that result from the procedure defined above. Using the effective counts reported in Table 6 for computation results in the ORs and confidence limits reported as adjusted values by Meliker et al. (2010), as expected.

Cumulative Incidence Cohort
Studies-For some analyses, the type of cohort study dealt with is called a cumulative incidence study (Rothman and Greenland 1998). In that case, the RR's are ratios of the proportions of subjects getting bladder cancer; the data might be presented as in Table 7 (no cumulative incidence cohort study was used for bladder cancer analyses, the examples herein use data from Chen et al. (2010) presented as if they were derived from such a study). The unadjusted RR's are simply equal to: Cases (0) Eq. 9 where i indicates the group number and, again, group 0 corresponds to the referent group. The equation for SE(log(RR)) is (Rothman and Greenland, 1998) SE(log(RR(i))) = 1

Eq. 10
Together, Eqs. 9 and 10 and define 2 × n equations for 2 × (n+1) unknowns ("unknowns" in the sense that the effective counts that we desire are unknown to us and will result in values for Cases(i) and N(i) for i = 0, …, n, where n is the number of dose groups. Two additional equations are required to specify the system. Extending our approach to the incidence rate studies discussed above, we will fix Cases(0) to be the same as reported in the data (e.g., 6). In addition, we will also fix the value of N(0) (e.g., the sample size in the referent group, 2534). The rationale for that choice is discussed in the next section.
The algebra that yields estimates for Cases(i) and N(i) is relatively straight-forward. These calculations are also illustrated in a spreadsheet provided in the Supplemental Material. As desired the RR's and the confidence limits computed using the effective counts reported in Table 7 match the corresponding "adjusted" values reported as the hypothetical "adjusted" RR values in Table 7, as desired. Coincidently, the effective cases in Tables 2 and 7 are very similar. That is the case because the effective count calculations in those two instances were based on the same data (just treated differently, according to the assumption about whether they are incidence rate or cumulative incidence data).

Discussion
The pre-analysis calculations shown here were part of a larger analysis that used Bayesian DRMA techniques to derive risk estimates (Figure 1). However, it is worth noting that we believe the pre-analyses demonstrated here provide a more reliable basis for analysis using summarized data derived from epidemiologic studies , whether they are analyses of individual studies or DRMA analyses of multiple studies, and whether the actual modeling is frequentist or Bayesian in nature.
When, as is most often the case in risk assessment contexts, one must deal with summarized, published data, the pre-analysis steps described above may be useful. While the lack of raw data is often the case in human health risk assessment for a variety of reasons (individual data not available, study authors unwilling to provide individual participant data), two obvious alternatives would be 1) for study authors to provide the slopes (and confidence intervals) for the association between the exposure and outcome of interest across the entire dose range; or 2) to obtain the original data and perform the dose-response analysis using it,. Even better would be that raw data are automatically provided (with some necessary precautions such as anonymizing personal identifiable information) as supplemental information to published epidemiologic studies. While that approach would be the most appropriate analysis with the least uncertainty surrounding it, methods development efforts must necessarily be cognizant of real-world limitations. One aspect of the reported methods that would retain its usefulness even in the light of fully available individual data is conversion of multiple distinct exposure metrics into a common metric for use in DRMA. However, even there some uncertainty in the analysis could be reduced by using the individual reported exposure, dietary, and study population characteristics data in the dose conversions.
Note that the methods we have used to derive effective counts differ from those proposed by either Greenland and Longnecker (1992) or Hamling et al. (2008). Our methods result in confidence limits that match those reported for adjusted values; the method of Greenland and Longnecker (1992) does not. Moreover, we use different constraints than used by Hamling et al. (2008). The reasons for the different constraints relate to the fact that we are adjusting counts so as to derive effective data that might have been collected had the covariate levels remained the same as in the referent group, for all groups. Thus, we have attempted to preserve the original referent group characteristics to the extent possible. As described below, there is some rationale for this particular method for calculating effective counts that does relate the fact that subsequent analyses will be Bayesian. If that is not the case, then the effective count derivations described by Hamling et al. (2008) should be sufficient. On the other hand, we see no reason that the approach described here should not be satisfactory, regardless of subsequent analysis plans, since it also yields effective counts that will match the reported adjusted RRs and their derived confidence limits.
The Bayesian analysis perspective that will be presented in the second, follow-on manuscript thus influenced our choices. For the incidence rate cohort studies, our adjustments preserve the original count of cases in the referent group, while still making adjustments to the counts in the other dose groups. We have made that choice based on the methodological decision to treat the referent group separate from the other cohort study groups with respect to definitions of prior distributions. As will be described in detail in the second part of this series, we derived priors for the Poisson expected value in the referent group using that group's observations. That group was then not included in the likelihood contributions during subsequent parameter updating. It is the referent group expected value (an uncertain quantity with a prior distribution to be updated), not the observed number in that group, that defines the RR for the non-referent groups. Therefore, the referent group data points are not "counted twice". The point here is that, since we derived priors for the referent group expected number based on the observed count in the referent group, we did choose to derive effective counts that preserve the original number of cases in that group.
For the case-control studies, the second manuscript will explain that the calculation of the log-likelihood contributions for the Bayesian analysis requires an approximation to the distribution of dose in the population. That approximation depends on the proportions of controls across the various dose groups. Those proportions are not something that is affected by covariate adjustment. Therefore, we desired that the method used for effective counts in case-control studies retained as much of the discretized control group distribution as possible. It turns out that it was only possible (with a few rare exceptions) to fix the proportion of controls in the referent group. Actually, the proportion of controls in any group could have been fixed, but the focus on the referent group makes sense when coupled with the fact that we also fixed the referent group number of cases. Thus, while retaining the fixed value of the referent group number of cases, we added the constraint that the proportion of controls in the referent group (relative to the total number of controls over all dose groups) be constant at the observed value. Informal observation for the studies considered in the bladder cancer DRMA suggested that fixing the number of cases and proportion of controls in the referent group did a good job of matching the other control proportions. That latter constraint is one of those that was applied by Hamling et al. (2008).
For the cumulative incidence studies, the characteristics of the referent group now include the group sample size as well as the number of cases. It was a natural extension of the option used with incidence rate studies to fix both the referent group cases and sample size. This, again, would lend itself to a Bayesian analysis where the referent group was treated separately with respect to setting priors and updating parameter distributions.
It is also worth noting that the approaches described here for dose estimation automatically factor in uncertainties. That includes uncertainties associated with summarized reporting (ranges for exposure groups are given but means for those groups are not). It also includes uncertainties associated with inter-individual variation in factors (e.g., drinking water intake rates) that affect conversions to a set of common dose metrics. While the approach herein accounts for uncertainties due to inter-individual variation in the factors used in the dose conversions, it is essential that the values for those factors are adequately identified and supported in order to lend credibility to the analysis. The ideal situation would be that the necessary factors for dose conversion are reported in the study under consideration and hence directly apply to the population investigated in the study. For example, for dose conversions of the Chen et al. (2010) study, the average duration of well-water exposure, reported duration of well exposure, and average age of diagnosis are derived from information reported in Chen et al. (2010) or studies investigating the same cohort (Hsieh et al. 2008). However, it is not always possible to derive all the necessary factors from the considered study. In cases such as this, other resources that provide data relevant to the general study population should be considered. For example, exposure factor handbooks can be a useful resource for locating data derived for the general population from which the study population is drawn. For example, the dietary intake value for the Chen et al. (2010) dose conversion was derived from information found in the Taiwanese Exposure Factors Handbook (TDOH 2007). Correspondingly, if the study under consideration were a United States population, EPA's Exposure Factors Handbook (2011) is an invaluable source of information on body weights, dietary habits, and exposure pathways for multiple subpopulations in the United States. Studies of different cohorts, but from the same geographic area and/or ethnic population, can also be useful sources of information regarding necessary conversion factors. It is important to note that while inter-individual variation is considered, we derive uncertainty representations appropriate for the means of groups that specifically have the sample sizes of the particular studies in question.
An additional consideration in identifying studies for use in a DRMA, and hence all necessary pre-analysis conversions, is what to do with studies that report associations between specific health outcomes and chemical exposures when the dose metric is a biomarker such as blood, urine, or tissue concentrations. In cases like this, a physiologicallybased pharmacokinetic (PBPK) model would be necessary to convert those biomarker metrics into an estimate of intake (e.g., oral consumption) for dose calculations. For arsenic specifically, we identified a number of studies that reported associations between health outcomes and urinary arsenic concentrations. To convert those urinary measures of arsenic exposure in to daily intake values, we used an arsenic PBPK model (El-Masri and Kenyon 2008). As PBPK models are typically chemical-specific, we do not focus on these conversions in this paper as it is intended to be as generally useful as possible and reports on manual dose conversions that can be easily applied to other chemicals. However, functionally, the use of a PBPK model would involve converting the urinary concentrations of arsenic into corresponding model estimates of daily intake and using MC methods similar to those described above to assess inter-individual variability. As mentioned in the Methods section, converting disparate exposure metrics (e.g., drinking water concentration, urinary iAs concentrations, cumulative exposure) into a common dose metric allows for a DRMA that uses a larger set of studies. Using more studies in the DRMA allows for a more precise estimate of the association between iAs exposure and the health outcome being modeled.
There are extensions and expansions of the ideas presented here that could be applied. For example, a lognormal assumption for describing exposure variability in the populations of interest was assumed. This is consistent with the approach of Moon et al. (2017), another meta-analysis of arsenic health effects. It is also consistent with the assumption of multiple, multiplicative effects contributing to exposure, as shown above in Equations 3 and 4 (Limpert and Stahel, 2011). However, we acknowledge that it may be just a first-level approximation. Indeed, such an assumption did not really describe well the observed proportions of the sampled individuals in every exposure range. A future refinement of this approach could be to consider a mixture approach where the population consists of a mixture of subpopulations (with proportions to be estimated) each of which has a separate exposure distribution. This requires just a little more sophistication with respect to estimating means for the exposure groups that are presented in published reports. Ultimately, for our preferred Bayesian approach, we would prefer to link the estimation of exposure distributions not only to the estimation of group-specific means, but also to the logistic regression approaches described in the next installment of this series. More will be said about that in the second paper of this series, where the Bayesian aspects of our analysis are presented in detail.

Conclusions
The analyses detailed herein describe a set of pre-analysis steps that were undertaken for an inorganic arsenic case study in order to account for uncertainty when calculating group means from interval censored reported exposures, calculation of a common intake metric (μg/kg day) from multiple studies using multiple exposure metrics (μg/L, μg/L·years, etc.), and calculating effective count response data that takes into account the individual study's controlling for confounders. It is important to note again that although these methods were developed in conjunction with the Bayesian DRMAs described in a companion paper, they may be applicable to any meta-analytic or DRMA analysis in that they allow for calculation of response measures that take into account confounders (and therefore are more useful than crude or raw effect measures) and allow for the use of exposure metrics presumed to be more closely related to the response in question.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material. The analyses covered in this paper are highlighted in blue whereas the analyses highlighted in orange and green are covered in the second paper of this series (Allen et al., 2020)  2.706 (a 95% bound on the high-group mean). LL is the log-likelihood for the lognormal distribution for the candidate parameter vector; MLL is the maximum log-likelihood. When a published study reports the mean or median values for each group, those values are used directly as the group-specific dose values, with no lognormal fitting. 2 The terminology "low-end," "high-end," and "best" estimates are used to avoid confusing the values with credible (or confidence) interval bounds having a specific numerical value (e.g., 95%). Combining the log-likelihood bounds for group-specific means, with percentiles from the Monte Carlo analysis does not allow determination that the bounding estimates have any identifiable associated "confidence level." They do, however produce reasonable semiquantitative limits on how uncertain the resulting estimates are.  Plot represents the distribution of the percentiles of the dose-estimation MCMC analysis for the low-and high-exposure groups. Frequencies were calculated in Excel using the "Frequency" function and user-defined bin values.    Table 1 for variable values and justifications 2 WE (well exposure, μg/L) = CE (cumulative exposure; μg/L-years) /RDWE (yr); CE for this table is the "Best" cumulative exposure calculated for the referent group: 157 μg/L-years (see Table 4) 3 Assumed distributions with associated means and standard deviations are sampled a number of times equal to the dose-specific Ns, for doses with N > 1000, 1000 samples are taken to ease computational burden Environ Int. Author manuscript; available in PMC 2021 September 01.    a The information in the first seven columns is "borrowed" from Chen et al. (2010) if we were to treat that study as a cumulative incidence study.
The last two columns are computed as described in section 2.5.3.
Environ Int. Author manuscript; available in PMC 2021 September 01.