Assessing Opioid Use Disorder Treatments in Trials Subject to Non-Adherence via a Functional Generalized Linear Mixed-Effects Model

The opioid crisis in the United States poses a major threat to public health due to psychiatric and infectious disease comorbidities and death due to opioid use disorder (OUD). OUD is characterized by patterns of opioid misuse leading to persistent heavy use and overdose. The standard of care for treatment of OUD is medication-assisted treatment, in combination with behavioral therapy. Medications for opioid use disorder have been shown to improve OUD outcomes, including reduction and prevention of overdose. However, understanding the effectiveness of such medications has been limited due to non-adherence to assigned dose levels by study patients. To overcome this challenge, herein we develop a model that views dose history as a time-varying covariate. Proceeding in this fashion allows the model to estimate dose effect while accounting for lapses in adherence. The proposed model is used to conduct a secondary analysis of data collected from six efficacy and safety trials of buprenorphine maintenance treatment. This analysis provides further insight into the time-dependent treatment effects of buprenorphine and how different dose adherence patterns relate to risk of opioid use.


Introduction
There have been nearly 500,000 overdose deaths from opioids in the United States alone in the last 20 years, with associated annual costs exceeding $1 trillion [1]. The treatment of opioid use disorder (OUD) is inherently complex, with clinician assessment of the patient, comorbidities, suitability for one of the three FDA-approved medications, psychosocial counseling and care for comorbidities [2]. One of the two more heavily utilized medications in OUD treatment is buprenorphine, an opioid partial agonist. Studies have been conducted to assess the effectiveness of various formulations and doses of buprenorphine on detoxification, retention in treatment, and on the elimination of illicit opioid use [3]. A meta analysis of clinical trials found that any dose over 2 mg of buprenorphine was useful at retaining patients in treatment but only higher doses (16 mg or more) reduced illicit opioid use [4]. It has been shown that illicit opioid use and the number of weeks abstinent from illicit opioid use are significantly associated with daily buprenorphine adherence [5]. Given the potential interaction between buprenorphine dosing and adherence, further investigations aimed at better understanding these interactions are warranted and will support translational clinical research that seeks to optimize the overall effectiveness of medications for OUD treatment (MOUD).
A challenge to being able to accurately assess the effectiveness of MOUD in clinical trials is non-adherence. For example, a multicenter, randomized clinical trial (CSP-999) considered the effectiveness and safety of four buprenorphine dose levels (1, 4, 8, 16, or 32 mg/day), which were administered daily in clinic. Due to the mode of delivery, adherence for this study was directly observed, i.e., patients were either present or absent from the clinic visit. Figure 1 provides a depiction of dose history over a 50 day period for four randomly selected CSP-999 trial patients. In particular, three of the selected patients were assigned to a 16 mg/day dose while the remaining patient was assigned to a 8 mg/day dose. Days when patients missed a dose (i.e., were non-adherent) are represented by a dose level of 0 mg/day. Induction and re-induction after a lapse in dosing can be seen by increasing dose levels from 0 mg/day to the assigned levels. Failing to account for the patterns in adherence depicted in Figure 1 when trying to relate assigned dose to opioid use will obscure the true effect of buprenorphine on illicit opioid use. In particular, not accounting for these patterns will lead to underestimation of the effect (or log-odds) of dose on reducing illicit opioid use. Given the challenge of non-adherence in substance abuse trials aimed at assessing efficacy of various treatments and dosing protocols, the goals of this paper are two-fold. First, we seek to develop and demonstrate a methodology that can be used to directly acknowledge and account for the effects of adherence when assessing treatment effects. Second, we seek to use our proposed model to assess the effectiveness of buprenorphine as a MOUD. To this end, we compiled, aligned, and harmonized publicly available data from six efficacy and safety clinical trials of buprenorphine maintenance treatment with detailed logs of patient buprenorphine dose. For more on the combined data set and merging steps see [6]. As a part of these trials, patients were expected to attend weekly follow-up clinic visits with opiate urinalysis testing. Dose administration varied across the six trials, with CSP-999 being the only trial requiring daily doses to be self-administered in clinic. Weekly dose adherence for the remaining trials was reconstructed using other available information, e.g., self reported non-adherence and returned pills. Additionally available were a collection of various sociodemographic and substance use variables that were included in the analysis to address potential confounding.
To analyze these data, we develop a generalized linear functional mixed-effects model. The proposed model views daily dose level as a functional covariate whose value reflects the mg/day dose taken, with a value of 0 corresponding to days when the subject is non-adherent. By construction, our model has several key features. First, we can extract an estimate of, and conduct inference about, the dose effect for individuals with strict adherence (i.e., 100% compliance with the prescribed dosing protocol). Second, we can assess the effect of different types of dose self-administration or medication adherence patterns on illicit opioid use. Our model makes use of random effects to account for across trial heterogeneity, and across subject heterogeneity within studies. To complete model fitting, we cast our model into the Bayesian paradigm and develop a custom Markov chain Monte Carlo (MCMC) posterior sampling algorithm.

Generalized Linear Functional Mixed-Effects Model
In what follows, we outline the salient features of our proposed model, which was designed to relate illicit opioid use to time-varying dose adherence while controlling for various demographic and drug use history factors purported to be related to the same. To this end, let Y ij , for j = 1, . . . , n i , and i = 1, . . . , m, be a binary indicator such that Y ij = 1 denotes the event that the ith individual has a positive urinalysis test during the jth clinic visit with urinalysis and Y ij = 0 otherwise. To relate the observed test data to the available covariates, we posit the following generalized linear functional mixed-effects model where g(·) is the logit link function which is used to relate the linear predictor, ν ij , to the probability of relapse, π ij = P(Y ij = 1). To elucidate the key feature of our model adopted to capture the effect of dose adherence, we note that the first term on the right-hand side of (1) is the functional component which consists of the time-varying buprenorphine dose curve (D ij (·), e.g., see Figure 1), the functional coefficient (β(·)), and the time frame (A ij ) leading up to the urinalysis clinic visit over which the dose levels are allowed to impact the probability of relapse. The remaining components of the model consist of x ij a P-dimensional vector of demographic and drug use history risk factors whose first entry is a one to allow for the usual intercept, α the corresponding vector of regression coefficients, γ 0i a subject-specific random effect entered into the model to account for the heterogeneity across subjects, γ 1k(i) = γ 1k if the ith subject is part of the kth trial, and γ 1k is a random effect specified to account for the heterogeneity across trials, k = 1, . . . , K. Herein, the random-effects distributions were taken to be and note, here the random effects are taken to be independent given the nesting of subjects within trials. A few comments on the form of the model are warranted. First, through adopting the functional regression framework, we are able to directly acknowledge and estimate the effect of time-varying dose adherence, whereas more traditional variable aggregation techniques (e.g., average dose) fail to acknowledge key trends in dose adherence, e.g., waning adherence from the point of care or weekly patterning. Second, the time window (i.e., A ij ) should be selected so that the upper bound is just before the jth clinic visit with urinalysis for the ith individual and that the length of the interval reflects the approximate elimination time for buprenorphine, i.e., buprenorphine doses taken prior to the lower bound are no longer present in the patient's system and therefore cannot impact the probability of opioid use. Generally speaking, it typically takes five half-lives for a drug to completely leave a subject's system. Thus, given that the elimination half-life of buprenorphine is 24 to 42 h, we specified a time window consisting of 15 days to more than adequately capture the relevant dose history. Lastly, given the form of the proposed model, we can extract dose effect for individuals with strict adherence to their prescribed dosing regime. To see this, we note that if a subject adheres to the dosing regime, then D ij (s) = D ij for all s. Thus, we would have Note, β * represent the usual increase in log odds associated with a one unit increase in buprenorphine dose level. Thus, by estimating β(·), we can also estimate β * .
Estimating the buprenorphine dose effect β(·) in model (1) is challenging from both a theoretical and computational perspective because of its infinite dimension. To reduce the number of unknown parameters needed to be estimated while also maintaining adequate modeling flexibility, we approximate β(·) using B-splines [7]. This leads to the following representation of β(·): where b (·) is a spline basis function and η is the corresponding spline coefficient, for = 1, . . . , L. The L spline basis functions are fully determined once the degree and knot set are specified, thus the only unknown parameters in (3) are the spline coefficients. In specifying the basis functions, the degree controls the overall smoothness of the basis functions and the number of knots determines the overall modeling flexibility; for further discussion see [7]. We suggest selecting a relatively large knot set (e.g., 5-6 knots) and then regularizing the estimation of the spline coefficients through the methodology outlined below.
Using the spline representation of β(·) depicted in (3), we can re-express the functional component in model (1) as follows where m ij is an L-dimensional vector whose th element is m ij = A ij D ij (s)b (s)ds and η = (η 1 , . . . , η L ) . Thus, the linear predictor of our model can be expressed as

Prior Specification
To facilitate both parameter estimation and inference, we cast our problem into the Bayesian paradigm. The first step in this process involves specifying prior distributions for all unknown parameters. Given the complexities of our problem, priors are chosen to regularize the estimation of the parameters. In particular, the prior for the spline coefficients is designed to encourage smoothness in the functional estimate while the prior for the regression coefficients is meant to "shrink" unimportant variables toward zero. In what follows, we briefly expand on these specifications.
To avoid overfitting issues and to encourage smooth functional estimates, herein we adopt a prior for the spline coefficients which leverages a covariance structure inspired by the usual roughness penalty [8]. This common penalty encourages smoothness by penalizing for abrupt changes in the function through the following: where β (2) (·) denotes the second derivative of β(·) and R is an L × L matrix with entries (2) (·) being the second derivative of b (·). Note, the spline representation adopted for β(·) is key to being able to represent this penalty as the quadratic form depicted above; for details of this derivation see [8]. Capitalizing on the structure of this penalty and the duality that exists between regularized estimation and prior distributions in the Bayesian paradigm, we specify the following smoothing penalty inspired prior distribution for η: In the prior specification above, the additional variance parameter λ governs the amount of smoothness and controls the trade off between over and underfitting the data.
To aid in variable selection, we adopt the generalized double Pareto shrinkage prior, proposed by [9], for all of the regression coefficients with the exception of the intercept, i.e., we specify where GDP(ψ, a δ ) refers to the generalized double Pareto distribution [9]. Under these prior choices, setting τ 0 to be large provides a vague prior on α 0 , while the hyperparameters a δ > 0 and b δ > 0 govern the amount of shrinkage. In particular, these parameters control the dispersion, with a δ controlling the heaviness of the tails of the distribution. A typical default specification, and the one adopted herein, is to set a δ = b δ = 1 which leads to Cauchy-like tail behavior which is known to have desirable Bayesian robustness properties [9].
Finally, we place inverse gamma priors on the variance components of the random effects, i.e., we specify This specification is common among the literature and it leads to a proper posterior [10]. Based on the prior specifications outlined above, we develop a Markov chain Monte Carlo (MCMC) sampling algorithm which facilitates both posterior estimation and inference.
In what follows, we provide a brief overview of this algorithm and its construction.

Data Augmentation and Posterior Sampling
With ease of implementation and computational efficiency in mind, herein we outline the construction of a posterior sampling algorithm that consists solely of Gibbs steps [11]. To accomplish this, we consider a two-stage data augmentation process. The first stage follows the work of [12], and introduces carefully constructed Pólya-Gamma latent random variables so that the logistic function can be hierarchically expressed as a scale mixture of normals, where the mixing distribution is Pólya-Gamma; for further details, see [12]. The second stage decomposes the generalized double Pareto shrinkage prior as a scale mixture of normals; for further discussion see [9]. For the specific details of this two-stage data augmentation process, see Appendix A.
The data augmentation scheme outlined above leads to the following full conditionals α|Y, w, η, γ 0 , γ 1 where the specific form of the parameters of these distributions are given in Appendix A. These full conditionals were used to construct an MCMC algorithm in the usual manner; for further discussion see [11].

Trial Data
Clinical trial data for this analysis were sourced from the Clinical Trials Network (CTN) at NIDA's Data Share resource (https://datashare.nida.nih.gov/ (accessed on 5 September 2017)). Using the search keyword opiate, we identified 10 efficacy and safety trials involving detoxification or maintenance treatment of DSM-IV opioid dependence. We selected six efficacy and safety trials focused on buprenorphine maintenance treatment for analysis. Detailed information on these trials are provided in Table A1.

Patient Characteristics
The data consist of 55,739 urinalysis results from 3022 subjects who participated in one of the six aforementioned clinical trials aimed at assessing the efficacy of buprenorphine for treating OUD. The number of urinalyses (i.e., urine drug screens for opioids) per subject ranged from 1 to 60 urinalyses, while the mean number of urinalyses per subject was 18.44 and the median was 18. The data were harmonized across the six trials and candidate predictors with a missingness greater than 25% were filtered out. This resulted in 18 demographic, sociodemographic, and substance use variables (excluding prescribed buprenorphine dose, handled by the functional component of the model). Missing demographic, sociodemographic, and substance use variables were imputed using the regularized iterative factorial analysis for mixed data (qualitative and quantitative variables) algorithm [13], implemented by the imputeFAMD function in the missMDA R package. Summaries of the retained variables (with imputed values included) are given in Table 1. The daily dose of buprenorphine taken by each patient was either reported (CSP-999) or inferred from alternate information. Daily dose could vary throughout time for a variety of reasons, e.g., adherence, induction, re-induction after lapse in dosing, and modification by a provider's clinical judgement.
Given the number of demographic and substance use variables considered, the reference group is specifically white men with a high school diploma who are employed full time doing skilled manual labor, married and living with a partner or child, and their primary mode of opioid use being intravenous, with a history of heroin, cocaine, alcohol and marijuana use and no history of methamphetamine, tranquilizer, or PCP use. The mean age, income, and years of opioid use are 36.14 years, $20,834 per year and 8.23 years, respectively (presented in Table 1), while the mean dose is 12.65 mg/day (presented in Table 2). When we discuss conditional probabilities of relapse, comparisons will be made with respect to this hypothetical individual in the reference group by changing specific variables as noted.

Functional General Linear Mixed Model
The outcome variable in this analysis was the urinalysis test result for illicit opioid use (1 = positive drug screen vs. 0 = negative drug screen). Through the model in (1), we relate the daily dose patterns leading up to the clinic visit with urinalysis, while controlling for the 18 demographic, sociodemographic, and substance use variables detailed in Table 1. For the functional dose component in model (1), the time trajectory was chosen to be the 15 days leading up to the current urinalysis clinic visit and, for the B-spline basis expansion of the coefficient function in (3), we specify the degree to be 3 to construct cubic basis functions. Two interior knots were placed at the 33.33th and 66.67th percentiles of our 15 day time range. This leads to five fully determined spline basis functions and, hence, five spline basis coefficients to estimate. For the priors outlined in Section 2.2, we take τ 0 = 1000 to specify a vague prior on the global intercept α 0 and let a 0 These hyperparameter values are chosen so to produce uninformative, proper prior specifications. For sampling, we retain 5000 MCMC iterates after a burn-in of 5000 samples. Convergence of the MCMC chains was assessed in the usual manner, i.e., trace plots. To summarize our analysis, we report the estimated posterior means (point estimates of the effects), estimated posterior standard deviations (measures of uncertainty), and 95% equal-tailed credible intervals. Figure 2 summarizes the estimated functional coefficientβ(t) (black solid line), which represents the buprenorphine daily dose effect for the 15 days leading up to a clinic visit with urinalysis. The dashed lines are the 95% credible interval limits. On the horizontal axis, if t is the day of the current urinalysis clinic visit, then t − 15 represents 15 days prior and t − 1 represents one day prior. Table 3 reports the demographic and substance use variables that were found to be significant. Of the 54 fixed effects, four were deemed to be important by the model (i.e., their estimated credible intervals did not contain zero). Table 3 summarizes these significant factors by reporting the estimated posterior means (point estimate of the effect), estimated posterior standard deviations (measure of uncertainty), and 95% equal-tailed credible intervals. The analogous results for the full set of demographic and substance use variables are provided in Tables A2 and A3 in Appendix C.2.  As previously stated, daily dose adherence was only directly recorded for patients in the CSP-999 trial. Specifically, while the assigned daily dose of buprenorphine was recorded as a part of the five other trials, adherence was not. For these trials, dose adherence was reconstructed using other available information, e.g., self reported non-adherence and returned pills. To examine how the buprenorphine dose reconstruction could have impacted our results, we reran our analysis on data from the CSP-999 trial only. A summary of these results can be found in Appendix C.
To extract an estimate of dose effect for subjects that were strictly adherent to their assigned dosing regime, we compute the following integral for each realization β(s), denoted β (g) (s), drawn from the posterior with β * (g) being a posterior realization of β * . Table 4 provides a summary of these results for both the full and reduced (CSP-999) analysis to include the posterior mean estimate (point estimate of the effect), estimated standard deviation of the posterior (measure of uncertainty), and 95% equal-tailed credible interval. Table 4. Analysis results: Summary includes the posterior mean estimate (Est), the estimated standard deviation of the posterior (ESE), and the estimated 95% equal-tailed credible interval (CI95) for the dose effect (i.e., β * ) for the full and reduced (CSP-999) analysis.

Discussion
The primary focus of our analysis is two-fold. First, we wish to demonstrate a novel approach to account for non-adherence that commonly arises in medication-assisted treatment trials; especially those targeting substance use disorders. Second we wish to refine the understanding of the effectiveness of buprenorphine as a MOUD, while accounting for the potential non-adherence of study patients. To accomplish both of these tasks, we investigated the influence of multiple demographic, sociodemographic, drug use history, and treatment variables on the risk of illicit opioid use with publicly available individual patient data from six Federally sponsored buprenorphine efficacy and safety trials. To acknowledge and account for patterns of non-adherence, we conceptualized the daily dose histories of the study patients as a functional covariate and we estimated an associated functional effect. A summary of this estimated functional is provided in Figure 2. From these results, we identify several key findings. First, these results suggest that buprenorphine, as an MOUD, significantly reduces the risk of illicit opioid use. This can be seen from the fact that the point estimates, and associated credible intervals, are all below zero, i.e., the integral over the the product of this functional and D ij (·) ≥ 0 results in a negative quantity. Second, we find that dose history extending to approximately 12.5 days prior to an opioid screening visit is significantly related to the risk of short-term lapses. Third, the risk of illicit opioid use is related to dosing adherence patterns throughout the considered 15 day window leading up to the urinalysis, although, recent patterns have more influence. This can be seen from the decreasing nature of the functional estimate, especially for the five (approximately) days before the urinalysis test. Fourth, based on our estimated functional, we are able to extract an estimate of dose effect for subjects that were strictly adherent to their assigned dosing protocol. Based on this approach, we estimate that the log-odds of short-term lapse decreases by 0.09 with every 1 mg/day increase in dose; see Table 4. This new assessment of the efficacy of buprenorphine as an OUD treatment is unobscured by the effects of non-adherence and leverages six NIDA-sponsored efficacy and safety trials to render its conclusions.
When examining the association between risk of illicit opioid use and the other demographic, sociodemographic, and drug use history variables, four of the 54 were found to be significant. In particular, we find that increasing age is protective, while being unemployed, a drug use history of using heroine, and a drug use history of using opioids intravenously are associated with an increased risk of illicit opioid use. A similar finding that increasing age is associated with no positive urine drug screen was recently reported in an analysis of Veterans Administration patients undergoing buprenorphine treatment [14]. The protective nature of employment for patients in recovery [15] is concordant with unemployment being identified as a risk factor for illicit opioid use. Older age, no heroine use history and no IV drug use have already been reported as protective with respect to successful opioid use outcomes (abstinent during week 24 and ≥2 of the previous 3 weeks) in a secondary data analysis of the Prescription Opioid Addiction Treatment Study (POATS or CTN-0300) [16,17], one of six CTN trials included in this study. Notably, the largest protective effect we observed was "primary mode of opioid abuse" with the log-odds of short-term lapse decreasing by 1.32 when the primary mode of abuse is oral. This could be attributable to the severity of the opioid use disorder, with intravenous use being a hallmark of more severe cases.
When examining the results of the sensitivity analysis (see Appendix C) of the CSP-999 trial only, we note several similarities and differences. Importantly, the full and CSP-999 analysis came to virtually the same conclusions with regard to the efficacy of buprenorphine as an OUD treatment. In particular, the estimates of β(·) are not statistically different from each other. However, the estimate from the full analysis is slightly attenuated toward zero when compared to the CSP-999 only analysis. This feature can also be observed in the effect estimate reported in Table 4. A plausible explanation for this would be that our approach to reconstructing dose histories for the study patients, though effective, was not perfect, and therefore introduced "measurement error" into this variable. A hallmark of measurement error is the attenuation of effect estimates toward zero, e.g., see [18]. When comparing the estimated effects of demographic, sociodemographic, and substance use variables we find that most are not statistically different, yet there are differences in those deemed to be significant by the two analyses. These differences are likely attributable to increased precision due to larger sample sizes in the full analysis and differences in the demographic distribution across the full and reduced data.
We also acknowledge the lack of other risk factors that could be used to better understand/predict short-term lapse. Inclusion of time-varying predictors such as current stress levels, occurrences of major life events (e.g., familial death and loss of job), and other psychological measures would undoubtedly enhance our model. However, the impact of not having these variables is mitigated by the adoption of subject specific random effects.
Importantly, this study was specifically aimed at estimating the effectiveness of recent buprenorphine treatment at reducing short-term lapse. With that being said, this analysis did not consider adherence and its impact on illicit opioid use over longer periods of time and the subsequent associations with OUD-related adverse outcomes, which is a far more complex problem. Studies aimed at these more long-term outcomes could reveal OUD treatment strategies that would be poised to positively impact public health. That is, there have been nearly 500,000 overdose deaths from opioids in the United States alone in the last 20 years [1]. Further, OUD-related mortality appears to be increasing. Specifically, the CDC estimates that overdose deaths from opioids increased to 75,673 in the 12-month period ending in April 2021, up from 56,064 in 2020 [19]. Moreover, less than one-third of patients enrolled in comprehensive health care with current OUD are being treated with one of three approved medications for OUD [20]. Extended MOUD treatment (>1 vs. ≤1 year) appears to reduce mortality [21]. Thus, more in-depth studies relating MOUD treatment to long-term outcomes has the potential to identify OUD treatment strategies that can be more effectively utilized to treat this epidemic and shift current clinical practice. Future research efforts will be aimed at studying these more complex topics related to dose, adherence and treatment outcomes and their association with OUD-related mortality rates.

Conclusions
Inspired by adherence issues in six efficacy and safety clinical trials of buprenorphine maintenance treatment, this work proposed a generalized linear functional mixed-effects model that can acknowledge and account for the effects of adherence when assessing treatment effect. The proposed model was applied to the six motivating clinical trials in an effort to better refine understanding about the time-dependent effect that buprenorphine has on treating OUD. In particular, we find that dose history approximately 12.5 days prior to an opioid screening visit is significantly related to the risk of short-term lapses, with the more recent history being more impactful. Further, we are able to extract an estimate of dose effect that is not obscured by adherence issues. That is, we estimate that the log-odds of short-term lapse decreases by 0.09 with every 1 mg/day increase in buprenorphine dose.  Using the search keyword opiate, we identified 10 efficacy and safety trials involving detoxification or maintenance treatment of DSM-IV opioid dependence. We selected six efficacy and safety trials focused on buprenorphine maintenance treatment for analysis (Table A1).
Conflicts of Interest: J.W.B. is a LLC member and employee of BioRealm LLC. BioRealm LLC offers data science services. All other authors declare no conflict of interest.

Appendix A. Posterior Distributions
We assume conditional independence given the covariate effects and random effects and observe that Y ij depends on the model parameters only through the linear predictor, ν ij . Hence, the likelihood can be expressed as where g(·) is defined to be the logit link function.
We develop a two-stage data augmentation process to construct a posterior sampling algorithm consisting only of Gibbs steps. In the first stage, we exploit a hierarchical representation of the proposed data model by introducing Póyla-Gamma latent random variables w ij ; for further details, see [12]. Under this specification, the joint density of the observed and latent data for the ith individual is given by }, and f (w ij |a, b) denotes the Pólya-Gamma density with parameters (a, b); for further details, see [12].
Attention is now turned to the second stage of the data augmentation process and the construction of the hierarchical representation of the joint posterior distribution. Recall from Section 2.2, we specify a generalized double Pareto shrinkage prior for all of the regression coefficients with the exception of the intercept, i.e., α 0 ∼ N(0, τ 0 ), α p ∼ GDP(ψ = b δ /a δ , a δ ), for p = 1, . . . , P − 1.
As noted by Proposition 1 of [9], the generalized double Pareto shrinkage prior can be represented as a scale mixture of normal distributions. Thus, for the regression coefficients, the following hierarchical representation provides the same prior specifications as those given above: The δ p s are the global shrinkage parameters, while the τ p s are the local shrinkage parameters and override the impact of the global shrinkage components for the variable fixed effects that are not near zero [9].
In deriving the full conditional distributions, for notational convenience, a dot · is used as shorthand for all the parameters one is conditioning on, e.g., we may write the posterior p (α Y, w, η, γ, τ) as p(α ·).
We begin by deriving the full conditional distribution for the spline coefficients η, based on the smoothing penalty inspired prior distribution outlined in Section 2.2. Letting Recognizing this as the kernel of a normal density, we have η · ∼ N(µ η , Σ η ).
Further, the full conditional distribution for the variance parameter λ associated with the smoothing prior for the spline coefficients η is derived as follows.
Moreover, the full conditional distribution for the global shrinkage parameters δ p is derived as follows. Exploiting the fact that the Laplace density is a scale mixture of normals with an exponential mixing density [22]: where a * δ p = a δ + 1 and b * δ p = b δ + |α p |. Recognizing this as the kernel of a Gamma density, we find that δ p · ∼ Gamma(a * δ p , b * δ p ), for p = 1, . . . , P − 1. We turn our attention to the random effects and derive the full conditional for the subject-specific random effects, γ 0 = (γ 01 , . . . , γ 0N ) T , where N is the number of participants in this study.
Recognizing this as the kernel of a normal density, we find that γ 0 · ∼ N(µ 0 , Σ 0 ).
The full conditional distribution for the variance component of the subject-specific random effects, σ 2 0 , is derived as follows.
The full conditional distribution for the variance component of the trial-specific random effects, σ 2 1 , is derived as follows.
where a * 1 = a 1 + K/2 and b * 1 = b 1 + 0.5 · γ T 1 γ 1 . Recognizing this as the kernel of a Gamma density, we find that σ 2 1 ∼ Inv-Gamma(a * 1 , b * 1 ). These full conditionals were used to construct an MCMC algorithm in the usual manner. To validate the proposed model and MCMC algorithm, an in-depth numerical study was conducted. This study was designed to emulate the primary features of the opioid data under study. The results (not shown) of this numerical study suggest that our proposed approach performs well and is appropriate for analyzing the motivating data.  The data consist of 15,983 urinalysis results from 654 subjects who participated in the CSP-999 trial. The number of urinalyses per subject ranged from 1 to 59, while the mean number of urinalyses per subject was 24.44 and the median was 17. Summaries of the demographic, sociodemographic, and substance use variables are given in Table A5. Several of the variables summarized in Table 1 are no longer available when only examining CSP-999 trial In particular, CSP-999 trial patients are either white, Hispanic, black, American Indian, or Asian and hence, the race categorized as "Other" was removed. None of the patients work for the military, so the work type category "Military" was removed. All of the patients have a stable living arrangement, so the "No Stable" living arrangement category was removed. All CSP-999 trial patients use heroine, so the "Heroine Use" categorical variable was removed. Finally, none of the patients' chosen mode of opioid abuse was sublingual, so the "Sublingual" category for the mode of opiate abuse variable was removed. For ease of comparison, the reference group is the same as the one used for the full analysis. The mean age, income, and years of opioid use are 36.16 years, $19,454 per year and 11.56 years, respectively (presented in Table A5), while the mean dose is 7.49 mg/day (presented in Table A4).

Appendix C.2. Functional Generalized Linear Mixed Model
Through the model in (1), we relate the daily dose patterns leading up to the clinic visit with urinalysis, while controlling for the 17 demographic, sociodemographic, and substance use variables detailed in Table A5. Note that, because all of our data comes from the CSP-999 trial, the trial-specific random effects γ 1k(i) are removed from the model. Figure A1 summarizes the estimated coefficient functionβ(t) (solid black line), which represents the buprenorphine daily dose effect for the 15 days leading up to a clinic visit with urinalysis. For comparison purposes, the coefficient function estimated from the full (all trials) analysis is also plotted (solid red line). The 95% credible intervals estimated from the full and reduced analyses are also plotted (red and black dashed line, respectively). Table A6 reports the demographic and substance use variables that were found to be significant. Of the 49 variable fixed effects, 5 were deemed to be important by the model (i.e., their estimated credible intervals did not contain zero). Table A6 summarizes these significant factors by reporting the estimated posterior mean (point estimate of the effect), estimated standard deviation of the posterior (measure of uncertainty), and 95% equal-tailed credible interval for each parameter. The analogous results for the full set of demographic and substance use variables are provided in Tables A2 and A3.   Table A5. Sociodemographic characteristics and drug use history for the individuals used in the CSP-999 trial analysis.