RealROC: a shiny based application for ROC curve study with covariate adjustment

Universidade do Minho Escola de Engenharia Departamento de Inform´ atica Francisco Lu´ ıs do Amaral Ribeiro Machado e Costa RealROC A Shiny based application for ROC curve study with covariate adjustment October 2020 Universidade do Minho Escola de Engenharia Departamento de Inform´ atica Francisco Lu´ ıs do Amaral Ribeiro Machado e Costa RealROC A Shiny based application for ROC curve study with covariate adjustment Master dissertation Master Degree in Bioinformatics Dissertation supervised by Ana Cristina Braga October 2020 DIREITOS DE AUTOR E CONDIC¸ ˜ OES DE UTILIZAC¸ ˜ A O D O TRABALHO POR TERCEIROS Este ´ e um trabalho acad ´ emico que pode ser utilizado por terceiros desde que respeitadas as regras e boas pr ´ aticas internacionalmente aceites, no que concerne aos direitos de autor e direitos conexos. Assim, o presente trabalho pode ser utilizado nos termos previstos na licen c¸ a abaixo indicada. Caso o utilizador necessite de permiss ˜ ao para poder fazer um uso do trabalho em condi c¸ ˜ oes n ˜ ao previstas no licenciamento indicado, dever ´ a contactar o autor, atrav ´ es do Reposit ´ oriUM da Universidade do Minho. Licenc¸a concedida aos utilizadores deste trabalho Atribuic¸˜ ao-N˜ aoComercial CC BY-NC https://creativecommons.org/licenses/by-nc/4.0/ i AGRADECIMENTOS Se o documento presente for julgado pelo leitor como sucinto e objetivo ser ´ a por esse mesmo impulso pessoal, de me manter breve e frio, que esta sec c¸˜ ao se torna a que mais tempo e reflex ˜ ao pede para a escrever, no que ser ´ a um eterno conflito, que adivinho partilhar com muitos colegas, do romance da minha l´ ıngua m˜ ae com a est´ eril realidade cient´ ıfica. Requer o estere ´ otipo de aqui catalogar fam ´ ılia e amigos para que possa agradecer o apoio cont ´ ınuo nos ´ ultimos anos. N ˜ ao posso deixar de rejeitar esta exig ˆ encia. O amor que sinto pelos mesmos ´ e melhor expresso diariamente e oralmente que num documento imut ´ avel por arte de escrita que de minha tristeza ser ´ a sempre inferior ` a eloqu ˆ encia da fala. Se porventura este gesto parecer radical direi preferir a compara c¸˜ ao aos antigos, dado ser (em parte) gra c¸ as ` a oralidade que mantemos at´ e hoje hist´ orias b´ ıblicas, filosofia socr´ atica e mitologia ass´ ıria. Libertando-me no entanto dos clich ´ es impostos de agradecer a pais e pa ´ ıs, em soberania de alma e em rigor telegr´ afico expresso um obrigado, ` A minha orientadora, Professora Doutora Ana Cristina Braga, por demonstrar um exemplo do humano e do acad ´ emico que individualmente merecem elogio, mas cuja conjuga c¸˜ ao exige louvor, Ao meu av ˆ o, Sr. Dr. Juiz Jos ´ e Machado e Costa, cuja eloqu ˆ encia de fala e escrita e procura da verdade tento replicar por diferentes caminhos, ` A minha av ´ o, Maria da Concei c¸˜ ao Amaral, a quem poderei entregar este documento como prova que as suas in´ umeras horas de ensino n˜ ao foram em v˜ ao. ii STATEMENT OF INTEGRITY I hereby declare having conducted this academic work with integrity. I confirm that I have not used plagiarism or any form of undue use of information or falsification of results along the process leading to its elaboration. I further declare that I have fully acknowledged the Code of Ethical Conduct of the University of Minho. iii RESUMO A curva ROC (Receiver operating characteristic) ´ e uma ferramenta anal ´ ıtica eficaz para testes cl ´ ınicos. A an ´ alise permite visualizar a varia c¸˜ ao de sensibilidade e especificidade para uma dada regi˜ ao de corte atrav´ es de um simples, mas robusto gr´ afico bidimensional. Num contexto biol ´ ogico, testes podem ser influenciados por m ´ ultiplas vari ´ aveis externas e como tal a analise ROC pode n ˜ ao ser a ideal ou gerar resultados incompletos. ´ E ent ˜ ao necess ´ ario saber que vari ´ aveis afetam determinado teste cl ´ ınico de forma a determinar os melhores par ˆ ametros para determinado teste ou at ´ e descartar determinada metodologia mediante a situa c¸˜ ao. O ajuste da curva ROC a covari ´ aveis permite a normaliza c¸˜ ao do efeito das mesmas ou diretamente ajustar a curva para os seus efeitos. Software direcionado ao ajuste da curva ROC ´ e, infelizmente, escasso e muitas vezes dif ´ ıcil de manusear por utilizadores n ˜ ao especializados. Recentemente o pacote AROC foi lan c¸ ado para R que disponibiliza v ´ arios recursos para estes ajustamentos, no entanto a dificuldade de utilizac¸˜ ao mant´ em-se. A combina c¸˜ ao deste pacote com a estrutura Shiny, um pacote que permite o desenvolvimento de aplica c¸ ˜ oes interativas, tem por objetivo a cria c¸˜ ao de um programa gr ´ atis e acess ´ ıvel que permita uma analise mais aprofundada dispon´ ıvel para todos os investigadores. RealROC foi capaz de replicar resultados de um caso de estudo que analisou a influ ˆ encia do sexo no sistema de pontua c¸˜ ao CRIB e respetiva previs ˜ ao de mortalidade, demonstrando a usabilidade e acessibilidade do programa que ser ´ a disponibilizado online e potencialmente contribuir para novos desenvolvimentos na ´ area. Palavras-chave: Curva ROC; AROC; Covari ´ aveis; Shiny; Bioestat ´ ıstica; Inform ´ atica M ´ edica; Classificac¸˜ ao estat´ ıstica; Software iv ABSTRACT Receiver operating characteristic (ROC) curves are a powerful analytical tool for clinical tests. The analysis allows the visualization of varying sensitivity and specificity for a given threshold through a simple, yet robust, two-dimensional plot. In a biological framework, tests can be influenced by multiple external variables, as such, standard ROC analysis may not be suitable or may provide incomplete data. It is then necessary to know which variables influence clinical test results to determine optimal conditions for trials or even to disregard a given method of evaluation in certain contexts. Adjusting for covariates allows ROC analysis to normalize the effects of the variable in question or to directly adjust the curve for its effects. Unfortunately ROC software that is able to conduct such an adjustment is sparse and proven difficult to use for non technical users. Recently, the AROC package for R was released and provides a robust resource for such adjustments however with he same usability problems previously stated. By combining this package with the Shiny framework, an R package that allows the creation of interactive applications, we hope to provide an accessible and free software that allows this extra depth of analysis to be available for all researchers. RealROC was able to mimic the results of a case study analysing the affects of sex to the CRIB score and resulting mortality rates that proving its practicality and will be made available online and hopefully contribute to the advancement of software in this field. Keywords: ROC curve; AROC; Covariates; Shiny; Biostatistics; Medical Informatics; Statistical classification; Software v CONTENTS Declarac¸˜ ao de Direitos de Autor i Agradecimentos ii Statement of Integrity iii Resumo iv Abstract v 1 introduction 2 1.1Context and Motivation 2 1.2Goals 3 1.3Document Organization 3 2 state of the art 4 2.1ROC Curve Concept 4 2.1.1ROC curve and applications 4 2.1.2Summary Indexes 5 2.1.3Binormal Model 7 2.1.4ROC curve for Ordinal Tests 8 2.2ROC Curve Estimation 9 2.2.1Empirical Estimation 9 2.2.2Modeling the ROC curve 12 2.3Bayesian Methods 14 2.3.1Bayesian Approach and General ROC Analysis 14 2.3.2Parametric and Non-parametric Bayesian Methods 16 2.4ROC Curve and Covariates 19 2.4.1The Need to Adjust for Covariates 19 2.4.2Adjusting the ROC curve 19 2.4.3AROC and Summary Indexes 22 3 software and packages 24 vi contents vii 3.1Introduction 24 3.2The AROC R Package 25 3.2.1Package Overview 25 3.2.2Bayesian Methods 25 3.2.3Kernel Methods 26 3.2.4Frequentist Method 26 3.2.5Other Methods 27 3.3R Shiny 28 3.4Comp2ROC 28 4 realroc web application 30 4.1Application Flow 30 4.2Application Modules 31 4.2.1Home Screen 31 4.2.2Classic ROC Analysis 32 4.2.3AROC 33 4.2.4Comparing adjustment 34 4.2.5Report Screen 35 4.2.6Advanced Options 36 4.2.7Help Section 36 5 case study -crib score and the sex covariate 37 5.1Neonatal Mortality and CRIB Score 37 5.2Dataset 38 5.3RealROC vs scripting 39 5.4Final remarks 44 6 conclusion 45 2.1. ROC Curve Concept 6 selection of results in indicated by the chance diagonal (area of the triangle of base and height 1) giving us the formal definition of AUC as AUC =Z1 0y(x)dx (5) or, by using equation 3 AUC =Z1 0ROC(t)dt. (6) A simple and commonly used interpretation of AUC is as an average value of sensitivity for all possible values of specificity taking values between 0and 1. An AUC value closest to 1indicates a test with 100% accuracy bringing its practical lower limit to 0.5or 50% accuracy which we refer to as the chance diagonal indicating a test relies on luck and is, fundamentally, not suitable (Park et al., 2004). A more formal interpretation states the AUC is equal to the probability that tests results from a randomly selected pair of diseased and non-diseased subjects are correctly ordered, i.e. P[YD>Y¯ D]. Often times a particular TPF is of interest, this is particularly the case in medical contexts where the FPF is particularly small ( <0.05 ), in these cases and for specifying a range of threshold values Partial Area Under the Curve (pAUC) is used. Working with equation 6, we define PAUC as, pAUC(t0) = Zt0 0ROC(t)dt. (7) Both AUC and pAUC are used regularly, however a few others deserve mentioning. The Youden Index (YI) is the maximum difference between TP and FP fractions, YI =max(tp −f p) = max(tp +tn −1), (8) The threshold at the point on the ROC curve corresponding to YI is often taken to be the optimal classification threshold. Another such summary index is the Max Vertical Distance (MVD) between the chance diagonal and the ROC curve, MVD =max|y(x)−x| a 2.1. ROC Curve Concept 7 particularly useful statistic for its equivalence to both YI and Kolmogorov-Smirnov statistic in the ROC curve domain (Krzanowski and Hand, 2009). 2.1.3Binormal Model In ROC analysis, the binormal model refers to the assumption of normal distributions of both populations. This is the cornerstone of ROC analysis and a standard by which other specialized analysis can be judged. Beginning with the aforementioned assumption YD∼N(µD,σ2 D),Y¯ D∼N(µ¯ D,σ2 ¯ D), then ROC(t) = Φ(a+bΦ−1(t)) (9) where Φ(·)is the normal cumulative distribution function (cdf) and, a=µD−µ¯ D σD,b=σ¯ D σD We see the definition of the binormal ROC curve in equation 9where we call a the intercept and b the slope for the curve where both values are positive if we abide by convention of larger values being indicative of disease (Pepe, 2003). Krzanowski and Hand (2009) further expanded this equation form y(t) = ΦµD−µ¯ D+σ¯ D×zt σD(10) where, zt=Φ−1[t(c)] = µ¯ D−c σ¯ D For the binormal ROC curve the AUC is AUC =Φa √1+b2, (11) an increasing function of a and decreasing of b . Partial AUC is not derivable from the previous expression and must be calculated using numerical integration or rational polynomial approximation. 2.1. ROC Curve Concept 8 As mentioned, the ROC curve is invariant to monotone increasing data transformations therefore to say YD and Y¯ D is to say that to for a strictly increasing transformation h , h(YD) and h(Y¯ D) have normal distributions, we can further impose that this function h transforms the data to normality. The assumption that such a function h makes both populations normally distributed exists is a weak one however empirical testing shows that for nonGaussian distributions the binormal model still holds (Pepe, 2003). This emphasizes the concept that ROC curves quantify relationships between distributions and have no link to a particular one. It is this standard model we will see adapted to allow covariate adjustment further ahead. 2.1.4ROC curve for Ordinal Tests So far, theoretical assessments have implied a numerical, continuous scale of data, however there are many cases where tests are not just discrete variables but are non numerical all together with subjective assessments that different assessors may ”grade” differently. This has been an issue since the ROC curve’s infancy and is a recurring problem in Radiology where, for instance, different categories (0-5) of assessments are given for the same image by different radiologists (Fig. 1). A common framework to tackle this problem is the latent variable framework. Say, L is an unobserved latent continuous variable corresponding to the assessor’s perception of the image. Much like mathematical thresholds, this assessor has its own, subjective, cut off points used to classify/rate the image. Y corresponds to the reported classifications, we say, Y=y⇐⇒ cy−1<L<cy,y=1, ...., P where c0=−∞ and cp=∞ . The reader classifies the image in the yth category if L falls within the interval corresponding to his implicit definition for the yth category (cy−1,cy). To model the ROC curve we use L , the latent decision variable, while not an observable variable we are able to identify P+1 points on the basis of Y and by interpolation the ROC curve for L. Since, Y≥y⇐⇒ L>cy−1 2.2. ROC Curve Estimation 9 Figure 1: Example of rater bias in ordinal tests where rater 2shows a more conservative threshold than rater 1for c (Pepe, 2003). We can identify TPF and FPF corresponding to the threshold cy−1 as P[Y≥y|D=1] and P[Y≥y|D=0] respectively (Pepe, 2003). If we are able to identify two non-degenerate points t1and t2we can apply the binormal model where aand bare given by a=Φ−1(ROC(t1)) −bΦ−1(t1) and b=Φ−1(ROC(t2)−Φ−1(ROC(t1)) Φ−1(t2)−Φ−1(t1). We will see this framework mentioned in the next subsection. 2.2 roc curve estimation 2.2.1Empirical Estimation Having explored the theoretical fundamentals of the ROC curve and its summary indexes we now turn to statistical methodology for inferring this curve from existing data rather than assuming the existence of sets of populations and classifiers that fit the mentioned criteria. Three distinct approaches can be considered for this estimation: 2.2. ROC Curve Estimation 10 1. Apply non parametric empirical methods to the data to obtain the empirical ROC curve, from which empirical summary indexes can be calculated; 2. Use statistical models for the distributions of cases and controls, parameters in these distributions are estimated and both induced ROC curve and summary indexes are calculated together; 3. Modeling the ROC curve, rather than the probability distributions, as a smooth parametric function. All three approaches have their respective strengths and drawbacks however in this subsection we will be focusing on the empirical methods followed by the modeling option in subsection 2.2.2for being the most popular when working with continuous and ordinal data respectively and for mirroring methods of covariate adjustment we will see in the following sections. Of course when dealing with statistical estimation, one must take into account the accuracy and precision of these estimates, as well as estimator biases and sampling variability to construct reliable Confidence Intervals (CI), this will be explored along with the estimation procedure for both ROC curve and AUC. Empirical estimation applies the ROC curve definition to the observed data thus empirical TPF and FPF are calculated as d TPF(c) = nD ∑ i=1 I[YDi≥c]/nD(12) d FPF(c) = n¯ D ∑ j=1 I[Y¯ Dj≥c]/n¯ D(13) The empirical ROC curve, d ROCe is a plot of d TPF(c) versus d FPF(c) for all c∈[−∞,∞] writable as, d ROCe(t) = b SD(b S−1 ¯ D(t)), (14) where b SD and b S¯ D are the empirical survivor functions of YD and Y¯ D respectively. An example of the resulting curve can be seen in Figure 2. 2.2. ROC Curve Estimation 11 Figure 2: Example of an empirical ROC curve. Sampling Variability and CIs Several approaches are available for sampling variability in the empirical ROC curve ranging from fixing threshold to fixing FP and TP fractions as well as considering the entire curve with confidence bands. For our purposes we will consider the first and latter options. Fixing the threshold, c , and calculating a joint confidence region using exact binomial or asymptotic methods or, alternatively, with bootstrap resampling if the samples are dependent, this method is particularly useful with knowledge of the threshold in advance, such as in blood tests where standard thresholds are already defined. Considering the entire curve using a confidence band is particularly applicable when no assumptions can be made regarding the test and is more useful for describing the curve. One approach to calculate these bands is to base the calculation on the distribution of sup|ROCe(t)−ROC(t)| , which can be calculated using two independent Brownian bridges or, alternatively, model the risk function P[D=1|Y] with logistic regression methods (Pepe, 2003). Empirical AUC As mentioned empirical indexes can be applied directly to the empirical curve thus, d AUCe=Z1 0d ROCe(t)dt (15) 2.2. ROC Curve Estimation 12 pd AUCe(t0) = Zt0 0d ROCe(t)dt. The area under the empirical ROC curve is the Mann-Whitney U-statistic (Pepe, 2003). Variability of AUCe Calculating variability in AUCe or any other summary index is often complicated, involving analytic expressions for asymptotic variance however in practice one simply uses bootstrap methods to calculate CIs. 2.2.2Modeling the ROC curve In the beginning of this section we described using a smooth parametric function to model the ROC curve as a popular option when working with ordinal data, however with a few adjustments we can also work with continuous data and we will describe both methods in the remainder of this chapter. Ordinal Tests For this approach we adopt the Latent Decision Variable (LDV) framework mentioned in section 2.1.4. Let L denote the underlying decision variable for a single reader, recall the binormal model of the ROC curve is ROC(t) = Φ(a+bΦ−1(t)) as seen in equation 9, we assume L∼N(0, 1) in the non diseased population, ¯ D . With these conditions we imply that in the diseased population D , L∼N(a/b,(1, b)2 , and calculate a and b using the already established relation of Y=y⇐⇒ cy−1<L<cy . We can then derive the probability of diseased and non diseased observation YDand Y¯ D P[Yb D=y] = Φ(bcy−a)−Φ(bcy−1−a). (16) P[Yb¯ D=y] = Φ(cy)−Φ(cy−1)(17) 2.2. ROC Curve Estimation 13 The log likelihood function is then constructed as nD ∑ i=1 logP[YD=YDi] + n¯ D ∑ j=1 logP[Y¯ D=Y¯ Dj](18) and maximize with respect to parameters {a,b,c1....cP−1} . The resulting ROC curve is a smooth curve that follows the binormal model replacing a and b with it’s estimators, b a and b b. The Standard Error (se) for b a+b bΦ−1(t)is se =var(b a) + (Φ−1(t))2var(b b) + 2Φ−1(t)cov(b a,b b)1/2 . (19) Corresponding confidence limits for a+bΦ−1(t)written as b a+b bΦ−1±Φ−1(1−α)se, (20) generate the confidence limits for the binormal ROC curve. Continuous Tests Metz et al. (1998) proposed dealing with continuous data by categorizing them into a finite number of pre-defined categories and applying the previously explained fitting methods for ordinal data. The asymptotic properties of the estimator require the categories to be predefined, however (Metz et al., 1998) proposed defining theses categories using observable data using vertical and horizontal jumps in the empirical ROC to define the categories. Having defined the categories we apply the LDV framework. This method is called LABROC and it’s most appealing feature is, by working with ranks, being distribution free, making it invariant to monotone increasing data transformations mirroring the ROC curve (Pepe, 2003; Metz et al., 1998). 2.3. Bayesian Methods 14 2.3 bayesian methods 2.3.1Bayesian Approach and General ROC Analysis In the previous sections we discussed what is known as the frequentist or classical framework for statistical inference. In this framework populations are represented by probability models whose parameters are treated as fixed but unknown quantities about which inferences are to be made and, thus has no scope for extraneous information integration such as previous experiments or subjective elements. The Bayesian approach treats population parameters as random variables with probability distributions reflecting the degree of belief the research has on the data, allowing the introduction and combination of prior knowledge and subjectivity to it. These prior distributions are combined with sample data to produce posterior distributions creating a basis for all inferences. Bayesian methods have been available for many years however they were hampered by their intractability, general unattractiveness by statistics and bottlenecks in computer processing, this however changed by the mid 90’s with the introduction of Markov-chain Monte Carlo methods (Krzanowski and Hand, 2009). This method generates sample values from the posterior distributions and approximates integrals by the average of sample values of a function f(·) , this features ensures each new proposal depends on the current one, and the sequence of proposals is guaranteed to converge to values from the desired distribution. Noticeable algorithms include MetropolisHastings, which uses joint distribution of new and current proposals and Gibbs sampling which uses a sequence of conditional distributions. We will see the latter mentioned in the next chapter. Bayesian methodology shines when studies have uncertainty about underlying quantities and ROC analysis has benefited greatly for this addition to the statistical repertoire. A notorious case of this uncertainty is in labeling as D or ¯ D of subjects from which a ROC curve is to be constructed is either fallible or not available, this is often the case in medical studies when the true disease status of each sample member is either equivocal or unknown and ways. 2.3. Bayesian Methods 15 This section will explore the fundamentals of Bayesian ROC analysis. For continuous data very few proposals for a standard method exist, Erkanli et al. (2006) suggests the binormal model, presented in section 2.1.3, is a poor model for analysis in Bayesian methodology for its pretense of normal distributed populations using a mixture of normals as a more flexible alternative. Ignoring for the time being population distinctions ( D and ¯ D ), say C denotes the number of components in the normal mixture and K is a random variable which indicates the operative component for a classification score, Y. The formal Bayesian approach is then, Y|K,θK∼N(µK,σ2 K)(21) where K and θk= (µK,σ2 K) are parameters and must be assigned prior distributions. The components of θK are assigned normal gamma baseline priors while K has an independent C-state multinomial distribution with probabilities w1, ..., wCspecified as, w1=R1,wk= (1−R1)(1−R2)...(1−Rk−1)Rkf or k =2, ..., C where Ri are independent Beta(1, α ) variables and RC=1 to ensure that the wi sum to 1. This model is termed a mixed Dirichlet process. With this model we can generate predicted values for Y given previous scores a density function, however the resulting integral is difficult to compute even with a small C , as such a Gibbs sampler is used to simulate observations and the expected value can be approximated by the average of these samples. Having approximated the posterior predictive density we then obtain cumulative distribution function from which the posterior predicted true positive and false positive rate are calculated, from there, varying the threshold yields the predictive ROC curve. Krzanowski and Hand (2009) dedicate a chapter to this model and its applications where this section was based on, for the sake of brevity some equations and proofs were excluded however this is a fundamental companion source of information to the original material of Erkanli et al. (2006) and should be consulted for further information. 2.4. ROC Curve and Covariates 22 The need for this flexibility and constraints has popularized the use of the generalized linear models (GLM) (McCullagh and Nelder, 1989), specifically the ROC - GLM model given by, h(y) = b(x) + βTX(28) where, •b(·) is an unknown baseline function monotonic on (0,1); •h(·) is the link function, specified as part of the model and also monotonic on (0,1); •X is the vector of covariates (however this time it’s a single vector rather than associated to a specific population); •βare the regression parameters associated with the covariates. Link functions, h(·), vary but the more common in GLM are •Probit: h(y) = Φ−1(y) •Logistic: h(y) = log y 1−y=logit(y) •Logarithmic: h(y) = log(y) The same methods of obtaining the mean of both populations can be as applied in equations 24 and 25 (Pepe, 2003) by replacing the individual covariates vectors, XD and X¯ D by X as well as and respective σD , σ¯ D standard deviations which speaks to the flexibility of the approach (Krzanowski and Hand, 2009). 2.4.3AROC and Summary Indexes In some cases, covariate adjustment interest lies on the summary statistics rather than the ROC curve itself where most of the theoretical attention is focused on the AUC. 2.4. ROC Curve and Covariates 23 If the ROC curve was adjusted indirectly then the effects on AUC are obtained by simply expressing AUC in terms of the model parameters and substituting estimates of the parameters in the binormal expression for xD and x¯ D for covariates XD and X¯ D respectively: AUC(xD,x¯ D) = φ[δ(xD,x¯ D)], (29) where δ(xD,x¯ D) = µD(xD)−µ¯ D(x¯ D) √σ2 D+σ2 ¯ D for µD(xD) = αD+βT DxD µ¯ D(x¯ D) = α¯ D+βT ¯ Dx¯ D. Estimation of these parameters by ordinary least squares and substituting into these expressions yields b δ(xD,x¯ D) and d AUC(xD,x¯ D) . This modulation will generally require strong or specific assumptions. This issue as well as specific cases where AUC values were available but not classification scores, prompted several authors to consider direct adjustment AUC. Dodd and Pepe (2003) adapted an early proposed model of A AUC calculation that would allow all types of covariates to be modeled with a binary estimation method. Once more, suppose there is a set of cD covariates XD associated with population D and a corresponding set c¯ D of X¯ D for population ¯ D and we write X as the aggregation of all covariates. We have the AUC regression model, E(h[d AUC]) = α+βTx(30) for parameters α , β and monotone increasing link funtion h(·) , with probit and logit again as natural link functions (Krzanowski and Hand, 2009). 3 SOFTWARE AND PACKAGES 3.1 introduction A large amount of software has been developed for ROC curve building and analysis across multiple platforms, enumerated in detail by a recent review by Obuchowski and Bullen (2018), some notable references include ”Metz ROC Software”, developed by the University of Chicago, ”Pepe Lab” for the Stata platform, ”Analysing ROC curves with SAS”, developed by Mithat G ¨ onen for SAS and the ”pROC” package for R. These packages while popular either offer a somewhat shallow approach to covariate adjustment or are published in a non readily available platform that can prove difficult to use for non technical users. To fulfill our goal of creating a user friendly application and bridging the gap between theoretical and practical application of the AROC in a widely available space the Shiny framework was selected. Shiny uses the R language to create an interactable application that can either be used on a local computer or online. However a framework in R requires packages to perform ROC adjustment. Recently, a package published by Rodr ´ ıguez- ´ Alvarez et al. (2018) on CRAN, the R Archive Network, goes to great lengths to implement several regression approaches for the inclusion of covariate information on the traditional ROC framework, however, it still suffers from a previously mentioned problem of difficulty in usability. By combining an interactable framework with this robust package for ROC curve adjustment the feasibility of our goal becomes clearer. Additionally, to provide the user an extra option of comparison between ROC curves to pre and post adjustment the Comp2ROC package was employed. 24 3.2. The AROC R Package 25 In this chapter we will take an indepth look at the AROC Package by Rodriguez-Alvarez and Inacio de Carvalho (2018), the Shiny framework as well as a quick overview of the Comp2ROC package by Braga et al. (2016). 3.2 the aroc r package 3.2.1Package Overview The AROC R-package developed by Rodr ´ ıguez- ´ Alvarez et al. (2018) implements different methods of computing covariate information in ROC curve construction and two additional methods of calculating marginal/polled ROC curve providing aditional tools for an underdeveloped sub-field. The methods include Bayesian, Kernel and Frequentist methods of adjustment and direct and indirect methods of regression, mentioned in the previous chapter, we will explore in the following subsections. 3.2.2Bayesian Methods The first methods presented by the package are two bayesian based estimations based on the Dirichlet process discussed in section 2.3. The AROC.bnp function estimates the AROC curve using the Bayesian nonparametric method. Working with equation 4we say, AROC(t) = Pr{YD>F−1 ¯ D(1−t|XD)} =Pr{1−F¯ D(YD|XD)) ≤t} =Pr(UD≤t), 0 ≤t≤1, (31) where UD=1−F¯ D(YD|XD) , a placement value of the test outcome in the diseased population or the the standardization of YD to the conditional distribution of Y¯ D making the AROC a cumulative distribution function of UD . This method first models the conditional distribution of test outcomes in the nondiseased group, F¯ D using a B-splines dependent Dirichlet process mixture of normals model followed by modeling UD and it’s cumulative distribution 3.2. The AROC R Package 26 using a non parametric regression model through Bayesian bootstrap (Rodr ´ ıguez- ´ Alvarez et al., 2018). The AROC.bsp method is very similar to the previous, non parametric method, in both construction and theory, where this models F¯ D using a normal linear regression, making it a counterpart to Janes and Pepe (2009) frequentist model and its AROC package method - AROC.sp (Rodriguez-Alvarez and Inacio de Carvalho, 2018). 3.2.3Kernel Methods The AROC.kernel, an earlier proposed model by Rodriguez-Alvarez et al (Rodr ´ ıguez- ´ Alvarez et al., 2011) is, as the name implies, a kernel based method. Test outcomes for the non diseased group are modeled with a location-scale regression model where both regression and variance functions are estimated using Nadaraa-Watson local estimators that in turn are used to compute standardised residuals to model UD and estimate the AROC curve, d AROC(t) (Rodr ´ ıguez- ´ Alvarez et al., 2018; Rodriguez-Alvarez and Inacio de Carvalho, 2018; Rodr´ ıguez- ´ Alvarez et al., 2011). The package author notes that, for now, this method, unlike the previous, can only handle a single continuous covariate. 3.2.4Frequentist Method This semiparametric frequentist method arguments and construction are fairly similar to the previous methods but, as previously hinted, uses a semiparametric location regression model for Y¯ D to estimate F¯ D and estimates outer probability empirically (Rodriguez-Alvarez and Inacio de Carvalho, 2018; Janes and Pepe, 2008) making it less computationally heavy, i.e faster in comparison. Another advantage is being able to provide direct insight on covariate influence over the test/marker with the fit model’s parameters. The frequentist method for covariate adjustment was first implemented in rocreg from Stata 2013 complementing other Stata ROC commands such as roctab and roccomp (StataCorp, 2013). It arguably remained one of the most in depth options for covariate adjustment available in the market until the release of the AROC package, it is however only available 3.2. The AROC R Package 27 through Stata which requires an yearly licence and can prove challenging for inexperienced users. 3.2.5Other Methods Posterior Predictive Checks (PPC) All previous methods are meant to construct and analyze the AROC curve and the behavior of the incorporated covariates, for the remainder of this section the methods focus on posterior predictive checks and pooled ROC estimation. Both predictive.checks.AROC.bnp and predictive.checks.AROC.bsp are implementations of PPCs on their respective Bayesian based method. The premise behind PPCs is evaluating the generated model on how well it is able to generate data similar to the data observed (Gabry et al., 2019) utilizing in this case the B-splines dependent Dirichlet process and Bayesian normal linear regression model for the AROC.bnp and AROC.bsp objects respectively (Rodriguez-Alvarez and Inacio de Carvalho, 2018). To exemplify these methods we use the previously generated objects. While the mandatory argument is solely the AROC object, the devnew was changed from the default TRUE argument to display all depicted graphics on the same device. These graphics are histograms of the desired test statistics, by default minimum, maximum, median and skewness and Kernel density estimates showing diagnostic test outcome in the nondiseased group as well. Pooled ROC The package also provides two methods of polled ROC estimation, pooledROC.emp an implementation of empirical estimation proposed by Hsieh and Turnbull (Hsieh and Turnbull, 1996; Rodriguez-Alvarez and Inacio de Carvalho, 2018) and pooledROC.BB for Bayesian bootstrap estimation proposed by Gu et al. (2008) (Rodriguez-Alvarez and Inacio de Carvalho, 2018). Both methods are similar to construct and have a similar output structure, mandatory 3.3. R Shiny 28 arguments identify diseased and non diseased groups in the test or marker column which can be achieved with a straightforward indexing in R. 3.3 r shiny Shiny is an R package which provides a framework to develop interactive web-based applications such as data summaries and queries to end users through a standard web browser. A recent addition to the R repository, its applications have provided an excellent resource for users to work with complex R packages or perform data analysis through a extensible visual framework based around HTML and CSS with further JavaScript and jQuery integration to extend the scope of possible applications. The package is based around reactive programming, a programming paradigm that facilitates the automatic propagation of change of dataflows. A practical understanding of this concept is, while working with several inputs generating a specific output, be it plots, tables or text, any modification on input will automatically generate and update the output without requiring a new command or refresh action on the user’s side. Another user friendly implementation is the customization of the shiny interface using widgets and code chunks promoted by RStudio and implemented across several other apps available on their archive The standard architecture of a shiny app is two scripts in the same directory, ui.R for the User interface (UI) and server.R for internal calculations and app behavior. These two files can also be joined in a single app.R however this option is best used for applications such as data visualization (Beeley, 2013). 3.4 comp2roc The final package required to build the RealROC app is Comp2ROC. To establish a statistical significant confounding affect of a given covariate on the ROC curve we have explored the possibility of using the fit model’s parameters as seen in section 3.2.4however another well established alternative is available in some scenarios, ROC curve comparison. When comparing ROC curves, Area Under the Curve (AUC) is frequently used to indicate greater performance. The Z-test, a non parametric approximation of the Wilcoxon-Mann- 3.4. Comp2ROC 29 Whitney test is often used for this comparison, however performance misreadings can occur due to curve crossing, further argued is the ROC curves inclusion of areas with very little interest that can influence the result perchance claiming that, for instance, two ROC curves have no statistical differences overall disregarding partial areas of interest (Braga et al., 2013). For a more robust comparison of the two ROC curves the Braga methodology (Braga et al., 2013) and package Comp2ROC from the same author was selected. The method uses a collection of sampling lines similar to multi-objective distinct optimization algorithms allowing ROC curve comparison in several regions of space evaluating statistical performance and generating confidence intervals with non parametric bootstrap re-sampling method Braga et al. (2013). Due to compatibility issues we cannot use this package to compare curves with and without accounting for a covariate affect however in instances where we can separate the populations assigned with a given covariate value, such as in the case of a categorical variable, this method can be used to calculate the differences between curves, inferring the covariate affect which, if used in tandem with the AROC method, can give further credence to the conclusions drawn when using RealROC. This dissertation did not thoroughly explore ROC curve comparison however both the package original article (Braga et al., 2013) and documentation (Braga et al., 2016) as well as Krzanowski and Hand (2009) provide valuable resources on this topic. 4 REALROC WEB APPLICATION 4.1 application flow RealROC follows a linear flow intended to replicate an analysis progression from building the standard ROC curve to adjusting it followed by a comparison of both and a report with all relevant statistics that can be consulted at any time. The user is also able to skip or return to any of the previous steps to change any parameters as needed. The app is comprised of seven distinct modules that will be individually explored in the next section. Figure 4: RealROC flowchart. 30 4.2. Application Modules 31 The application is available online at https://frmachadoecosta.shinyapps.io/RealROC/ , additionally it can be downloaded or run locally by downloading it at https://github.com/ frmachadoecosta/RealROC where the source code is also available. 4.2 application modules 4.2.1Home Screen When RealROC opens, the user is presented with an option of uploading the desired data, in csv or xml format, or using the sample data provided by the app. The sample data, endosim.csv, is present in another package developed by Rodriguez-Alvarez and Javier RocaPardinas (2017) and used Body Mass Index (BMI) to detect patients having a higher risk of cardiovascular problems, with age and biological sex as covariates. This data was selected not just because BMI is a well established index, making both values and results easier to understand but also for the presence of both a continuous and binary variables, age and sex respectively, to allow for different approaches when adjusting the ROC curve. Figure 5: Home Screen options flowchart. Either option will trigger the appearance of a data-table where the user can consult their data and adjust the options for import if necessary as can be seen in Figure 6. 5.2. Dataset 38 base deficit, minimum and maximum appropriate fraction of inspired oxygen (Fi O2 ), all measured in the first 12 hours) resulting in a score between 0and 24 where higher values denote a higher risk of death. This scoring system is a staple in NICUs for both statistical and practical reasons, being not just an overall better scoring system than most of its contemporaries (Braga et al., 2013; Bastos et al., 1997) but also for the ease of data collection and score calculation making the test last only 5minutes per infant (Dorling et al., 2005). The CRIB score was derived using data from infants admitted to UK tertiary neonatal units from 1988 to 1990, this lead to Parry et al. (2003) raising concerns over poor calibration to contemporary data, the inclusion of Fi O2 also warranted criticism since it was considered a subjective entry determined by the care team rather than a physiological measurement. Similarly the addition of data up to 12h after admission also lead to concern over early treatment bias (Parry et al., 2003). These issues eventually resulted in an updated CRIB score, CRIB-II. CRIB-II score system maintains a relatively low number of variables needed for calculation, maintaining the advantage of its predecessor, gestational age, birth weight, admission temperature and base excess are used to predict mortality. This new prediction tool was met with some skepticism with reports of studies comparing both systems, showing no statistical difference (Gagliardi et al., 2004; Felice et al., 2005) nevertheless CRIB-II is now a known and recognized neonatal scoring system used in some hospitals. 5.2 dataset The dataset used is part of the Portuguese National Registry on low weight newborns between 2013 and 2018 and was made available for research purposes. The original data included possible confounding observations such as repeated ids and twins that were removed to ensure no unaccounted variable. After ensuring all id’s were unique, these were promptly removed along with any possible identifiable features to abide by European Union’s anonymity and data protection standards. The resulting dataset composed of 3823 unique entries registering gestational age in weeks, the mothers age, in years, biological sex of the infant (1-Male; 2-Female), CRIB score (0-21), survival (0-Survival; 1-Death) and other possibly relevant covariates were used for the remainder of the study. An abridged dataset 5.3. RealROC vs scripting 39 with relevant information was published and is available for consultation (Machado E Costa, 2019). 5.3 realroc vs scripting In both methodologies the dataset must be imported to the environment, this is achieved by choosing a file in the application of by using the read.csv() function in R, note that server-side RealROC is performing this same command however it simplifies the process and replaces R syntax with simple button presses. A successful data import can be seen in Figure 12 with a table output where the user can consult their data unlike in native R. Figure 12: Successful data import with table output. Moving to the Classic ROC section, we can begin selecting the parameters and ploting several different ROC curves, in R this requires loading and potentially installing the desired libraries however these come preloaded in the application. Once again the need to know each package correct syntax is replaced by simple data inputs that can display the empirical or pooled ROC curve as seen in Figure 13. 5.3. RealROC vs scripting 40 Figure 13: Pooled empirical ROC curve for CRIB score. Extra information is even present in the ”Population Distribution” tab displayed in Figure 14 where we can see the population densities for controls and cases, giving a broader view of the data. Figure 14: Population densities of mortality by CRIB score. Moving to the AROC section, the user will see all previous selected inputs are saved in memory and all that is needed is to correctly select the covariate, in this case ”sex”. The curve type or method of adjustment selected is the frequentist method mentioned in section 3.2.4for a more in depth summary on the effect of the covariate. While the AROC curve can 5.3. RealROC vs scripting 41 be seen, the population distribution tab will not display the correct result, this is because the data was originally submitted with the ”sex” column as numeric rather than a factor, this however can be easily fixed in the Advanced section where the nature of the column can be changed. Results for both outputs can be seen in Figures 15 and 16. Figure 15: AROC curve for CRIB score adjusted for sex. Figure 16: Boxplot of covariate specific CRIB controls and cases. 5.3. RealROC vs scripting 42 Already we can see the behaviour of the covariate and might opt to go straight to the Report section to get our summary however we can go further and directly compare the two curves in the Comparison section. Given the nature of how covariate we can use both AROC and Comp2ROC methods mentioned in Chapter 3. In the Comparison section the gap of usability between app and script is far more noticeable, the AROC method has no method for superimpose two curves generated by the package so for scripting the user would require extensive knowledge of R plot syntax to be able to manually collect each curve points and plot them. This, like other features of the application, is achieved with simple button presses after selection of the method of comparison, displaying a ggplot object with CI present as can be seen in Figure 17. For the Comp2ROC method, the package requires a particular data structure to operate, and the user would need to manually rearrange the dataset and re-import it, in RealROC however this is seamless to the previous method and results can be seen in Figure 18. Figure 17: AROC and ROC curve comparison. With the data thoroughly analysed, the user can now see the summary of their findings in the Report section that can be seen in Figure 19. This section displays several summaries for the many operations made, and a equivalent in the script method would be several different summary commands. The section organizes the many summaries offered by each package in an easily digestible format containing all relevant information. 5.3. RealROC vs scripting 43 Figure 18: Comp2ROC comparison of covariate specific ROC curves. Figure 19: Report section for sex adjusted CRIB. In this instance we can see that the AROC method shows a p-value >0.05 which coupled with the Z statistic p-value >0.05 and the sum of global areas crossing 0 along with the CI intersections clearly shows no statistical significance to the sex covariate in the CRIB score system. 5.4. Final remarks 44 5.4 final remarks Results from both covariate adjustment and ROC curve comparison methods clearly indicate sex is not a confounding covariate in the CRIB score meaning the sex of the infant has no statistical relevance to the score system or the mortality outcome. These findings validate previous research on the matter (Terzic and Helji ´ c, 2012; Mour ˜ ao et al., 2014) from more contemporary sources but also come from a substantially larger pool of data. RealROC displays an intuitive option selection menu and clear output which allows any user to perform a thorough analysis without R code syntax details. It is, of course, worth noting that all results displayed by the app can be achieved with R by an experienced user however the time to do so would be substantially longer than using the app. 6 CONCLUSION ROC curves have aided statistical analysis and decision statistics for over half a century, their simple classification model and straightforward summary statistics have been paramount in the virtually limitless two-class prediction problem. The curve’s history is filled with new coefficients, derivations, summary statistics and interpretations for existing theory to better fit such an ubiquitous tool to specific applications. Covariate specific and covariate adjusted ROC curves are a relatively newer addition to the existing theory hence the lack of specific tools mentioned in section 3.1. This addition however widens the scope for this tool even further with works like the ones from Janes and Pepe (2008) and Rodr´ ıguez- ´ Alvarez et al. (2018) only serving to expand its potential. The goal for this dissertation was to cement these new found applications for the ROC curve via the development of a specialize tool. RealROC was developed as a way to introduce users to the AROC concept and apply it to their database and hopefully lead them to more in depth conclusions about their data and the relations between variables. While this dissertation focused on the bioinformatics and medical applications for this tool, RealROC is built to allow the same diversity of data the classic ROC analysis is able to compute. The case study presented in the previous chapter demonstrated an intuitive software that simplifies the AROC analysis and allows users to perform their intended studies in a shorter amount of time and effort, and while undeniable that a pure R code can provide greater freedom to an experienced user in customization and specific calculations the app greatly reduces the know-how entry barrier. 45 46 The application is already released and available at https://frmachadoecosta.shinyapps.io/ RealROC/ however new improvements are expected such as introducing several covariates to the confounding study as well as specifying the relation between covariates and between each covariate and the marker, introducing the ability to compare markers in the Comparison module, new additions to the Advanced tab to allow smoother data manipulation inside the app and any and all performance and bug fixes reported by users through github. BIBLIOGRAPHY Bastos, G., Gomes, A., Oliveira, P., and Da Silva, A. T. (1997). Compara c¸˜ ao de quatro escalas de avalia c¸˜ ao da gravidade cl ´ ınic a (CRIB, SNAP, SNAP-PE, NTISS) em rec ´ em nascidos prematuros. Acta Medica Portuguesa,10(2-3):161–165. Beeley, C. (2013). Web application development with R using Shiny. Packt Publishing. Bradley, A. P. (1997). The use of the area under the ROC curve in the evaluation of machine learning algorithms. Pattern Recognition,30(7):1145–1159. Braga, A. C., Costa, L., and Oliveira, P. (2013). An alternative method for global and partial comparison of two diagnostic systems based on roc curves. Journal of Statistical Computation and Simulation,83(2):307–325. Braga, A. C., Frade, H., Carvalho, S., and Santiago, A. M. (2016). Comp2ROC: Compare Two ROC Curves that Intersect. R package version 1.1.4. Brito, A. S., Matsuo, T., Gonzalez, M. R. C., de Carvalho, A. B. R., and Ferrari, L. S. L. (2003). CRIB score, birth weight and gestational age in neonatal mortality risk evaluation. Revista de sa´ ude p´ ublica,37(5):597–602. Choi, Y.-K., Johnson, W. O., Collins, M. T., and Gardner, I. A. (2006). Bayesian inferences for receiver operating characteristic curves in the absence of a gold standard. Journal of Agricultural, Biological, and Environmental Statistics,11(2):210–229. Dodd, L. E. and Pepe, M. S. (2003). Semiparametric regression for the area under the receiver operating characteristic curve. Journal of the American Statistical Association,98(462):409–417. Dorling, J. S., Field, D. J., and Manktelow, B. (2005). Neonatal diseases severity scoring systems. Archives of Disease in Childhood: Fetal and Neonatal Edition,90(1):11–16. Egan, J. P. (1975). Signal detection theory and ROC-analysis, volume 1. ISBN 0122328507. Erkanli, A., Sung, M., Jane Costello, E., and Angold, A. (2006). Bayesian semi-parametric roc analysis. Statistics in Medicine,25(22):3905–3928. Faraggi, D. (2003). Adjusting receiver operating characteristic curves and related indices for covariates. Journal of the Royal Statistical Society: Series D (The Statistician),52(2):179–192. 47