Statistical analysis

The subcohort (n = 236) was generated with a random sampling fraction of 0.45 and included 54 cases. Cases that were not in the subcohort (n = 59) entered via delayed entry. The time axis was follow-up time in study. Cox proportional hazards was used for analysis, with robust variance and Barlow weights to account for the case cohort design [36].

Descriptive analysis was conducted in the subcohort using box plots, histograms, correlation matrices with Spearman’s correlation coefficient, scatter plots and Q-Q plots. Values of the markers that were near the limits of detection were assigned the lowest detectable value.

Predictors were selected based on their availability to physicians, strength of correlation with each other [37], magnitude of associations in univariable analyses, ability to improve model fit as indicated by the Akaike Information Criterion (AIC) and ability to maximize discrimination as measured by Harrell’s C [38]. For immune markers, variables were also included in models if they captured a different stage of the underlying etiology of fibrosis development or were linked to fibrogenesis in the literature.

Univariable and multivariable analyses were conducted on untransformed variables, as well as after log-transformation or using median or quartile distributions of the immune markers. Continuous variables were centered at their mean values. Log transformation was used for all included continuous variables such as baseline APRI and immune markers IL-8, sICAM-1, RANTES, hsCRP, and sCD14 but not for age, which was modeled as a restricted cubic spline with 3 knots at the 10th, 50th and 90^th percentiles corresponding to ages 33, 44, and 54. HIV viral load was dichotomized (undetectable or not at ≤ 50 copies/ml), as were alcohol use (currently drinking or not); HCV genotype (3 vs. non-3, i.e. types 1, 2, and 4); and host IFNL genotype rs8099917 (TT vs. non-TT).

Proportional hazards were assessed using the using Stata command–stphtest, detail- which uses scaled Schoenfeld residuals to check if proportional hazards holds globally and for included predictors.

Discrimination, calibration and changes in reclassification were compared between Model 1 (selected clinical predictors only) and Model 2 (clinical predictors from Model 1 plus selected genetic and immune markers) for predicting 3-year risk of significant liver fibrosis.

Discrimination was measured with a weighted [39] Harrell’s C or concordance index using Stata command–somersd- with robust jackknife estimator for standard errors [38]. A C-index value of 1 indicates perfect discrimination, while 0.5 means no better than random guessing. Calibration was assessed statistically (Hosmer-Lemeshow statistic and the Gronnesby and Borgan (GB) test) [40], and graphically with the Stata command–stcoxgrp using imputed data [41].

Change in reclassification was measured by the net reclassification improvement (NRI) summary index [42]. We calculated both the category-based and the continuous NRIs. For the category-based NRI, we used 3 clinically relevant risk categories: low risk, < = 10%; medium risk, >10–25%; and high risk, >25%. The categories were determined based on estimates of mortality from liver disease in those with chronic HCV infection from published reports [43] as well as opinions of knowledgeable hepatologists and clinicians. For the continuous NRI, no categories were needed and any upward or downward movement in risk was considered, regardless of magnitude.

The models were internally validated using bootstrapping. All analyses were conducted using Stata 12.

Missingness for plasma samples and other variables was assumed to be at random. Multiple Imputation by Chained Equations (MICE) was used on the full cohort to account for all missing data, using all the predictors in the final models, all the immune markers, as well as variables that were possibly related to the reasons for missingness.

We compared predictive accuracy after using unweighted Cox proportional hazards regression on the imputed full cohort data [44].