Estimating Patient-Specific Relative Benefit of Adding Biologics to Conventional Rheumatoid Arthritis Treatment

Key Points Question Should biologics be added to conventional treatment for specific patients with rheumatoid arthritis considering the short-term benefit? Findings In this meta-analysis of individual participant data from 3790 patients, the addition of certolizumab, a tumor necrosis factor α inhibitor, to conventional rheumatoid arthritis treatment was associated with an increased probability of reaching low disease activity in general. However, an interactive application based on a 2-stage model, which can visualize the estimated results for individual patients, showed that the absolute benefits differed in patients with different baseline characteristics. Meaning These findings suggest that the benefit of adding biologics to conventional rheumatoid arthritis treatment is uncertain for specific patients, and the interactive application based on a 2-stage model may help with treatment selection in practice.

-Different participants (e.g., patients who have been using CTZ, or who have achieved low disease activity before randomization, etc.) (n = 11) -Different interventions (e.g., CTZ used as a monotherapy, or together with other biologics) (n = 5) -Different controls (e.g., the control group also received CTZ or other biologics) (n = 9)   Calibration slope 1.0101 0.9699 1.5037 0.9751 1.1222 Abbreviations: BMI: body mass index; csDMARD: conventional synthestic disease modifying anti-rheumatism drug; bDMARD: biologic disease modifying anti-rheumatism drug; MTX: methotrexate; NSAID: nonsteroidal anti-inflammatory drug; CRP: C-reactive protein; ESR: erythrocyte sedimentation rate; RF: rheumatoid factor; HAQ-DI: Health Assessment Questionnaire-Disability Index; AUC: area under the curve. a Parameters were combined from 10 imputed datasets according to the Rubin's rule. b LASSO models were estimated using two selected tuning parameter : _MAX to maximize the AUC, _1SE by the one-standard-error rule. c These variables were transformed to resolve skewness before model development.
Then, is the log odds of the outcome for the control group in trial , assumed independent across trials; is the treatment effect in terms of log OR in trial , assumed exchangeable across trials; hence, is the summary estimate of the log-odds ratios for the intervention versus the control arm, and 2 is the heterogeneity of across trials. We used R2Jags package to fit this model 1 .

Model: stage one
We fitted three penalized logistic regression models to estimate the probability of remission or low disease activity using baseline characteristics regardless of treatment (i.e., treatment was not included as a covariate in the model). Penalization methods to shrink the coefficients were used to avoid extreme estimations 2 . Before developing the model, we transformed the continuous variables which displayed noticeable skewness and checked the correlations between the continuous covariates (results are shown in eAppendix section 4. The models were: (1) LR22+PML: Logistic regression of 22 covariates + penalized maximum likelihood (PML) shrinkage methods, where the tuning parameter was selected as the one that can maximize a modified Akaike's information criterion (AIC) 2 3 ; the lrm and pentrace function from the rms package were used 4 .
(2) LR22+LASSO: Logistic regression of 22 covariates + LASSO (least absolute shrinkage and selection operator) method, where two tuning parameter λ values were selected based on 10-fold cross validation AUC: one to maximize the AUC, and another by the one-standard-error rule (i.e., the most parsimonious model whose AUC is no more than one standard error lower than that of the best model) 5 6 ; the glmnet command from glmnet package was used 7 . The logistic regression part of all the models had the same form: ~( ) is the probability of the outcome for patient from trial at baseline; is the th prognostic factor in study for patient . Then, 0 and are the intercept and the regression coefficient for the th prognostic factor in study ; they are fixed across studies.
The outcome variable was measured repeatedly during the follow-up. Patients who had used rescue therapy before a follow-up visit were defined as non-responders at that visit, despite the measured disease activity scores. Then, we checked the missing pattern over time in the outcome variable for each study. We considered patients who had no outcome data at consecutive visits starting before three months until six months were unlikely to achieve response, thus we defined them as nonresponders at the 3-month visit. Other than this situation, patients who had missing outcome data at the 3-month visit were treated as missing values. Finally, we used multiple imputation chained equation (mice function from the mice package 9 ) to handle missing values in the covariates and outcomes for each model based on the missing-at-random assumption 10 11 . The variables used in the imputation model were the same as in the substantive model. Predictive mean matching was used for continuous variables, while logistic regression imputation for categorical variables. We imputed 10 datasets for each model and checked the convergence. We fitted the above three models on each imputed dataset, and combined the coefficients based on Rubin's rule.
We conducted bootstrap internal validation to estimate the optimism-corrected AUC, calibration intercept and slope for three stage-one models 2  Then, is the log odds a patient whose baseline expected probability of the outcome equals to the mean probability in the control group, and is assumed to be independent across trials; To solve it, we imputed the baseline expected probabilities by incorporating a Bayesian model, which took the reported aggregate-level means and standard deviations or proportions of these variables, to assume the prior distributions.
In order to estimate a new patient who is not from any trials and has a logit baseline expected probability of the outcome ( ) estimated from the stage-one model, we used the following equations: , 0 and were estimated in stage-two, ( ) was estimated as the mean of ( ) over all the individuals from all included studies; and was estimated as the mean log-odds in all the control arms. Then we estimated the individual probability of the outcome if receiving CTZ+csDMARDs or placebo+csDMARDs respectively, and estimated the absolute risk difference between two groups. We used R2Jags package to fit this model 1 .
Note that this stage-two model could not be validated for its risk difference estimation like the stage-one model. It is because stage-two model aimed to estimate the risk difference of adding CTZ vs no CTZ. In order to evaluate whether the estimation is correct or not, we need the actual outcome of each patient. In this case, both the outcome when the patient received CTZ and the outcome when he/she did not receive CTZ were necessary. However, as a specific patient was treated by either CTZ or no CTZ, we could not acquire necessary data to assess the model performance. New methods have been proposed to evaluate the performance of models aiming to estimate personalized treatment effect between two groups very recently 13 , which could be considered in future research. (2) Age (years): *Note: Only age band was provided for each patient in the dataset, and the cut-offs were different from study to study. Therefore, we took the median age in the band for each patient.   There are two pairs whose correlation coefficients are higher than 0.6: patient's global assessment and patient's assessment of pain: 0.85; tender joint count and swollen joint count: 0.62. Considering these two variables were both critical items in most disease activity indices, and models using the PML method did not perform variable selection, we kept both variables in all the models.