Item response theory (multivariate qualitative to quantitative latent variable transformation)

In the application of the IRT method in handling multivariate phenotypic data, it was assumed that each individual classified with RA had a latent RA variable. Such a latent feature for RA was the result of a statistical summary of information from nine observed risk factors. These factors were qualitative anti-CCP, qualitative IgM, gender, smoking status, and five severity classes of left and right hand erosions for RA, which were transformed to individual severity dummy variables as 1/0 for each corresponding combined class (see Data Input). Applying IRT, each member of the sample was assigned a latent (z _i ) score. Our analysis is constraint to one latent variable z_1i. The probability density function of y _i given the latent variable z can be written as $p(y_i|z,x,{\theta }_i,{\varnothing }_i)=e\left\{\frac{y_i{\theta }_i-{\beta }_i({\theta }_i)}{{\alpha }_i({\varnothing }_i)}+d_i(y_i,{\varnothing }_i)\right\}$, where θ _i and ∅ _i are the location and dispersion parameters, respectively. In the IRT, the item characteristic curves represent the positive response probability of each risk factor j relative to the latent RA score. Formally, this relationship is expressed as $p\left(y_{ij}=1|z_i\right)=\frac{1}{1+e\left\{-{\alpha }_i(z_i-{\beta }_i)\right\}}$, and graphically it represents the importance of each risk factor toward the latent RA. The discrimination parameter α _i is viewed as how well a risk factor discriminates among subjects with opposite extreme latent RA. The difficulty parameter β _i characterizes the difficulty level of individual risk variables. In this study we applied the two-parameter IRT model (α _i , β _i ), which was constrained to one latent variable (z _i ). Software for IRT analysis in R language is developed by Rizopoulos (ltm package version 0.5–1 for Linux OS and R version 2.3.0) [5, 6].

Blom score (quantitative univariate to quantitative rank score transformation)

In the second method we used the Blom transformation. Blom scores represent rank approximations of the exact order of a normal distribution. In the group of rank transformed quantitative variables, one extracts a Blom score by applying the following formula on the anti-CCP.

$Blo{m}_i=\frac{{\varphi }_i^{-1}(r_i-\frac{3}{8})}{n+\frac{1}{4}}$, where ϕ _i ^-1 is the inverse of cumulative normal function, r _i is the rank of observation i, and n is the number of non-missing observations. The Blom scores of anti-CCP represented a better normal distribution than the original values. Software for the Blom transformation is available within the PROC RANK of SAS, v 9.1.3 for Linux OS.

Linear mixed effects models (association tests)

The LME model used in our study to test the association between an SNP and changes in the response variable (IRT RA latent variables/anti-CCP Blom transformed variables) follows in a matrix form: Y= XB+ ZU+ ε, where Yis an m × 1 vector of responses; Xis an m × p design matrix of the fixed effects; Bis the parameter p × 1 vector of fixed effects; Zis an m × q incidence matrix of random effects, and Uis a q × 1 vector of random effects with E(U) = 0, and covariance matrix G; 0is an m × 1 vector of random effects with E(0) = 0 and covariance matrix R. In the fixed effects we included SNP genotypes recoded as additive effects (-1 for one homozygote, 0 for the heterozygote, and 1 for the other homozygote genotype), and gender. In the random effects we included the family identification number. We tested whether the SNPs additive effects are different from zero, and especially we identified the highest significances, considering that multiple comparisons as well as correlated SNPs tests because of linkage disequilibrium are present. The LME in SAS (v. 9.1.3 for Linux OS) was applied with PROC MIXED using the option for EMPIRICAL variance, which computes the estimated variance-covariance matrix of the fixed-effects parameters by using the following asymptotically consistent estimator ${\left(X{{}^{‵}V}^{-1}X\right)}^{-}\left({\displaystyle {\sum }_{i=1}^s{X}_i{{}^{‵}V}_i^{-1}{\in }_i{\in }_i`{{}^{‵}V}_i^{-1}{X}_i}\right)\left(X{{}^{‵}V}^{-1}X\right)`$, where var(Y) = V= ZGZ' + R[7].

Results of association tests in 50 replications of the simulated data were analyzed with two definitions: sensitivity, defined as the proportion of the true causative SNPs that had a positive association result (p-value < 0.05); and specificity, defined as the proportion of true non-causative SNPs that had non-significant associations (p-value ≥ 0.05).

Data input

NARAC phenotypes in our analysis included anti-CCP, rheumatoid factor (RF), gender, severity of left and right hand side erosions on a scale of 1 to 5 (severityLH and severityRH, respectively), whether or not the subject smoked cigarettes (HxCigN). Qualitative dummy variables were created for each of the following variables, where the subject with anti-CCP > 49 [8], RF > 170 (to capture approximately the RF upper quartile in a very skewed RF distribution); smokers (HxCigN = 2); gender (if female); and whether severityLH or severityRH was 1, 2, 3, 4, or 5 were considered as affected. In the corresponding nine dummy variables (accp, igm, smok, sex, sev1, sev2, sev3, sev4, and sev5) affection was codes as 1, unaffected as 0, otherwise as missing. The same nine dummy variables were created also for the RA GAW15 simulated data. These dummy variables were analyzed individually for the two sources of data with the ltm package to extract the corresponding RA latent factors. Also, the anti-CCP variables were independently Blom-transformed within groups defined by gender and each age decade. LME model was performed to evaluate association of these traits, with 404 Illumina SNPs on chromosome 6 for NARAC data, and with 17,820 dense simulated SNPs on chromosome 6 for 50 replications on the simulated GAW15 data. In the NARAC data the sample size was 1499 subjects and 757 families, and 1340 subjects and 757 families for anti-CCP Blom transformation and IRT latent variable, respectively. In the simulated data the sample contained a mean of 3497 ± 20 subjects and 1500 nuclear families for both traits.