Data were collected from January 2005 through December 2009 from 36,928 enrolees in 85 OTPs located in 34 states. Not all programs began study participation in January of 2005; 70 (82%) OTPs participated for 2 or more years; 12 (14%) participated for at least one year but less than two years; and 3 (4%) for less than one year. All programs were federally approved opioid agonist treatment programs and followed federal methadone treatment protocols that require an opioid-dependence diagnosis and an addiction history of at least one year . The research protocol was approved by the Institutional Review Board of the National Development and Research Institutes, Inc., and oral informed consent was obtained from the study participants. The respondents in this study include treatment seeking persons who reported abusing POs or heroin in the past 30 days and were not in methadone treatment in the previous 30 days. Details of OTP participation and subject recruitment protocol can be found in an earlier publication by our group  and are only briefly reviewed here.
OTPs were selected based on our knowledge on which OTPs would most likely participate in the study and to represent regions in the U.S. where prescription opioid (PO) abuse was believed to be prevalent, e.g., non-urban areas, especially those in the Appalachian region. Some OTPs were located in major metropolitan areas such as San Francisco and New York City, where PO abuse among OTP enrolees is believed to be relatively less prevalent. A total of 7,814 respondents were not included in analysis due to missing data on age (n = 1037), age less than 16 years (n = 5), age greater than 70 years (n = 24), missing data on type of opioid use in the month before enrolment (n = 1713), or use of neither heroin nor POs in the month before enrolment (n = 5271), resulting in a sample for analysis of 29,114 OTP enrolees.
Enrolees were asked to indicate their age in whole years. For summary and analysis, seven age categories were formed: 16-19; 20-25; 26-29; 30-39; 40-49; 50-59; and 60-70. In analysis to determine associations between enrolee age and the probability of each type of opioid use, Helmert contrasts were applied to these age categories so that each age group except for the first was compared with all younger enrolees. The mean age of patients in each OTP and deviations of individual enrolee age from the OTP mean age were calculated.
United States region
OTPs were classified as falling into one of four U.S. regions. These regions, as designated by the U. S. Census Bureau , are: Northeast (Maine, New Hampshire, Vermont, Massachusetts, Rhode Island, Connecticut, New York, New Jersey, Pennsylvania); Midwest (Ohio, Indiana, Illinois, Michigan, Wisconsin, Minnesota, Iowa, Missouri, North Dakota, South Dakota, Nebraska, Kansas); South (Delaware, Maryland, District of Columbia, Virginia, West Virginia, North Carolina, South Carolina, Georgia, Florida, Kentucky, Tennessee, Alabama, Mississippi, Arkansas, Louisiana, Oklahoma, Texas); and West (Montana, Idaho, Wyoming, Colorado, New Mexico, Arizona, Utah, Nevada, Washington, Oregon, California, Alaska, Hawaii).
Beale urbanicity code
OTPs were in counties coded as high density areas (population > 1 million); moderately populated counties (≥250,000 and <1 million residents), and low populated counties (<250,000 residents); these three categories were determined by a modified version of the Beale Urbanicity Code .
Recent heroin use
Enrolees were asked whether they had used heroin in the 30 day period before admission (Yes/No).
Recent prescription opioid use
Enrolees were asked whether they had used a prescription opioid to get high in the 30 day period before admission (Yes/No).
To account for clustering of the 29,114 patients in 85 OTPs, multinomial logit mixed model analysis with a random intercept  was used to predict 1) heroin use only; 2) PO abuse only; and 3) both types of use for individual patients. The xtlogit and gllamm procedures of the Stata program [12, 13] were used to fit these mixed models. In all multinomial logit mixed models, the reference category for the dependent variable of type of opioid use was heroin only. Patient-level predictors included Helmert contrasts of the seven age groups (where each contrast compares one age cohort with all younger respondents) and contrasts of interview year (2005-2009) where the initial year of the survey, 2005, was the reference category. Program-level predictors included urbanicity, U.S. region, and the program average of the age Helmert contrasts. Including the program average of individual age group contrasts disentangles the effects of individual patient age and the effects of the average age of patients in an OTP . While the program average of the contrasts are not of direct interest, including them allows the enrolee-level age contrasts to be interpreted unambiguously as pooled within-OTP regression coefficients. This approach of breaking down the age predictor into within-OTP and between-OTP components is used to estimate the effects of individual age within OTPs . Point-biserial correlations between age and type of opioid use were estimated using the within and between analysis (WABA) approach . The WABA approach decomposes associations between age and type of opioid use into within-OTP and between-OTP components. The Mplus modelling software  was used to estimate covariances between age and type of opioid use and variances at each level. These covariances were tested for significance using z-tests and converted to correlations for presentation.