|Home | About | Journals | Submit | Contact Us | Français|
Wrote the paper: SRF BT ERP JEB HLFC HIH AL BSW. Conceived of the study: SRF. Directed the study and led the writing and analysis: SRF BT. Contributed to the analyses and interpretation of the data: SRF BT ERP CMC JEB.
People who inject drugs (PWID) have increased risk of morbidity and mortality. We update and present estimates and trends of the prevalence of current PWID and PWID subpopulations in 96 US metropolitan statistical areas (MSAs) for 1992–2007. Current estimates of PWID and PWID subpopulations will help target services and help to understand long-term health trends among PWID populations.
We calculated the number of PWID in the US annually from 1992–2007 and apportioned estimates to MSAs using multiplier methods. We used four types of data indicating drug injection to allocate national annual totals to MSAs, creating four distinct series of component estimates of PWID in each MSA and year. The four component estimates are averaged to create the best estimate of PWID for each MSA and year. We estimated PWID prevalence rates for three subpopulations defined by gender, age, and race/ethnicity. We evaluated trends using multi-level polynomial models.
PWID per 10,000 persons aged 15–64 years varied across MSAs from 31 to 345 in 1992 (median 104.4) to 34 to 324 in 2007 (median 91.5). Trend analysis indicates that this rate declined during the early period and then was relatively stable in 2002–2007. Overall prevalence rates for non-Hispanic black PWID increased in 2005 as compared to other racial/ethnic groups. Hispanic prevalence, in contrast, declined across time. Importantly, results show a worrisome trend in young PWID prevalence since HAART was initiated – the mean prevalence was 90 to 100 per 10,000 youth in 1992–1996, but increased to >120 PWID per 10,000 youth in 2006–2007.
Overall, PWID rates remained constant since 2002, but increased for two subpopulations: non-Hispanic black PWID and young PWID. Estimates of PWID are important for planning and evaluating public health programs to reduce harm among PWID and for understanding related trends in social and health outcomes.
Injection drug use continues to account for a substantial proportion of new Human immunodeficiency virus (HIV) diagnoses in the United States, and is the third most frequently reported risk factor for HIV infection, after male-to-male sexual contact and high-risk heterosexual contact , . Overall, people who inject drugs (PWID) represented 9% of new HIV infections in 2009 and 17% of those living with HIV in 2008 , . Yet, the most alarming feature of HIV among PWID in the U.S. is racial/ethnic disparities. Disparities have been apparent in HIV among PWID since early in the epidemic , , ,  and still are very marked , , , . More than 50% of PWID living with a diagnosis of HIV infection at the end of 2009 were non-Hispanic black, 27% were Hispanic, and 21% were non-Hispanic white , , . Non-Hispanic blacks who inject drugs are ten times as likely to be diagnosed with HIV as non-Hispanic white injectors , .
PWID are at high risk for HIV, hepatitis B and C, and many socially related problems. PWID experience high rates of morbidity and mortality, often from drug overdose, endocarditis, cellulitis, and abscesses , , , ,  – in addition to, increased rates of exposure to violence and injury , , . Further, PWID experience poor health outcomes due to either lack of or delayed access to effective treatment, continuation of illicit drug use, and depression and negative life events , , , , , .
Drug use patterns and prevalence reflect notable differences in the U.S. , , . U.S. figures, as reported by SAMHSA , for 2002–05 show injection was the primary route of administration for a higher proportion of heroin users (42%) than for methamphetamine (9%), stimulant (5%) or cocaine users (3%).
Although difficult to ascertain, geographic-specific data over time on PWID prevalence rates are important, as they may help policy makers allocate resources and establish public policy to reduce harm among PWID and PWID subpopulations , , . Additionally, such data could provide a foundation for the design, implementation and evaluation of structural interventions and service coverage, such as the expansion of Opioid Treatment Programs (OTPs) and Medication-Assisted Treatment (MAT) facilities in areas of need , , . Trend data on PWID prevalence may help forecast which metropolitan areas may be at greater risk for outbreaks of drug injecting and blood-borne infections associated with drug injection. Data on PWID populations also allows study of patterns of change of PWID prevalence in relation to social, economic and political predictors in metropolitan areas , .
Drug use–particularly illicit and injection drug use –carries a heavy stigma. PWID are near the bottom in terms of social tolerance in the hierarchy of client groups , , and most are reluctant to divulge any illegal drug use or needle use. Thus, despite the need for data on numbers of PWID, several factors make it difficult to assess the actual number of PWID over time and across U.S. geographic areas , , , , , , .
Brady and colleagues (2008)  created estimates of PWID population prevalence for the same 96 Metropolitan Statistical Areas (MSAs) that are studied here for 1992–2002. These data indicated an overall decreasing trend in prevalence across the 96 MSAs until 2000, after which there was a slight increase; they also identified substantial variation across MSAs in these trends .
Analyses have also discovered substantial variation in PWID trends across subpopulations over time. Cooper and colleagues (2008)  reported that non-Hispanic black PWID prevalence declined from a median of 279 injectors per 10,000 adults in 1992 to 156 injectors per 10,000 adults in 2002. PWID prevalence for non-Hispanic white adults remained relatively flat over time (median values ranged between 86 and 97 injectors per 10,000 adults). Hispanic PWID prevalence rate also declined significantly (1992 mean=192, median=133; 2002 mean=144, median=93) .
Chatterjee and colleagues (2011)  found the mean proportion of young PWID (ages 15–29) across MSAs increased by almost 20%, from a mean of 103 per 10,000 in 1996 to 122 per 10,000 in 2002. Young PWID increased significantly in nearly thirty-four separate MSAs, and declined significantly in approximately ten.
In this paper, we revise and extend the PWID estimates for the period 1992–2007 for both the overall PWID population and for specific subpopulations of PWIDs. PWID subpopulations are defined by: (1) sex (male and female); (2) age (youth 15–29 years and older 30–64 years); and (3) race/ethnicity (non-Hispanic white, non-Hispanic black, and Hispanic). Our estimates capture current injectors as opposed to ever-injected (i.e., since 1977). Measures and estimates of current injectors provide valuable information for needs assessment, public health program planning, program evaluation, and development of public policy to reduce harm among PWID. Current estimates among PWID are of particular importance for linking PWID to services and care.
In the following sections, we describe our method of estimating PWID prevalence rates for adults living in large MSAs each year between 1992 and 2007 both overall and for specific subpopulations; validate these estimates; and report temporal and spatial variations in these estimates. Lastly, we close with a discussion of the PWID estimation methods and limitations, and discuss possible causes and consequences of the observed trends in PWID and PWID subpopulations.
In this analysis, we update and amend previously published methods of estimating current PWID prevalence among adult residents of large U.S. MSAs , . For each MSA, four component sets of estimates are calculated. Each of these allocates a proportion of the national total of PWID to an MSA. The first set of estimates represents the MSAs proportion of the national total of PWID based on the proportion of HIV Counseling and Testing (Centers for Disease Control and Prevention [CDC] HIV Counseling and Testing Services database [CTS] ) events in which PWID were clients. The second uses the MSAs proportion of the national total of PWID in drug abuse treatment in the Uniform Facility Data Set (UFDS) , ,  and Treatment Entry Data System (TEDS) from the Substance Abuse and Mental Health Service Administration (SAMHSA) . These estimates both depend on the extent to which services are provided for PWID–which can vary by budgetary, siting and other decisions made by funders or service providers. They thus need to be balanced by estimates that might be higher where under-servicing is particularly high. We accomplish this by using a third component estimate reflecting: (1) the extent to which HIV-positive PWID develop new AIDS diagnoses (data from CDC’s National HIV Surveillance System on AIDS diagnoses among PWID) ; and (2) adjusting this with our previous estimates of HIV prevalence among PWID in the MSA . A fourth estimate is based on extrapolation of estimates from earlier years to more recent years , , creating four distinct series of estimates of the number of injectors in each MSA and year. The final estimates are the result of smoothing each U.S. apportioned data series over time using loess regression and taking the mean value of the four series as a best estimate of PWID for that MSA and year.
In order to compute MSA-level estimates from the national estimates of PWID for each year, we used multiplier methods to allocate the data. In the past, such methods have been widely used to estimate the numbers of global or local problem drug-using populations , , . National estimates were attributed to MSAs using the four data series discussed above (Figure 1). As in Friedman et al (2004)  and Brady et al (2008) , we hypothesized that apportioning the national number of PWID to each MSA using different data sources helps to balance biases and offset potential limitations of individual datasets. For example, increases in drug treatment funding and treatment slots for PWID could make it appear as if the number of PWID was increasing. By combining our series measures, some of which are either not subject to service bias or inversely related to services, we mitigate such service bias effects.
The resulting estimates were smoothed and averaged to create a PWID-specific estimate for each MSA and year. These MSA estimates were validated using yearly cross-sectional correlations with variables theoretically related to drug injection (i.e., unemployment, HCV mortality, and drug-poisoning mortality). Concordance between these constructs and our estimates lends credence to the validity of our estimates . Our methods and steps used to create national and MSA-specific PWID prevalence estimates from 1992 to 2007 are described below and are outlined in Figure 1.
Lastly, we estimate PWID prevalence rates for three sets of subpopulations: (1) sex (male and female); (2) age (youth ages 15–29 years, older ages 30–64); and (3) race/ethnicity (non-Hispanic white, non-Hispanic black, and Hispanic). Our subpopulation estimates use a different set of procedures than that of the overall PWID estimates. Here, to estimate subpopulation prevalence we multiplied yearly estimates of the proportion of PWID who were in the subpopulation in each MSA by yearly estimates of the total number of PWID in each MSA, and divided by their respective populations. The proportions of PWID who were in each subpopulation were estimated by combining three data series based on indicators of PWID use of health services and AIDS surveillance data. The PWID subpopulations are combined by using predicted values from a binomial mixed-effects regression to produce MSA-specific rates. The PWID subpopulation estimates are also validated using yearly cross-sectional correlations with variables theoretically related to drug injection. Figure 2 includes a flow chart of procedures for estimating PWID subpopulations prevalence.
This analysis estimates the annual number of current PWID in 96 MSAs for each year from 1992 to 2007. The June 30, 1993 Office of Management and Budget  MSA definitions, which were based on the application of 1993 metropolitan area standards to 1990 census data, were used in these analyses, except in New England, where we instead used New England County Metropolitan Areas (NECMA) definitions , .
MSAs are defined to represent a central population in a large city and surrounding urban and suburban areas which have strong economic and social ties , . The MSA is a particularly useful geographic unit for studying PWID and HIV epidemiology among PWID for several reasons. Drug distribution networks often parallel networks of distribution of legal commodities , . In addition, social networks of drug users often include individuals from inside and outside central cities, but include individuals from multiple MSAs less often , . Most importantly, MSAs are large enough to have discrete HIV outbreaks and unique data on drug use and HIV while still being small enough to be sensitive to relatively small differences in HIV outbreaks and small changes in social and economic contexts and in population characteristics over time.
We analyzed three data series (described below) to estimate the number of PWID nationwide for each year of the study period. Our first step was to create an index for each data series by dividing the number of PWID in the U.S. for that year and the data series by the average number of PWID over the study period for that series. Indices serve as a way to compare values of a variable over time by relating each value in a time series to a reference or standard value. For each year, index values were averaged across data sources to create an overall index for each year. The first set of national PWID estimates was the ratio of the overall index in a given year to that in 1992 multiplied by Holmberg’s 1992 estimate of the number of PWID . The second set of estimates was the ratio of the overall index in a given year to that in 1998, multiplied by Friedman’s 1998 estimates . The final set of national estimates was created by averaging the first and second set of approximations (Figure 1).
To create the national estimates for PWID for 1992–2007, we utilized three different data series to estimate the number of PWID nationwide from 1992–2007: (1) data from the CDC’s CTS which monitors the utilization of counseling and testing services in CDC-supported sites; , ; and (2) the SAMHSA’s TEDS  to determine the number of injectors entering treatment nationwide each year of the study period; and (3) the Uniform Crime Reporting Program County-Level Detailed Arrest and Offense time series data, produced and distributed by Inter-university Consortium for Political and Social Research (ICPSR) , to calculate the number of arrests between 1992 and 2007 for possession of opium, cocaine, and their derivatives (i.e., “hard” drug arrests). These data include arrests of non-injecting users and of PWID.
In 2005–07, CTS data reported lower national testing numbers which may be attributable to fewer health departments reporting test level data for those years based on a new reporting form starting in 2005. In 1999–2004 approximately 48–50 health departments were reporting test level data whereas those numbers have been decreasing 2005 (43), 2006 (34), 2007 (31). To correct for this, we predict US CTS data for year 2005–2007 using a quadratic polynomial model based on year and using years 1992–2004.
In creating U.S annual PWID estimates, we utilized the method of Brady et al  by first creating index values  for each of the three data series described above. We use the average value for the 16 years of our study period as our reference value for each data series. For example, the index value for TEDS, in a given year, is the number of PWID entering treatment in that year divided by the average number of PWID entering treatment 1992–2007, multiplied by 100. The indices for the three series were averaged across data sources to create an overall index for each year which was less sensitive to year-to-year fluctuation from each data source.
We then multiplied the index by (1) Holmberg’s 1992 PWID estimate  and by (2) Friedman et al’s 1998 PWID estimate , thus creating two sets of national PWID estimates for each year of the study period. The final set of national estimates was created by averaging the Holmberg-based and Friedman-based ratios to calculate a “best” national PWID estimate (Figure 3).
As in Brady et al. (2008)  and Friedman et al. (2004) , multiplier techniques were utilized to allocate national estimates of PWID to each MSA. To estimate the prevalence of PWID in each MSA, we calculated the number of PWID in each MSA and year from each of four data series as described in Table 1. We refer to these data series as the “component series”; the number of PWID from a given component series, MSA and year as “component series counts”; and MSA-level estimates derived from a given component series and the national estimates as its “component estimate.” Each component estimate required a unique methodological approach.
We use the following formula to apportion the national PWID estimate for each year to each MSA for the CTS and drug treatment-based data series independently.
Formula 1: Calculating the database-specific PWID in each MSA and year.
Ait – Component estimate for an MSA i in year t;
Bit – Component series count for an MSA i in year t per 10,000 population 15–64 years;
IUSt – Estimated number of PWID for the entire U.S. in year t;
CUSt – Component series count for the entire U.S. in year t;
RUSt – Ratio of the estimated number of U.S. injectors to the component series count for U.S. in year t.
For the first component estimate using CTS data  we applied Formula 1 to allocate the number of people tested for HIV who identified injection drug use as a risk factor in the CDC’s HIV CTS data. These testing data also included men who have sex with men (MSM) and who also reported injecting drugs. For each year we allocated our national PWID estimate to MSAs by taking the number of PWID tested for HIV per population aged 15–64 years in the MSA times 10,000 and multiplying this number by the ratio of the national number of PWID to the number of PWID tested in the U.S.
Identical methods were used to create estimates based on drug treatment. Our second component estimate uses the drug treatment data (SAMSHA’s UFDS/N-SSATS , ,  and TEDS ); we utilize Formula 1, and for each year, our national PWID estimates were allocated by taking the number of drug users (UFDS/N-SSATS) multiplied by the proportion of drug users entering treatment who injected in that MSA (TEDS) and multiplying this number by the ratio of the national number of PWID to the number of drug users in treatment in the U.S. (UFDS/N-SSATS).
Our third component series was created using data on the yearly number of incident AIDS case diagnoses ,  with transmission categories noted as injection drug use or both male-to-male sexual contact and injection drug use. We have adapted these data to use time-varying estimates of HIV prevalence rates in each MSA and in the US as a whole (Appendix S1 describes in detail our methods for estimating time-varying HIV prevalence rates in each MSA and in the US as a whole among PWID for years 1992–2008). We base our AIDS component estimates on the assumption that there is a constant proportionality of HIV prevalence rate five years before to the percent of PWID who are HIV-positive and who develop AIDS. We use a five-year lag because an untreated person takes approximately 5–10 years from HIV infection to progress to AIDS infection , , , and because using a lag larger than five would cut into the number of years for which we can make an adjustment. The assumption of constant proportionality is counterfactual but benign for this purpose. It is benign to the extent that a large part of the variation of HIV infection and progression to AIDS is due to poor access to Highly Active Antiretroviral Treatment (HAART) and utilization of HIV counseling and testing or drug treatment services for PWID–which implies that errors in this adjustment should counter-balance errors in estimates due to changing service-levels and service encounter biases. Formula 2a thus assumes that each MSA’s share of the total number of injection-related AIDS cases diagnosed nationwide in year t equals that MSA’s share of the total number of injection-related HIV cases nationwide in year t −5.
Inserting the adjustment into equation 2a, we get:
Iit – Estimated number of PWID in the MSA i in year t;
IUSt – Estimated number of PWID in the U.S. in year t;
Ait – Number of AIDS cases diagnosed in PWID in the MSA i in year t;
AUSt – Number of AIDS cases diagnosed in PWID in the U.S. in year t;
Hit5– Estimated number of HIV cases among PWID in the MSA i in year t−5;
HUSt5– Estimated number of HIV cases among PWID in the U.S. in year t−5.
Solving equation 2b for the number of PWID in an MSA in a given year t gives us equation 2c.
Lastly, to complete the AIDS adjusted component estimate, for years 1992–1996, we use data calculated from Brady et al (2008) .
Our fourth and final component estimate updates and extends the 1992–2002 Holmberg-Friedman PWID estimation series for years 2003–2007. We predict data for year 2003–2007 using a polynomial regression model using quadratic and cubic measure of time based on year and using years 1992–2002. The Holmberg-Friedman PWID 19922002 estimate series were calculated by Brady and colleagues 43 by interpolating yearly MSA-specific PWID estimates based on Holmbergs 1992 56 and Friedmans 1998 55 estimates, and then extrapolating the number of PWID in each year from 19992002 by using the yearly change in PWID between 1992 and 1998 (see Brady et al (2008) for a detailed discussion of the Holmberg-Friedman PWID estimate series calculation).
We correct for missing values in each data series when calculating the PWID component estimates. Missing values occurred for a variety of reasons, although most arose because the source database obtained by the project lacked data on either CTS HIV testing or drug treatment entry data for specific MSAs and years under study. Thus, missing data were considered under the following circumstances:
When data were missing for any of these reasons, we imputed a value for this component equal to the average of the values for the non-missing components. For this, single linear regression imputation was used to estimate the predicted value of these observations. For each component series, missing values were predicted by the other component series separately. The missing values were replaced with imputed predicted values for that MSA based on regression models results for each component series with non-missing.
The values of each of the component estimates for each MSA across time were smoothed using loess regression, which fits curves to noisy data and smoothes data in a manner similar to computing a weighted moving average on time series data.
The predicted values of the component estimates resulting from the loess regressions for each MSA and year were averaged to create our final estimates (Appendix S2). We chose a smoothing coefficient, between 0.4 and 1 by intervals of 0.1, for each MSA based on the smallest AICC1 criterion . Statistical analyses were performed by using SAS software, version 9.2 .
The following sections describe our methods for calculating the proportion of PWID in three subpopulations across MSAs and year. We stratify each three subpopulation into specific groups: (1) sex (male and female); (2) age (youth 15–29 years and older 30–64 years); (3) and race/ethnicity (non-Hispanic white, non-Hispanic black, and Hispanic). Figure 2 describes our methods and procedures for calculating the subpopulation prevalence estimates. We further assess the reliability and validity of both the overall PWID estimates and PWID subpopulation estimates, and discuss trend results for both.
Our estimates of MSA-specific subpopulation prevalence for PWID are built on previous methods and estimates by Cooper et al, (2008) , Pouget et al (2012) , and Chatterjee et al, (2011)  although no previous prevalence estimates exist for men and women who inject drugs. The subpopulation prevalence estimates utilize our current PWID estimates as described above. The proportions of PWID who were in each subpopulation were estimated by combining three data series: (1) CDC CTS ; (2) SAMSHA TEDS ; and (3) CDC’s data on PWID diagnosed with AIDS . The CTS data had missing values for four of the 96 MSAs for the entire study period (Honolulu, HI, Little Rock–North Little Rock, AR, San Juan-Caguas-Arecibo, PR, and Springfield, MA). We also excluded existing data for one to six years from 14 MSAs in CTS in order to avoid proportion estimates that could be unreliable due to small total numbers of PWID (fewer than 20). In addition, data on proportions of each subpopulation were available from at least 86 MSAs in each year from the TEDS treatment data series. Lastly, AIDS data were not reported according to racial/ethnic category in some states. This included two MSAs (Cincinnati, OH–KY–IN and Louisville, KY–IN) within our study boundaries.
To reduce potential bias due to variation among MSAs in relative HIV prevalence, we adjusted the proportion of PWID who were in each subpopulation from the AIDS data (PAIDS) for HIV prevalence among PWID calculated from CTS data (see Appendix S1 for adjustment formula). Database-specific estimates were combined by using predicted values from a binomial mixed-effects regression to produce MSA-specific rates.
To estimate subpopulation PWID prevalence we multiplied yearly estimates of the proportion of PWID who were in the subpopulation (averaging proportions from three data series) in each MSA by yearly estimates of the total number of PWID in each MSA, and divided by their respective populations. The age range used for population subpopulation denominators was 15–64, except for estimates of young PWID, where the 15–29 age range is applied. The calculation is summarized below.
Subpopulation ij=Estimated prevalence of PWID among subpopulation residents in study year i, MSA j in database k.
PropSubpopij=Estimated proportion of PWID who were in the subpopulation in study year i, MSA j in database k.
PWIDNij=Estimated total number of PWID, aged 15–64, in study year i, MSA j in database k.
PopSubpopij=Number of residents aged 15–64 who were in the subpopulation, in study year i, MSA j and database k.
Missing data for calculating subpopulations estimates were imputed separately. CTS data on the proportion testing positive for HIV among PWID subpopulations were incomplete due to inconsistent reporting, removal by the CDC of testing results from small cells (fewer than 5 positive results) to protect participant confidentiality, and exclusion of data on the proportion positive that would have been based on fewer than 20 PWID.
Missing CTS data for computing the proportion testing positive for HIV among PWID subpopulations were imputed in two steps. First, in a binomial mixed-effects regression, values were imputed as a function of: (a) cubic polynomial terms for time; (b) the proportion of total PWID testing positive for HIV; and (c) the proportion of PWID tested who were in each subpopulation.
Mixed-effects regressions were performed using the SAS procedure PROC GLIMMIX . Intercepts and time parameters were set to vary randomly using residual pseudo-likelihood estimation. Next, in MSAs where PWID subpopulation HIV data were missing, we imputed missing values using predicted values from a mixed-effects Poisson regression using a cubic polynomial model of time. We imputed values to fill in missing CTS HIV data for total PWID in a parallel manner, using year, the proportion of all CTS clients testing positive for HIV, and the proportion of CTS clients tested who were PWID as predictors for the mixed-effects binomial regression. Thirteen MSAs had missing data for all study years; as a result we replaced missing cells with the average values as predicted by the binomial regressions.
We created a complete data set of the proportions of PWID who were in each subpopulation from the three data sources–CTS, TEDS and AIDS (adjusted for HIV prevalence)–using predicted values from a binomial mixed-effects regression, with one exception: in the San Juan-Caguas-Arecibo, PR MSA we assume that all of the PWID were Hispanic because CTS data were not available, TEDS data were not collected by race/ethnicity, and AIDS data indicated that few PWID with AIDS were non-Hispanic).
Data were stacked so that each MSA and year was represented by as many observations as there were data sources (1, 2, or all 3 sources). The model used all the available data to estimate parameters, even where some data were missing, under the assumption that the data were missing at random, conditional on the observed data. Dummy codes were used to compare data sources. We used a cubic polynomial model of time to maintain consistency with procedures used to produce the overall PWID estimates. Intercepts and time parameters were set to vary randomly using residual pseudo-likelihood estimation. The resulting complete sets of proportions were averaged to create single yearly estimates of the proportion of PWID who were in each subpopulation. These proportions were then multiplied by estimates of the number of total PWID to produce estimates of the PWID subpopulations. The resulting numbers of PWID in each subpopulation were divided by each subpopulation, respectively, to produce subpopulation prevalence rates.
To test the reliability of our overall PWID estimates (i.e., the extent to which our component estimates seem to be capturing the same underlying construct), we examined the correlation between each of our component estimates for each year. These correlations describe the extent to which our component estimates produce similar results.
To assess the reliability of our subpopulation estimates were compared the proportions of PWID in each subpopulation across sources of data. We used CTS, TEDS, and AIDS data for racial/ethnic subpopulation estimates. For estimating sex and age PWID subpopulations we used CTS and TEDS data only because AIDS data showed less consistency and produced weak correlations with the subpopulation proportions in CTS and TEDS data. As such, we considered the AIDS data less reliable as a measure for these specific subpopulations estimates. By definition, young PWID have less time to become HIV-infected and progress to AIDS, and thus we expect low correlations because AIDS data among youth are sparse, reflecting less time-at-risk to develop AIDS among youth. Additionally, HIV is generally found to be lower among females, which reflect less stable numbers of AIDS diagnosis , , .
To assess the validity of the overall PWID prevalence estimates, and the PWID subpopulation estimates we compared these with mortality related to drug-poisoning and HCV mortality. A minimum of three deaths was used as a criterion for calculating subpopulation mortality rates. We do not present HCV mortality analyses for young PWID since they would have had less time to develop fatal liver disease caused by HCV acquired by injecting drugs. Additionally, numerous studies have found that economic conditions are associated with rates of substance use and/or HIV prevalence , , , , . Thus, we further utilize a number of socioeconomic metropolitan area characteristics, including income inequality, related to the subsequent prevalence of drug injection in the population.
We include drug-poisoning mortality as a key validator due to the high rates of overdose deaths that occur among drug-using populations. For example, in the U.S., drug overdose is the second leading cause of unintentional injury death behind motor vehicle crashes, but is the leading cause of injury death among persons 35 to 54 years of age . We refer to drug-poisoning mortality as those deaths that are directly caused by an accidental or unintentional drug-poisoning death. Since injection drug use cannot be determined in the multiple causes of death files , , , , for both ICD-9  and ICD-10 , ,  coding, we included a range of substances that correspond to injection practices. For the ICD-9, we had included a more expansive list of substances; however, for ICD-10 coding, we decided to adapt our coding scheme to include only opioids (heroin, methadone, opium, other opioids, and other synthetic narcotics), cocaine and psychostimulants. The motivation for this change was both a desire to simplify the coding scheme and to reduce inconsistencies in case definitions inherent in trying to identify overdose deaths related to injection drug use. Drug-poisoning deaths were included if the MSA in which the death occurred matched one of our study MSAs. Our algorithm for drug-related mortality is adapted from the European Monitoring Centre for Drugs and Drug Addiction (EMCDDA) ICD-9 and ICD-10 coding scheme , .
HCV virus (HCV) was selected as a validator because it is a blood-borne infection most commonly transmitted in the U.S. via PWID . The HCV was first identified in 1989 and widespread screening and testing of the virus began in 1993. Thus the ICD-9 coding conventions used to identify HCV might result in an underestimation of the contribution of HCV to mortality in the early years of the study period , , , . We utilize the ICD-10 coding , ,  scheme to identify HCV mortality, defined as deaths recorded as either acute (B171) or chronic (B182) HCV ; accordingly the correlations values only reflect data from 1999 to 2007.
We analyze trends over time for overall PWID prevalence rates and subpopulation prevalence rates by utilizing mixed-effects polynomial models with PROC MIXED using linear, quadratic and cubic terms for time . We centered time at the mid-point of the study period, year 2000. Nested models were compared using likelihood ratio tests .
This is a study of 96 MSAs that are the largest metropolitan areas in the U.S. Therefore, it was a study of a population rather than of a sample. This means that there was no sampling error (although there was measurement error). Statistical analyses are primarily descriptive rather than inferential. Thus, the use of p values in this analysis is a guide to the relative “importance”, and not the significance of possible associations. Some researchers conducting studies with similar populations use P values as a heuristic device to avoid over-interpreting model parameters , ,  (We refer to these as pseudo-Ps.) Other analysts would consider the population to be a random sample of “possible universes”; under this interpretation, pseudo-Ps have a probabilistic interpretation.
We assess the reliability of our PWID estimates by comparing each component series estimate across data sources. Figure 4 presents the correlations across MSAs of the component series estimates by year for the overall PWID prevalence estimates. In general, the intercorrelations were moderately strong and mostly stable over time, which is expected because the component estimates aim to measure the same underlying construct. As expected CTS and TEDS component series correlate more strongly with each other than with AIDS component series, since counseling and testing, and drug treatment represent services that can help to prevent HIV infection or delay progression to AIDS for those who are infected. These was an overall declining correlation between the drug treatment counseling and testing and AIDS-based component estimates, and to some extent, between the counseling and testing and AIDS-based component estimates which may suggest that the performance of the AIDS adjustment formula that relates HIV prevalence to AIDS cases deteriorates over time, or that our estimates are affected by service bias to a greater extent over time.
Convergent validity was assessed by evaluating if the PWID prevalence estimates correlated with theoretically relevant predictors, including unemployment, drug-poisoning mortality and mortality from HCV. Table 2 displays the cross-sectional correlations between PWID prevalence and theoretical constructs. Correlations between PWID prevalence and unemployment varied little over time ranging from 0.42 in 1992 to 0.35 in 2003 where correlations declined starting 2004. Correlations between PWID prevalence and drug- poisoning mortality held steady over time ranging from 0.51 in 1992 to 0.40 in 2003– correlations declined starting 2004. Correlations between PWID prevalence and HCV mortality were consistent across years, ranging from 0.48 in 1999 to 0.43 in 2007.
Additionally, correlations of the estimate of the prevalence of PWID with theoretically related indicators, such as percent living in poverty and income inequality were analyzed. The correlations between percent in poverty in 1990 and our PWID estimates per 10,000 population in 1992; between poverty in 2000 and our PWID estimates per 10,000 population in 2000 were 0.33 and 0.36 respectively; and finally between income inequality 1990 and 2000 and our 1992 and 2000 PWID estimates were 0.24 and 0.23, respectively (each of these is significant at p=0.05 level). This suggests that for the given years our estimates were similarly accurate to those estimates presented by Brady and colleagues (2008) .
Pearson correlations values for the subpopulation estimates are presented in Table 3. Reliability analysis on the subpopulation proportions among the data series showed good consistency for most comparisons as well. CTS and TEDS proportions were moderately correlated for young PWID and for male PWID, and correlated somewhat less for female PWID. For all three racial/ethnic groups, CTS and TEDS proportion data correlated more strongly with each other than with AIDS proportions.
Female and male PWID proportions in CTS and TEDS data may not as accurately represent the proportion of females and males among PWID in the MSA populations relative to proportions for the other subpopulations. Weak correlations found for female PWID subpopulation may reflect the fact that female PWID are much more likely to be stigmatized by society than their male counterparts because drug use violates social norms of behavior as it diverges from the traditional expectations of women as nurturers of families [86, 87). As such, female PWID may be more inclined to conceal injecting behavior and hesitant to utilize counseling and testing or drug treatment services , , , , .
Subpopulation validity was assessed by comparing subpopulation PWID prevalence rates with rates of drug-poisoning deaths and HCV deaths series showed moderate to strong associations for most comparisons (Table 4). Correlations with HCV deaths begin in 1999, when ICD-10 implemented explicit coding for HCV. Across subpopulations, correlations between PWID rates and drug-poisoning deaths were higher in earlier years. Correlations averaged 0.52 for drug-poisoning deaths from 1992–1999, and 0.35 for 2000–2007. Correlations with HCV deaths were 0.30 or above for all subpopulations in all years, averaging 0.50 across subpopulations. Correlations with PWID subpopulations and drug-poisoning deaths generally were reduced over time. Some values were below 0.20 for Hispanic and male subpopulations in some years. We conjecture that the downward trend in these correlations may reflect an increase in the contribution of non-injecting opioid analgesic users to drug-poisoning deaths over the study period .
Appendix S2 shows the estimates of the number of PWID per 10,000 population aged 15–64 years for each of the 96 largest MSAs in the US for each year from 1992 to 2007. Mean PWID prevalence rates and standard deviations are presented in Table 5. Overall, the number of PWID per 10,000 persons aged 15–64 years varied from 31 to 345 across MSAs, median 104.4 (mean 127.4; standard deviation 66.7; percentile range 76–162) in 1992 and from 34 to 324 across MSAs, median 91.5 (mean 103.6; standard deviation 56.4; percentile range 61–125 ) in 2007 indicating an overall decline in PWID prevalence across MSAs.
Figure 5 shows the overall trajectory of the PWID prevalence rates based on the multilevel model. Trend analysis of the overall results is consistent with a decline in the early study period, followed by an increase in 2000–02, and then remaining stable thereafter over time. On average there has been very little change since 2002 (mean 105.0) to 2007 (mean 103.6). Overall, across the 96 MSAs the mean PWID prevalence mostly decreased during our study period, as did the dispersion of estimates over time.
Trend analysis results for the overall PWID prevalence are in Table 6. The linear model indicates the average number of PWID across MSAs was about 110 PWID per 10,000 population age 15–64 years in year 2000 where each year from 2000, the average number of PWID across MSAs declined slightly by about 1.4 per 10,000 population. Quadratic and cubic terms are also presented. Overall linear and quadratic terms for PWID were significant; the somewhat smaller linear decline was offset over time by the quadratic term. Although the cubic term in Model 3 was not significant, it significantly improved the fit of the model. All models indicate a decline in PWID prevalence, increasing in years 2000–02; thereafter, the trend indicates our PWID prevalence rates stabilize after 2002 and vary little through the end of the study period.
We calculated percent change values of the earliest (1992–1994) and latest (2005–2007) 3 years of data. MSAs that experienced increased rates in PWID prevalence reside mainly in the northeast (Baltimore, MD; Springfield, MA; Boston, MA–NH; Providence–Fall River–Warwick, RI–MA; Harrisburg–Lebanon–Carlisle, PA). Sarasota–Bradenton, FL; Stockton–Lodi, CA; Youngstown–Warren, OH; Salt Lake City–Ogden, UT; Tulsa, OK MSAs also experienced increases in PWID prevalence rates over time.
Table 7 shows the results of test for trends over time for each subpopulation and are also graphically displayed in Figure 5. The intercept column represents the modeled subpopulation PWID prevalence rate in the year 2000. All intercepts were significantly different from zero. For non-Hispanic white, only the quadratic term for time was significant, in the positive direction, reflecting a small increase since 2000 (7.5% in total, or an average of 1.1% per year). For the non-Hispanic black subpopulation, the linear and quadratic terms were significant; the strong linear decline seen in the early 1990s was offset over time by the smaller quadratic term, reflecting an increase in each year since 2004 (0.6% in 2005, 2.1% in 2006 and 4.8% in 2007). For the Hispanic subpopulation, linear and quadratic terms were also significant; the somewhat smaller linear decline was also offset over time by the quadratic term, averaging a yearly decline of 4.7% in 1993–1997, 2.6% in 1998–2002, and 1.5% in 2003–2007. For the male subgroup, linear and quadratic terms were also significant, with a more modest linear decline offset by the quadratic term, reflecting an average yearly decline of 2.6% in 1992–2000, and 0.5% in 2001–2007. For the female subpopulation, no time effects were significant.
Overall, the number of PWID per 10,000 persons aged 15–64 years varied from 31 to 345 across MSAs in 1992 (mean 127.4) and from 34 to 324 across MSAs in 2007 (mean 103.7). Trend analysis indicates that declines during the early period (1992–2002) leveled off, and rates were relatively stable in 2003–2007. Mean rates for non-Hispanic black PWID were the highest of all racial/ethnic groups in 1992 (337.0). After an early steep decline in the mean, non-Hispanic black PWID prevalence leveled off in 2000, and on average rates increased after 2005. Rates for non-Hispanic whites increased on average since 2003, although the rate was only about half that of non-Hispanic blacks. Hispanic PWID prevalence, in contrast, has consistently declined, although the decline has slowed across time. Mean male PWID prevalence shows a decline that decreased over time. Mean female prevalence shows no clear trend over time, although mean values since 2000 are somewhat lower than those before 2000.
Most notably, our results indicate an overall increasing trend in young PWID. For the young subpopulation, linear, quadratic and cubic terms were significant; positive linear and quadratic increases were offset with a negative cubic term over time, reflecting a general increase since 1996. After decreasing on average by 11.0% from 1992–1996 (2.8% per year), the average rate of young PWID increased by 40.0% from 1996–2006 (3.1% per year). In contrast, for older PWID, only the linear term was significant, reflecting a continuous linear decline.
Our national estimates indicate that PWID prevalence remained relatively stable from 1992 to 2007 nationwide. The overall MSA-level population prevalence of PWID decreased through the 1990s (as reported in Brady (2008) ). Trend analysis showed an overall decline in the early study period, followed by an increase in 2000–02; thereafter PWID prevalence rates remained stable. Trends among subpopulation prevalence rates varied over time. Importantly, results show a worrisome trend in young PWID prevalence since HAART was initiated – the average increase across MSAs for 1996 to 2007 was 35%.
Our validation test showed relatively high to modest correlations between each series of MSA-level component estimates across years, suggesting that the component estimates measure the same underlying construct and that our data series support each other (Figure 4). Multiple sources of data on PWID showed consistency, and comparisons with data on drug-poisoning mortality and HCV mortality supported the validity of the estimates.
PWID prevalence among MSA residents aged 30 and older closely matched PWID prevalence overall. However, prevalence was slightly higher among older MSA residents from 1992–2002, and lower from 2003–2007, reflecting a decreasing average age of older PWID. Overall, declining PWID prevalence rates among those 30 and above may stem from increases in overdose and/or HCV mortality and from age-associated mortality in aging cohorts of PWIDs. Another explanation for the downward trend in PWID prevalence in metropolitan areas is that drug users have been employing alternative modes of drug administration, or discontinue drug use in general. There are numerous reasons to choose modes of administration other than injection and discontinue use, including fear of HIV or HCV.
This might be balanced in some MSAs due to declines in HIV-related mortality due to improved treatments and decreases in HIV-associated mortality in MSAs where HIV among PWID is prevalent. Drug treatment programs or other processes may have led some PWIDs to stop drug use or at least injecting. In addition, the rising use of methamphetamine and prescription drugs, which are often administered by means other than injection, may mark a shift away from injection drug use. Increases in non-injection opioid use may have influenced changes in PWID prevalence in contradictory ways, perhaps lead some youth not to inject their drugs but leading others to start injecting .
Prevalence rates for young PWID are consistent with previously published estimates , and confirm and extend the previously observed increase that began in 1996. There are at least two potential explanations for this increase. First, young MSA residents may have avoided injecting in the early 1990s due to the well-publicized risks of HIV/AIDS, and morbidity and mortality associated with injecting drugs. Potentially, this reluctance may have been reduced over time since the HAART era due to HIV treatment optimism, similar to explanations some have hypothesized for increased sexual risk behavior among MSM , .
Moreover, recent increases in non-medical use of prescription opioids may have led some new opioid users to inject prescription opioids or transition to heroin or methamphetamine injecting , . A recent study by Jones (In Press) , found that over 77% of people using both prescription opioid and heroin in the past year report initiating a prescription opioid prior to heroin initiation. Among the three age groups examined in the study, 18–25 year olds had the highest rate of past year heroin use – of those in this age category 12.3% also used heroin in the past year, a 95% increase compared to 2002–2004. Additional studies by Church et al., (2011)  and Stanley et al., (2012)  reportedly linked new HCV infections in young adults to the use of prescription opioids and transition to heroin or methamphetamine injection.
If young PWID engage in high levels of risk behavior, this could lead to widespread increases in HIV transmission that could parallel recent increases in HIV incidence among young MSM. Young MSM increasingly assume the heaviest burden of HIV infection, presenting more than half of new infection in 2010 , . A critical component of risk among young MSM is the intersection of the HIV epidemic, substance use (which includes both injection drug use risk and non-injection risks), and unprotected sexual activity with multiple partners , , . The increasing practice of injecting methamphetamine among MSM leads to greater risk of HIV, STIs and hepatitis, offering new challenges for prevention and reducing harm among MSM injectors , , .
An early steep decline in non-Hispanic black PWID prevalence leveled off, and on average, rates increased since 2005. Mean rates for non-Hispanic black PWID were the highest in 1992 (mean 337.0), and lowest in 2004 (mean 195.0). Although representing less than 14% of the total population in the US, non-Hispanic blacks are disproportionate affected by HIV–accounting for 44% of HIV infection for 2010 , . Recent increases in non-Hispanic black PWID prevalence warrant further investigation, given that the average annual rate of new HIV infection diagnosis per 100,000 population during 2004–2007 was 11.0 for non-Hispanic black PWID, 4.9 per 100,000 for Hispanic PWID, and 0.9 per 100,000 for non-Hispanic white PWID , . As such, HIV infection among PWID subpopulations continues to be a public health issue – where non-Hispanic blacks are disproportionately affected by HIV , .
Knowledge of these findings can assist research that investigates the causes and consequences of increasing PWID prevalence among subpopulations, including how risk exposures differ among subpopulations. For example, young injectors may be less aware of the dangers of injecting drugs and how to reduce their risk, and as such, more likely to share syringes and drug preparation equipment increasing their risk of HIV, hepatitis B and C, as well as fatal overdose and homicide , , , , , , . In contrast, non-Hispanic black injectors are more likely than other PWID subpopulations to experience drug-related arrests and incarceration , . Fear of arrest has been associated with higher rates of risk behaviors, with difficulty in using prevention programs, and injection-related infections , , , , .
Resource availability and the allocation of resources can vary greatly, affecting where and when prevention and treatment services are implemented. Access to a broad range of comprehensive health services, medical care and drug treatment services, such as those offered at MAT facilities, are fundamental to achieving and reducing harm among PWID , . Effectiveness of MATs may have played a role in the overall decreases in PWID rates we observed and decreasing HIV transmission rates among PWID. Provisions of MAT for PWID and PWID subpopulations remain central to efforts for the prevention and management of HIV infections and health care among PWID; however, access and coverage of MAT might not be sufficient in some U.S. localities, further impacting PWID rates and HIV and HCV transmission, particularly among young PWID and non-Hispanic black PWID.
One of the most practical applications of our PWID estimates in terms of targeting primary and secondary HIV prevention initiatives would be the computation of population-based HIV and HVC transmission rates among PWID and subpopulations. Furthermore, trend data on PWID prevalence and subpopulations may provide a foundation for the design, implementation and evaluation of structural interventions and service coverage, such as the expansion of OTPs and MAT facilities in areas of need , , . Efforts to characterize PWID populations and PWID subpopulations have been complicated by prejudices about the stigmatizing nature of injection drug use , . Prejudice about drug use, and in particular injection drug use, makes it difficult to assess the actual numbers of PWID over time and across geographic areas , , , , , , . Our estimates on overall PWID population prevalence and PWID subpopulation prevalence have attempted to correct for this inability to measure and estimate PWID populations.
Data sources were chosen so that biases of the different data sets, in theory, would cancel one another out , . For example, the CTS data and the TEDS data both reflect health service provision: if funding for health service declines, so will the number of PWID noted for that data set. In contrast, a high component estimate based on PWID being diagnosed with AIDS, in conjunction with low service-based component estimates, may reflect inadequate services for PWID, such as having poor access to HAART. It is important that our various component estimates are not subject to the same biases, so that when they are averaged, even if one component estimate declines because of funding cuts, another component estimate will be unaffected by this event or will be affected in the opposite way.
Despite the fact that using different data sources helps to balance biases, each data series has strengths and weaknesses and are subject to several limitations. These limitations have been discussed in detail elsewhere , , , , , and are briefly reviewed here.
First, with regard to TEDS and UFDS/N-SSATS, during our study period SAMSHA eliminated questions from UFDS about the number of PWID in treatment programs after 1998. We therefore multiplied the proportion of drug users entering treatment who inject drugs (from TEDS) in each MSA and year by the total number of drug users in treatment as reported by both UFDS/N-SSATS. Second, these data sets differ in what they count: TEDS counts each admission in a given year, so an individual who enters drug treatment twice or more in a year is counted as two or more independent cases. In contrast, UFDS/N-SSATS client count numbers are from a point-prevalence survey–those in treatment on one specific day. It gives a snapshot of the substance abuse treatment system on a given day. Consequently, if PWID differ from non-PWID in the ratio of admissions to those remaining in treatment, our estimates will be biased.
Additionally, UFDS/N-SSATS data sets were not available for 1992, 1994, 1999, and 2001. We did not impute UFDS/N-SSATS for the years 1992, 1994, 1999, and 2001 because we did not believe there was a need to estimate the variance of the individual series, as the treatment component estimate would be smoothed and averaged with the other component estimates. As such, our final PWID estimates are calculated without using UFDS/N-SSATS data in the above years.
The counseling and testing data series has several limitations. First, data only exist where there are publicly funded CDC HIV counseling and testing sites. Since CTS systems only collect data from subsets of United States locations where publicly funded CDC HIV counseling and testing take place our estimates may be biased to the extent that the demographic distributions of PWIDs using these sites differ from those at sites that were not included. For example, the HIV CTS data as a whole are considered to only represent approximately 10–15% of all HIV tests in the U.S. , . Therefore, HIV testing of PWID occurring at a health maintenance organization or at a blood donation center is not included in these data and may lead to an underestimation of the number of PWID. Second, reluctance to identify and report illicit and otherwise stigmatized behaviors may lead to misclassification and underreporting of PWID in the counseling and testing data.
Additionally, the counseling and testing data classify injectors as anyone injecting since 1978, although definitions used to classify injectors may vary by health department e.g., some people were classified as PWID if they reported injecting in their lifetimes, while others were classified as an PWID only if they reported injecting in the past 30 days, or since the last HIV test , . This means we may not be capturing recent injectors in this PWID prevalence measure and that comparability across health department jurisdictions may be limited.
Included in our PWID estimate model was an adjustment for time to the advent of highly active antiretroviral therapies HAART in 1996. Since then, HAART has become the standard of care for the treatment of HIV-infected individuals. As HAART has become more widely used, AIDS researchers have seen an increase in the time from HIV infection to progression of AIDS, and an increase in survival time from a diagnosis of AIDS to death , , , . Thus, as survival time increases prevalence rates in the research estimates we needed to adjust for this in our equation. During our study period, access to HAART changed over time and took time to reach PWID–thus, we assume year of 1996 as a uniform time for which PWID had access to HAART. To adjust for this in our AIDS component estimation series, we assume that PWID in all MSAs have equal access to HAART, starting 1996.
The component series based on the interpolation and extrapolation of Holmberg and Friedman’s estimates assumes a linear change on the log scale for each MSA. The weakness in this assumption is that the changes in number of PWID may not have changed linearly (on the log scale) over time in each MSA. Additionally, this component series is further subject to all the limitations, biases, and weaknesses that Holmberg and Friedman highlighted in their papers.
Finally, our PWID estimates do not extend beyond 2007. Completion of our estimates is based on the availability and release of secondary data; as such data for producing PWID estimates beyond 2007 were not readily available. In this paper, we presented estimates of the prevalence of people who injected in the last year. As such, these estimates do not provide information about the frequency of injecting or about risk behaviors among injectors.
With proper acknowledgement of their limitations, these data can be used for structural analyses of the correlates and predictors of the population density of PWID and subpopulations in metropolitan areas. Specifically, the correlations found between the PWID estimates and theoretically relevant variables suggest that further research should be done to determine whether socioeconomic and health policy affect the population density of PWID.
These PWID prevalence and subpopulation prevalence estimates can help local agencies plan local drug-related interventions and policies. For example, where both our broader estimates and local data suggest that PWID prevalence is steady or rising, local health departments working together with syringe exchange programs might intensify outreach efforts targeting this population. Our estimates can also be used to measure the extent to which health systems and social service systems are providing services to PWID and PWID subpopulations.
For example, before now no adequate data have existed on patterns of female PWID prevalence over time within geographic areas. The lack of such prevalence data makes it difficult to assess whether at a metropolitan area has gender disparities in access to essential services , . Thus the sex-specific PWID prevalence estimates reported here can inform planning efforts for health services for female PWIDs, and help design comprehensive health programs for female PWID. In addition, in the U.S., overdose is the most frequent cause of mortality among PWID , ,  and the second leading cause of accidental injury death overall . Our PWID prevalence estimates and subpopulation estimates can deepen our understanding of the causes of determinants and changes over time in overdoses deaths and HCV mortality among PWID populations. Finally, with more young people injecting, drug-related morbidity, including overdose, HIV and HCV infection, is likely to increase in this age group. Funding for harm reduction and other health-related services targeting young PWID should be increased.
In summary, it is clearly of very high importance to conduct research to understand what it is about some, but not other, MSAs that leads to increasing PWID prevalence, and to monitor future changes in injection by youth and, in general, changes in PWID prevalence and subpopulation prevalence rates. Our estimates can be used to inform public health officials about population vulnerability to future blood-borne epidemics, and will help local policy makers differentiate why PWID prevalence and subpopulation prevalence varies over time. Additionally, our estimates may be used in research into how socioeconomic conditions, public health programs, policies, and other socio-political factors shape the prevalence of PWID in a population.
A1. AIDS cases as adjusted for HIV prevalence Figure A1: Estimation of Annual U.S. HIV Prevalence Rates among People Who Inject Drugs, Observed and Predicted Values A2. Adjustment of AIDS diagnoses for the proportion of PWID testing positive for HIV.
The Estimates of the Number of PWID per 10,000 population aged 15–64 years for each of the 96 largest MSAs in the US for each year from 1992 to 2007.
Disclaimers: The findings and conclusions in this paper are those of the authors and do not necessarily represent the views of the Centers for Disease Control and Prevention.
The research described in this paper is supported by the National Institute of Drug Abuse grant # R01 DA13336. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.