auditor: an R Package for Model-Agnostic Visual Validation and Diagnostic

Machine learning has spread to almost every area of life. It is successfully applied in biology, medicine, finance, physics, and other fields. The problem arises if models fail when confronted with the real-world data. Therefore, there is a need for validation methods. This paper describes methodology and tools for model-agnostic audit. Introduced techniques facilitate assessing and comparing the goodness of fit and performance of models. In addition, they may be used for analysis of the similarity of residuals and for the identification of outliers and influential observations. The examination is carried out by diagnostic scores and visual verification. Presented methods are implemented in the auditor package for R. Due to the flexible and consistent grammar, it is simple to validate models of any classes.

In this paper, we introduce the auditor package for R, which is a tool for diagnostic and visual verification. As it focuses on residuals 1 and does not require any additional model assumptions, most of the presented methods are model-agnostic. A consistent grammar across various tools reduce the amount of effort needed to create informative plots and makes the validation more convenient and available.
Diagnostic methods have been a subject of much research (Atkinson, 1985). Atkinson and Riani (2012) focus on graphical methods of diagnostic regression analysis. Liu et al. (2017) present an overview of interactive visual model validation. One of the most popular tools for verification are measures of the differences between the values predicted by a model and the observed values (Willmott et al., 1985). This includes Root Mean Square Error (RMSE) and Mean Absolute Error (MAE) (Hastie et al., 2001). Such measures are used for well-researched and easily interpretable linear model as well as for complex models like a random forest (Ho, 1995), an XGBoost (Chen and Guestrin, 2016) or a neural network (Venables and Ripley, 2002).
Figure 1: Anscombe Quartet data sets are identical when examined using simple summary statistics. The difference is noticeable after plotting the data.
However, no matter which measure of model performance we use, it does not reflect all aspects of the model. As Breiman (2001) indicates, the linear regression model validated only on the basis of R 2 may lead to many false conclusions. The best known example for this issue is the Anscombe Quartet (Anscombe, 1973). It contains four different data sets constructed to have nearly identical simple statistical properties like mean, variance, correlation, etc. These measures directly correspond to the coefficients of the linear models. Therefore, by fitting a linear regression to the Anscombe Quartet we obtain four almost identical models (see Figure 1). However, residuals of these models are very different. The Anscombe Quartet is used to highlight that the numerical measures should be supplemented by graphical data visualizations. In Chapter 5 we show more advanced example related to this issue.
The diagnostic analysis is well researched for linear and generalized linear models. This is typically done by extracting raw, studentized, deviance, or Pearson residuals and examining residual plots. Common problems with model fit and basic diagnostic methods are presented in Faraway (2002) and Harrell (2006) (see Figure 4). Model validation may involve both, checking the overall trend in residuals and looking at residual values of individual observations (Littell et al., 2007). Gałecki and Burzykowski (2013) discussed methods based on residuals for individual observation and groups of observations. Diagnostic methods are commonly used for linear regression (Faraway, 2004). Complex models are treated as if they were assumption free, that is why diagnostic of them is often ignored. Considering the above, there is a need for more extensive methods and dedicated software for model auditing. Many of diagnostic tools like plots and statistics developed for linear models are still useful for exploring any machine learning model. Applying the same tools to all models makes it easy to compare them. Yet another approach to auditability of models was proposed in Becker and Chambers (1988).
The paper is organized as follows. Chapter 2 summarizes related work and state of the art. Chapter 3 contains an architecture of the auditor package. In Chapter 4 we provide the notation and discuss four main aspects of the model audit. In Section 4.1 we describe methods that help to assess the goodness-of-fit and may be useful to check whether the model does not miss relevant information. In Section 4.3 we present tools that help to assess which model has better performance. In Section 4.2 we describe graphical visualizations which allow checking the similarity of the models behavior. Section 4.4 contains methods for identification of influential observations. A use-case of model auditing is in Chapter 5 and conclusions are in Chapter 6.

Related work
In this chapter, we overview common methods and tools for auditing and examining the validity of the models. There are several attempts to validation. This includes diagnostic for predictor variables before and after model fit, methods dedicated to specific models, and model-agnostic approaches.

Data diagnostic before model fitting
The problem of data diagnostic is related to the Objective 2 presented in Introduction, that is, the identification of problems with observations. There are several tools that address this issue. We review the most popular of them.
• One of the tools that support the identification of errors in data is the dataMaid package (Petersen and Ekstrom, 2018). It creates a report that contains summaries and error checks for each variable in data. Package lumberjack (van der Loo, 2017) provides row-wise analysis. It allows for monitoring changes in data as they get processed. The validatetools (de Jonge and van der Loo, 2018) is a package for managing validation rules.
Figure 2: Page from an exemplary report generated with dataMaid package. Here we have a table that summarizes checks for variables. This includes, among others, identification of missing values and outliers.
In the further part of the report, there are summary tables and plots for each variable separately. The Figure  comes from the report generated with the code from the dataMaid documentation.
• The datadist function from rms package (Harrell Jr, 2018) computes distributional summaries for predictor variables. The characteristics include the overall range and certain quantiles for continuous variables, and distinct values for discrete variables. It automates the process of fitting and validating several models, due to storing model summaries by datadist function.
• While above packages use pipeline approaches, there are also tools that focus on a specific step of data diagnostic. The package corrgram (Wright, 2018) calculates a correlation of variables and displays corrgrams. Corrgrams (Friendly, 2002) are visualizations of correlation matrices, which help to identify the relationship between variables. The plot is taken from Friendly (2002) and generated with the code from vignette "Examples for the corrgram package".

Diagnostic methods for linear models
As linear models have a very simple structure and do not require high computational power, they have been and still are used very frequently. Therefore, there are many tools that validate different aspects of linear models. Below, we overview the most known of them that are implemented in R packages.
• The stats package (R Core Team, 2018) provides basic diagnostic plots for linear models. Function plot generates six types of charts for "lm" and "glm" objects. A plot of residuals against fitted values, a scale-location plot of |residuals| against fitted values and a normal quantile-quantile plot. These visual validation tools may be referred to the Objective 3 of diagnostic, related to the examination of model by analyzing the distribution of residuals. The other three plots, that are: a plot of Cook's distances, a plot of residuals against leverages, and a plot of Cook's distances against leverage 1−leverage may be referred to the identification of influential observations (Objective 1). Cook's distance and leverages are introduced in details in Section 27. Figure 4: Six diagnostic plots for linear models created with generic plot function. Upper row contains plots related to the distribution of residuals, lower row is related to an influence of observations. Plots are generated with the code included in the documentation of function plot.lm.
• Package car (Fox and Weisberg, 2011) extends the capabilities of stats by including more types of residuals, such as Pearson and deviance residuals. It is possible to plot against values of selected variables and grouping residuals by levels of factor variables. What is more, car provides more diagnostic plots. Among others, partial residual plot (crPlot), index plots of influence (infIndexPlot) and bubble plot of studentized residuals versus hat values (influencePlot). These plots allow to check both, the effect of observation and the distribution of residuals, which relates to the Objective 2 and the Objective 3 respectively.
Figure 5: Two diagnostic plots generated with the car package. Left plot is a partial residual plot, right is a bubble plot of studentized residuals versus hat values. The areas of circles represent Cook's distances. Plots were generated with code from the documentation of car package.
• A linear regression model is still one of the most popular tools for data analysis due to its simple structure. Therefore, there is a rich variety of methods for checking its assumptions. For example, the normality of residual distribution and the homoscedasticity of the variance. The package nortest (Gross and Ligges, 2015) provides five tests for normality: the Anderson-Darling (Anderson and Darling, 1952), the Cramer-von Mises (Cramer, 1928;Von Mises, 1928), the Kolmogorov-Smirnov (Stephens, 1974), the Pearson chi-square (F.R.S., 1900), and the Shapiro-Francia (Sanford Shapiro and S. Francia, 1972) tests. The lmtest package (Zeileis and Hothorn, 2002) also contains a collection of diagnostic tests: the Breusch-Pagan (Breusch and Pagan, 1979), the Goldfield-Quandt (Goldfeld and Quandt, 1965) and the Harrison-McCabe (Harrison and McCabe, 1979) tests for heteroscedasticity and the Harvey-Collier (Harvey and Collier, 1977), the Rainbow (Utts, 1982), and the RESET (Ramsey, 1969) tests for nonlinearity and misspecified functional form. A unified approach for examining, monitoring and dating structural changes in linear regression models is provided in strucchange package . This includes methods to fit, plot and test fluctuation processes and F-statistics. The gvlma implements the global procedure for testing the assumptions of the linear model, find more details in Peña and Slate (2006). The Box-Cox power transformation, introduced by Box and Cox (1964), is a way to transform the data to follow a normal distribution. For simple linear regression, it is often used to satisfy the assumptions of the model. Package MASS (Venables and Ripley, 2002) contains functions that compute and plot profile log-likelihoods for the parameter of the Box-Cox power transformation. • The broom package (Robinson, 2018) provides summaries for about 30 classes of models. It produces results, such as, coefficients and p-values for each variable, R 2 , adjusted R 2 , and residual standard error.

Other model-specific approaches
There are also several tools to generate validation plots for time series, principal component analysis, clustering, and others.
• Tang et al. (2016) introduced the ggfortify interface for visualizing many popular statistical results. Plots are generated with ggplot2 (Wickham, 2009), which makes them easy to modify. With one function autoplot it is possible to generate validation plots for wide range of models. It works for, among others, lm, glm, ts, glmnet, and survfit objects. The autoplotly (Tang, 2018) package is an extension of ggfortify and provides functionalities that produce plots generated by plotly (Sievert et al., 2017). This allows for both, modification and interaction with plots. However, ggorftify and autoplotly do not support some popular types of models. For instance, random forests from randomForest (Liaw and Wiener, 2002) and ranger (Wright and Ziegler, 2017) packages. • The hnp package (Moral et al., 2017) provides half-normal plots with simulated envelopes. These charts evaluate the goodness of fit of any generalized linear model and its extensions. It is a graphical method for comparing two probability distributions by plotting their quantiles against each other. The more detailed description is provided in the Section 4.1.7. There is a possibility to extend the hnp for new model classes. However, this package provides only one tool for model diagnostic. In addition, plots are not based on ggplot2, which makes it difficult to modify them.

Model-agnostic approach
Tools, presented above, target specific model classes. The model-agnostic approach allows to compare different models.
• The DALEX (Descriptive mAchine Learning EXplanations) (Biecek, 2018) is a methodology for exploration of black-box models. Main functionalities focus on understanding or proving how the input variables impact on final predictions. There are also two simple diagnostic diagnostic: reversed empirical cumulative distribution function for absolute values of residuals and box plot of absolute values of residuals. As methods in the DALEX are model-agnostic, they allow to compare two or more models. Figure 8: Model explanation plots generated with the DALEX package. The upper plot is a Partial Dependence Plot (Greenwell, 2017). The lower plot is an empirical cumulative distribution function of absolute values of residuals. Plots were taken from Biecek (2018).

Model-agnostic audit
In this paper, we present the auditor package for R, which fills out the part of model-agnostic validation.
As it expands methods used for linear regression, it may be used to verify any predictive model.

Package Architecture
The auditor package works for any predictive model which returns a numeric value. It offers a consistent grammar of model validation, which is an efficient and convenient way to generate plots and diagnostic scores. A diagnostic score is a number that evaluates one of the properties of a model. That might be, for example, an accuracy of model, an independence of residuals or an influence of observation. Figure 9 presents the architecture of the package. The auditor provides 2 pipelines for model validation.
First of them consists of two steps. Function audit wraps up the model with meta-data, then the result is passed to the plot or score function. Second pipeline includes an additional step, which consists of calling one of the functions that generate computations for plots and scores. These functions are: modelResiduals, modelEvaluation, modelFit, modelPerformance, and observationInfluence. Further, we call them computational functions. Results of these functions are tidy data frames (Wickham, 2014). Both pipelines for model audit are compared below.

model %>% audit() %>% computational function %>% plot(type=. . . )
We recommend this pipeline. Function audit wraps up a model with meta-data used for modeling and creates a "modelAudit" object. One of the computational functions takes "modelAudit" object and computes the results of validation. Then, output may be printed or passed to functions score and plot with defined type. We describe types of plots in the Chapter 4.
This approach requires one additional function within the pipeline. However, once created output of the computational function contains all necessary calculations for related plots. Therefore, generating multiple plots is fast.

model %>% audit() %>% plot(type=. . . )
This pipeline is shorter than the previous one. The only difference is that it does not include computational function. Calculations are carried out every time a generic plot function is called. Omitting one step might be convenient for ad-hoc model analyses.
Implemented types of plots are presented in Table 1. Scores are presented in Table 2. All plots are generated with ggplot2. This provides a convenient way to modify and combine plots.  Table 2: Columns contain as follows: name of a score, name of a computational function, value for type parameter of function the score, indications whether the score can be applied to regression and classification tasks.

Model audit
Diagnostic allows to evaluate different properties of a model. We divide them into four groups. Each of them is related to one of the following aspects and corresponding questions.
1. First aspect is assessing the goodness-of-fit and whether the model does not miss relevant information (Does the model fit data?). See subsection 4.1.
2. Second one is examining similarity of models (How similar models are?). See subsection 4.2.
3. Third aspect is an evaluation of model performance. ( These aspects are directly related to the diagnostic objectives described in the Introduction. First two of them are related to the examination of distribution of residuals, which was proposed as the Objective 3. The third aspect is an evaluation of a model performance (Objective 1). The last one refers to influential observations (Objective 2).
We introduce notation to follow throughout the paper.
Let us use the following notation: A single observation we denote as a pair (y i , x i ) and n is the number of observations. Let denote a model as a function f : X → R. Predictions of the model f for particular observation we denote as Row residual or, simply, residual is the difference between the observed value y i and the predicted valueŷ i . We denote residual of particular observation as

Aspect: Does the model fit data?
A good fitted model predicts response variable well. In this subsection, we present methods to assess the goodness-of-fit and check whether the model does not miss any relevant information.
To illustrate applications of the auditor we use an artificial data set apartments available in the DALEX package. First, we fit two models: simple linear regression and random forest.

Residuals Plot
Residuals vs Fitted Values is the most frequently used plot in model validation. It is a scatter plot of residuals r i on the y axis against fitted valuesŷ i on the x axis. Example plot is presented in Figure 10. On alterations of this plot on the x-axis are values of the variable.
This plot is used to detect dependence of errors, unequal error variances, and outliers. For appropriate model, residuals should not show any functional dependency. Expected mean value should be equal 0, regardless ofŷ values. Structured arrangement of points suggests a problem with the model. It is worth looking at the observations that clearly differ from the others. If points on the plot are not randomly dispersed around the horizontal axis, it may be presumed that model is not appropriate for the data.
This plot is generated by the plot function with parameter type = "Residual" or by plotResidual function.
Other variants of the Residual Plot may be obtained by parameter variable.

Residual Boxplot Plot
Residual boxplot shows the distribution of the absolute values of residuals r i . They may be used for analysis and comparison of residuals. Example plots are presented in Figure 11.
Boxplots (Tukey, 1977) usually consists of five components. Box corresponds to the first quartile, median, and third quartile. The whiskers extend to the smallest and largest values, no further than 1.5 of Interquartile plot(lm_res_fitted, rf_res_fitted, type = "ResidualBoxplot") Figure 11: Boxplots of absolute values of residuals. Dots are in the similar places, so RMSE for both models is almost identical. However, the distribution of residuals of this two models is different. For the linear model (blue), most of the residuals are around the average. For the random forest (red), most residuals are small. Nevertheless, there is also a fraction of large residuals.

Autocorrelation Plot
Autocorrelation Plot is a tool for checking whether there is a relationship between residuals on the lag 1. Example chart is presented in Figure 12. It is a scatter plot of i-th vs i + 1-th residuals, ordered by fitted values or by values of one of the variables.
We expect that model residuals are independent. Therefore, points should be randomly dispersed. The structured arrangement of residuals suggests a problem with the model. This plot is generated by plot function with parameter type = "Autocorrelation" or by function plotAutocorrelation.
plot(lm_res_fitted, rf_res_fitted, type = "Autocorrelation") Figure 12: Autocorrelation plot of i-th residual vs i + 1-th residual. The residuals of the linear model (blue) are divided into 4 groups, which indicates that they are strongly non-random. Residuals of random forest (red) show an increasing trend, which also suggests a problem with the model structure.
As results of plots may be ambiguous, there are two scores implemented in the auditor: Durbin-Watson (DW) and Runs. They help to asses the autocorrelation of residuals.
The Durbin-Watson (DW) statistic is used to detect the autocorrelation in the residuals at lag = 1. It is introduced by Durbin and Watson (1971). Definition 4.1. The Durbin-Watson test statistic is given by It detects only linear autocorrelation and only first order effects.
We treat Durbin-Watson statistic as a score. Its values lie between 0 and 4. DW around 2 means no autocorrelation. Small value may be an evidence for positive correlation and large value for negative correlation.
Runs score is the second tool for checking the independency of residuals.
A run is a series of increasing or decreasing values. The U is an observed number of runs in residuals ordered by fitted values or values of one of the variables. TheŪ is an expected number of runs and s U is the standard deviation of the number of runs. Runs test statistic was defined by Wald and Wolfowitz (1940). Definition 4.2. The Runs test statistic is given by For random sequence, the distribution of Z follows the N (µ, σ 2 ).
We consider values of Runs statistic as scores.
Both presented scores may be obtained by score function with type = "DW" or type = "Runs".

Autocorrelation Function Plot
Autocorrelation Function (ACF) is a tool for finding patterns in the data. In particular, the correlation between observations. Figure 13 shows an example of the Autocorrelation Function Plot of residuals.
There are values of ACF on the y axis against lags on the x axis. As the value of autocorrelation of lag 0 always equals 1, it is skipped. Blue horizontal lines are confidence intervals.
We assume that (X 1 , X 2 , ...X n ) ∈ R n are observations of a time series X t with mean µ. Covariance function γ t = Cov(X t+τ , X τ ) and correlation function ρ t = Cov(X t+τ , X τ ) do not depend on τ . The following estimators come from Venables and Ripley (2013). Definition 4.3. The sample autocovariance function is given bŷ and the sample autocorrelation function is given bŷ Now, let us consider residuals r i as elements of the time series X t . Autocorrelation Function Plot visualizes the sequence of (ρ(t)).
This plot is generated by plot function with parameter type = "ACF" or by function plotACF.

Scale-Location Plot
Scale Location Plot is used to visualize the variance of the residuals. Example chart is presented in Figure 14. This plot is similar to the Residuals Plot, but it uses the square root of the standardized residuals instead of raw residuals. Letr be mean of residuals,r Definition 4.4. The standardized residual r std i is the i-th residual divided by the residual standard deviation: On the x axis, there are responses fitted by a model. On the y axis, there are square roots of the absolute value of the standardized residuals. Blue line is an estimate of the conditional mean function. Black dots represent peaks explained below. An alternative version of this plot contains values of one of the variables on the x axis.
The presence of any trend suggests that the variance depends on fitted values, which is against the assumption of homoscedasticity. Not every model has an explicit assumption of homogeneous variance, however, the heteroscedasticity may indicates potential problems with the goodness-of-fit. This plot is generated by plot function with parameter type = "ScaleLocation" or by function plotScaleLocation. Variable on the x axis may be specified by variable parameter.
plot(rf_res_fitted, lm_res_fitted, type = "ScaleLocation") Figure 14: Scale Location plot. The trend in residuals of random forest (red) suggests that variance of residuals changes with the change in fitted values. The linear model (blue) seems to have no problems with variance. However, residues formed into separate groups suggest a problem with model structure.
A tool for assesing the homoscedasticity of variance is Peak test introduced by Goldfeld and Quandt (1965). A peak is an occurrence of observation i such that |r i | ≥ |r j | for all i > j.
Let us note that P ∈ (0, 1]. If the residuals r i are heteroscedastic, the Peak score is close to 1.
Scores may be obtained by score function with type = "Peak" or scorePeak function.

Residual Density Plot
Residual Density plot detects the incorrect behavior of residuals. Example is presented in Figure 15. On the plot, there are estimated densities of residuals. Their values are displayed as marks along the x axis.
For some models, the expected shape of density derives from the model assumptions. For example, simple linear model residuals should be normally distributed. However, even if the model does not have an assumption about the distribution of residuals, such a plot may be informative. If most of the residuals are not concentrated around zero, it is likely that the model predictions are biased.
This plot is generated by plotResidualDensity function or by plot function with parameter type = "ResidualDensity".
plot(rf_res_observed, lm_res_observed, type = "ResidualDensity") Figure 15: Residual Density Plot. The density of residuals for the linear model (blue) forms two peaks. There are no residuals with values around zero. Residuals do not follow the normal distribution, which is one of the assumptions of the simple linear regression. There is an asymmetry of residuals generated by random forest (red).

Half-Normal Plot
The Half-Normal Plot is one of the tools designed to evaluate the goodness-of-fit of a statistical model. It is a graphical method for comparing two probability distributions by plotting their quantiles against each other.
Points correspond to ordered absolute values of model diagnostic (i.e. standardized residuals) plotted against theoretical order statistics from a half-normal distribution. Lines mark the simulated envelope. The simulation process consists of simulating response variables using the model matrix, error distribution and parameters of an original model. Then the same model is fitted to data with new simulated response variable. The envelopes are formed on the basis of the 2.5 and 97.5 percentiles of the extracted simulated residuals at each value of the expected order statistic. Further details on the half-normal plots with simulated envelopes may be found in Moral et al. (2017).
If residuals come from the normal distribution, they are close to a straight line. However, even if there is no certain assumption of a specific distribution, points still show a certain trend. Simulated envelopes help to verify the correctness of this trend. For a good-fitted model, diagnostic values should lay within the envelope.
This plot is generated by plot function with parameter type = "HalfNormal" or by function plotHalfNormal. The parameter quant.scale = TRUE transforms y axis into quantile scale.
lm_mf <-modelFit(lm_audit) plot(lm_audit, type = "HalfNormal") Figure 16: Half-Normal plot with simulated envelope. Most of the residuals lie outside the envelope. They do not behave as we would expect from the simulations. That suggests a poor model fit.
A useful tool to compare goodness-of-fit of two models is a Half-Normal Score.
Let us consider n observations and m simulations for each observation. We use the following notation: r i is a residual of i-th observation, r simj i is the residual of j-th simulation for i-th observation. Definition 4.6. The score of i-th observation is given by S i is a number of simulated residuals for observation i that are greater or equal than r i .
Value of S i close to 0 or m means that the original residual stands out from the simulated ones. The closer it is to m 2 , the less it deviates from the simulated results. Definition 4.7. We define Half-Normal score as HN is the sum of the deviations of S i from m 2 .
Let us note that HN ∈ [0, nm 2 ] and lower value of score means better model fit. HalfNormal score is calculated by score function with parameter type = "HalfNormal" or by function soreHalfNormal. Scores are calculated on the basis of simulated data, so they may differ between function calls.

Aspect: How similar models are?
In this subsection, we present methods for assessing the similarity of models, understood as the closeness of the residuals structures. The following methods may be used to identify such a similarity or its absence.

Model PCA Plot
Model PCA plot can be used to assess the similarity of the models in terms of residuals. Example plot is presented in Figure 17.
The idea of PCA is reducing the dimension of the data set matrix by creating a set of linearly uncorrelated variables called principal components. At the same time, keeping as much variance as possible Jolliffe (1986).
Model PCA plot is a biplot introduced by Gabriel (1971). On axis of the plot, there are first two principal components. Grey dots represent observations. Arrows are pointing in the direction of the models projected into a two-dimensional space. The interposition of arrows provides information about the similarity of models in terms of residuals. If they are close to each other, it indicates similar residuals structures. This plot is generated by plot function with parameter type = "ModelPCA" or by function plotPCA. rf_mp <-modelPerformance(rf_audit) lm_mp <-modelPerformance(lm_audit) plot(rf_mp, lm_mp, type = "ModelPCA") Figure 17: Model PCA plot. The diagram shows the arrangement of the residuals in two groups. The arrow corresponding to the linear model is almost perpendicular to the groups direction. Therefore, such a division can also be seen in the structure of the residuals of this model. The arrow corresponding to the random forest is almost parallel to the direction of the residuals. Probably the division into groups does not appear in the residuals of this model.

Model Correlation Plot
Model Correlation plot, presented in Figure 18, is a pairs plot of observed response and fitted values of different models. Pairs plot first appeared in Hartigan (1975). In the diagonal, there are estimated densities, in the upper triangle, the correlations between the different models and between models and observed response. In the lower triangle, there are scatter plots of these values. This plot is generated by plot function with parameter type = "ModelCorrelation" or by function plotCorrelation. Taking the parameter values = "Residuals" residuals will be plotted instead of fitted values.

Aspect: Which model has better performance?
We consider the model performance more widely than the value of one measure. The values and structures of the residuals also help to assess the quality of the model. In this subsection, we present visual methods of evaluating models performance.

Predicted Response Plot
Predicted Response Plot gives an additional information about model performance. It is a scatter plot of predicted valuesŷ i on the y axis against observed response y i on the x axis. Example plot is presented in Figure 19. On alterations of this plot on the x-axis are values of one the model variables. Ideally, all points are close to a diagonal line. For the appropriate model, points should not show any functional dependency. This plot is generated by plot function with parameter type = "Prediction" or by plotPrediction function. Other variable on the x axis may be specified by variable parameter. plot(rf_res_observed, lm_res_observed, type = "Prediction") Figure 19: redicted Response Plot. The patterns for both models are non-random around the diagonal line. The points corresponding to a random forest (red) show the tendency to underprediction for large values of observed response and over-prediction for small values. Points for linear model (blue) are arranged into two separate groups. That suggests problems with the structures of the models.

Plots: ROC Curve and LIFT Chart
In this subsection, we present methods for assessing model performance for classification problem. Most common methods for regression models were inspired by tools for classification.
Let us consider a binary classification problem, in which the responses are labeled either as positive (P ) or negative (N ). True Positive (T P ) is a case that the prediction is P and the actual value is also P . If the prediction is P and the actual value is N then it is a False Positive (F P ).
First, we fit two binary classification models on the PimaIndianDiabetes data set from the mlbench package.

Receiver Operating Characteristic (ROC) curve (See
Let t ∈ [0, 1] be a threshold. We introduce parametric definition of ROC curve. Definition 4.8. The ROC curve is a plot of the True Positive Rate (T P R) against the False Positive Rate (F P R) on a threshold t. and Each point on the ROC curve represents values of T P R and F P R of different thresholds.
The diagonal line y = x corresponds to a classifier that randomly guess the positive class half the time. Any model that appears in the lower right part of plot performs worse than random guessing. The closer the curve is to the the left border and top border of plot, the more accurate the classifier is.
This plot is generated by plot function with parameter type = "ROC" or by function plotROC.
plot(au_glm_class, au_svm_class, type = "ROC") Figure 20: Receiver Operating Characteristic curve. The svm (red) curve is mostly above glm (blue). In the remaining part they overlap. This indicates that the support vector machines performs better than linear model.
ROC curves usually intersect. Therefore, it is not possible to choose a better model explicitly. A good tool to help assess the quality of a classifier is the Area Under the Curve (AUC) (Bradley, 1997). The AUC is a part of the area of the unit square. Therefore, its value is always be between 0 and 1. The AUC of a classifier is equal to the probability that the classifier will rank a randomly chosen positive instance higher than a randomly chosen negative one (Fawcett, 2006). AUC may be obtained by score function with type = "ROC" or scoreROC function.
score(au_glm_class, type = "ROC") score(au_svm_class, type = "ROC") LIFT charts (Witten et al., 2011) also evaluate performance of classification models. Vuk and Curk (2006) introduced parametric definition of LIFT chart. Let t ∈ [0, 1] be a threshold. Definition 4.9. The LIFT chart is Rate of Positive Prediction (RP P ) plotted against True Positive (TP) on a threshold t. y = T P (t) (14) and Each point on the LIFT chart represents values of T P and RP P of different thresholds.
As for ROC curve, LIFT (See Figure 21) illustrates varying the model performance for different thresholds. LIFT chart shows how good the classifier distinguishes between two classes.
A random and ideal models are represented by black and orange curves, respectively. The closer the LIFT get to the orange curve, the better a model is.
This plot is generated by plot function with parameter type = "LIFT" or by function plotLIFT.
plot(au_glm_class, au_svm_class, type = "LIFT") Figure 21: LIFT Chart. Red model (svm) The red model (glm) is above the blue one, which suggests that it has a better performance.

REC Curve Plot
Regression Error Characteristic (REC) curve (See Figure 22) is a generalization of Receiver Operating Characteristic (ROC) curve for binary classification described in subsection 4.3.2.
REC curve estimates the Cumulative Distribution Function of the error. On the x axis of the plot there is an error tolerance. On the y axis there is an accuracy at the given tolerance level. Bi and Bennett (2003) define the accuracy at tolerance as a percentage of observations predicted within the tolerance . In other words, residuals larger than are considered as errors.
Let consider pairs (y i , x i ) as at the beginning of the Chapter 4. Bi and Bennett (2003) define an accuracy as follows.
Definition 4.10. An accuracy at tolerance level is given by REC Curves implemented in the auditor are plotted for a special case of Definition 4.10 where the loss is defined as The shape of the curve illustrates the behavior of errors. The quality of the model can be evaluated and compared for different tolerance levels. The stable growth of the accuracy does not indicate any problems with the model. A small increase of accuracy near 0 and areas, where the growth is fast, signalize bias of the model predictions.
This plot is generated by plot function with parameter type = "REC" or by plotREC function.
plot(rf_res_fitted, lm_res_fitted, type = "REC") Figure 22: REC curve. Curve for linear model (blue) suggests that model is biased. It has poor accuracy when the tolerance is small. However, once exceeds the error tolerance 130 the accuracy rapidly increases. The random forest (red) has a stable increase of accuracy when compared to the linear model. However, there is s fraction of large residuals.
As often it is difficult to compare models on the plot, there is an REC score implemented in the auditor. This score is the Area Over the REC Curve (AOC), which is a biased estimate of the expected error for a regression model. As Bi and Bennett (2003) proved, AOC provides a measure of the overall performance of regression model. Scores may be obtained by score function with type = "REC" or scoreREC function.

Two-sided ECDF Plot
Two-sided ECDF plot (See Figure 23) shows an Empirical Cumulative Distribution Functions (ECDF) for positive and negative values of residuals separately.
Let x 1 , ..., x n be a random sample from a cumulative distribution function F (t). The following definition comes from van der Vaart (2000). Definition 4.11. The empirical cumulative distribution function is given by Empirical cumulative distribution function presents a fraction of observations that are less or equal than t.
It is an estimator for the cumulative distribution function F (t).
On the positive side of the x-axis, there is the ECDF of positive values of residuals. On the negative side, there is a transformation of ECDF: Let n N and n P be numbers of negative and positive values of residuals, respectively. Negative part of the plot is normalized by multiplying it by the ratio of the n N over n N + n P . Similarly, positive part is normalized by multiplying it by the ratio of the n P over n N + n P . Due to the scaling, the ends of the curves add up to 100% in total.
The plot shows the distribution of residuals divided into groups with positive and negative values. It helps to identify the asymmetry of the residuals. Points represent residuals. This plot is generated by plotTwoSidedECDF function or by plot function with parameter type = "TwoSidedECDF". plot(rf_res_fitted, lm_res_fitted, type = "TwoSidedECDF") Figure 23: Two-sided ECDF plot. The plot shows that majority of residuals for the random forest (red) is smaller than residuals for the linear model (blue). However, random forest has also fractions of large residuals.

RROC Curve Plot
The Regression Receiver Operating Characteristic (RROC) curve for regression shows model asymmetry, which is an inequality in the number and values of positive and negative residuals. It may be useful for problems, where there is an asymmetric cost of errors. An example of RROC plot is presented in Figure 24.
The RROC curve was introduced by Hernández-Orallo (2013). The base of plot is a shift s. It is an equivalent to the threshold for ROC curves presented in subsection 4.3.2. For each observation we calculate new prediction:ŷ s =ŷ + s. Hernández-Orallo (2013) define over-and under-estimation as follow. Definition 4.12. Over-estimation depending on shift s is given by Definition 4.13. Under-estimation depending on shift s is given by The RROC plot consists of under-estimation on the y axis against over-estimation on the x axis. The shift equals 0 is represented by a dot.
The shape of the curve illustrates the behavior of errors. The quality of the model can be evaluated and compared for different shifts. The AOC measures the model insensitivity to asymmetric costs of errors.
This plot is generated by plot function with parameter type = "RROC" or by plotRROC function.
plot(rf_res_fitted, lm_res_fitted, type = "RROC") As often it is difficult to compare AOC on the plot, there is an RROC score implemented in the auditor package. This score is the Area Over the RROC Curve (AOC). As Hernández-Orallo (2013) prove AOC equals to the variance of the errors multiplied by a n 2 2 which is independent of the model. Lower values for AOC are better. Score may be obtained by score function with type = 'RROC' or scoreRROC function.

Model Ranking Plot
In this subsection, we propose a Model Ranking plot, which compares models performance across multiple measures (see Figure 25). The implemented measures coincide with the scores listed in the Chapter 3.
In previous subsections, we introduced the following scores: Additional scores implemented in the auditor are: • MAE (Mean Absolute Error) • MSE (Mean Squared Error) • RMSE (Root Mean Squared Error) Model Ranking Radar plot consists of two parts. On the left side, there is a radar plot. Colors correspond to models, edges to values of scores. Score values are inverted and rescaled to [0, 1].
Let us use the following notation: m i ∈ M is a model in a finite set of models M, where |M| = k, score : M → R is a scoring function for the model under consideration, that higher value means worse model performance. The score(m i ) is a performance of model m i . Definition 4.14. We define the inverted score of model m i as Models with the larger invscore are closer to the centre. Therefore, the best model is located the farthest from the center of the plot.
On the right side of the plot there is a table with results of scoring. In the third column, there are scores scaled to one of the models.
Let m l ∈ M where l ∈ {1, 2, ..., k} be a model to which we scale. Definition 4.15. We define the scaled score of model m i to model m l as As values of scaled l (m l ) are always between 0 and 1, comparison of models is easy, regardless of the ranges of scores.
This plot is generated by plot function with parameter type = "ModelRanking" or by function plotModelRanking. The scores included in the plot may be specified by scores parameter.
rf_mp <-modelPerformance(rf_audit) lm_mp <-modelPerformance(lm_audit) plot(rf_mp, lm_mp, type = "ModelRanking") There is also a possibility to add custom scores ( Figure 26). They may be provided by parameter new.score. It requires a named list of functions that take one argument: object of class "modelAudit" and return a numeric value. The measure calculated by the function should have the property that lower score value indicates better model.

Aspect: Influence of Observations
In this subsection, we focus on the impact of individual observation on a model.

Cook's Distances Plot
Cook's distances plot presented in Figure 27 is a tool for identifying observations that may negatively affect the model. They can be also used for indicating regions of the design space where it would be good to obtain more observations. Data points indicated by Cook's distances are worth checking for validity. Let us extend the notation provided iat the beginning of the Chapter 4. We denote the number of predictors as p and the mean squared error as s 2 . We consider recalculated model that is fitted on the original data set with removed j-th observation. Letŷ i(j) be the prediction of such model calculated for i-th observation. Cook (1977) defined the influence of a single observation as follows. Definition 4.16. Cook's distance of i-th observation is Cook's Distances are calculated by removing the i-th observation from the data and recalculating the model. It shows an influence of i-th observation on the model. Cook (1977) proved that for linear models Cook's distance may be computed in an alternative and computationally convenient way. Let X be the design matrix and H = X(X X) −1 X be a projection matrix. The leverage h i is the i-th diagonal element of H. D i is obtained by This plot is generated by plot function with parameter type = "CooksDistance" or by function plotCooksDistance. For model of class "lm" and "glm" the distances are computed from the diagonal elements of the hat matrix. For other models they are computed directly from the definition.

Example of model audit
This section contains a use-case of model audit. We show an example, where the choice of a model on the basis of a single measurement is not obvious. We illustrate, how to use auditor to analyze and compare models. In addition, we present how to use residual analysis to identify outliers and better understand a model structure.

The artificial data set
In order to present the auditor we created an artificial data set auditorData. It consists of 2000 observations. We generated first 1998 as follows. Let 2 , ..., x 2 , ..., x for i ∈ 1, 2, ..., 1998.
Let us simulate a response Y = (y 1 , y 2 , ..., y 1998 ) as a function of five arguments: X 1 , X 2 , X 3 , X 4 , , where The 1999-th and 2000-th observations are added manually. They are meant to be outliers. Their values are 1999 , x 1999 , x 1999 , x 2000 , x 2000 , x 2000 , x First four of simulated variables are treated as continuous while the fifth one is categorical. In the Table 3 we present first six observations from the auditorData included in the auditor package.

Model audit
At first, we generate four diagnostic plots. They are: Model Ranking Plot, Predicted Response Plot, PCA of Models, and Residuals Plot. We present results in Figure 28 and Figure 29.
plot(au_rf, au_svm, au_lm, type = c("ModelRanking", "Prediction"), variable = "Observed response", smooth = TRUE) On the Model Ranking plot (see Figure 28) we can see that random forest has the best performance in terms of M SE and RROC sores. Support vector regression is the best in terms of M AE and REC scores. Linear regression appears to have the worst performance in every aspect. As none of the models performs best in terms of all measures, it is clear that Model Ranking plot is not enough to evaluate models. There is a need for further analysis of residuals.
Prediction Response plot (see Figure 28) indicates that there are observations that may be outliers and have significant influence on model's structures. Most of the residuals of all models are arranged along a black line that shows the ideal trend. While, there is a group of points that clearly stands out from the rest. Further plots allow to take a closer look at these observations.
PCA of Models (see Figure 29) confirms that there are 2 observations that have high influence on the structures of all models. Residuals plot with labeled points shows the numbers of outliers. These are 1999-th and 2000-th observations. The same that we added artificially.

Removing identified outliers and models improvement
In previous section, we identified two outliers in auditorData. In this section, we remove those outliers, fit models to a new data set, and create "modelAudit" objects.
plot(au_rf, au_svm, au_lm, type = c("ModelRanking", "Prediction"), variable = "Observed response", split = "model") This time, the Model Ranking plot clearly indicates that support vector regression performs best. However, to better exploration of models, it is still worth to carry out a further audit. The Predicted Response Plot shows that residuals of support vector regression are closest to the diagonal line than residuals of other models. Residuals of linear model have the largest dispersion, while random forest residuals indicate over-prediction for small values of observed response and under-prediction for large values.

Extended model audit
In this section we take a closer look on the residuals. We generate Residual Boxplot and Two-sided ECDF plot for better comparison of models in aspect of residuals.
In Figure 30, we noticed that support vector regression has the smallest residuals. By analyzing boxplots (see Figure 31), we can additionally see that random forest has smaller residuals than linear model. Two-sided ECDF plots show that conclusions from boxplots are, in general, correct. In addition, we can see that growth of small negatives residuals for random forest and linear model is similar. Now, we audit residuals due to the model variables. We use Residual Density to understand residual behaviour for different values of the X4 variable.
plot(au_rf, au_svm, au_lm, type = "ResidualDensity", variable = "X4") In the Figure 32, we see that structure of residuals for random forest and support vector regression do not vary due to the value of X4 variable. In contrast, there are differences in the shape of the densities for linear model's residuals. In the data set, we included the X4 variable in the interaction with X1. Due to its structure, the linear model do not catch this interaction. We use Residuals Plot to examine behaviour of residuals due to second variable included in the interaction.
In the Figure 33, we see that residuals form three groups. Based on the previous analysis, we conclude that groups correspond to the levels of the X4 variable.

Summary
We fitted three different models to the data set auditData. We used Predicted Response, Model PCA and Residuals plots to identify two observations, that highly influenced model's structures. Then we refitted models on data without outliers. The audit showed that support vector regression has the best performance and we did not identify any serious problems with the structure of the residuals.

Conclusion and future work
In this article, we presented a number of diagnostic scores and plots. We discussed the existing methods of model validation and proposed their adaptations to any model. We also specified three objectives of model audit, proposed relevant verification tools, and demonstrated their usage by incorporating examples.
• First objective of model diagnostic was the examination of a model performance. It involves using scores and supplementing them by visual validation. In Section 4.3 we showed tools that may be used for that purpose. These include scores such as Mean Squared Error, Area Over the RROC, and Area Over the REC. In addition, we presented plots that enrich information about model performance. They are, among others, Model Ranking Plot, Receiver Operating Characteristic Curve, and Regression Error Curve. • Second objective was the identification of observations that are worth further examination. The detection of influential observations is possible with Cook's Distances Score and related Cook's Distances Plot presented in Section 4.4. • Third considered objective was the recognition of problems that might be identified by analysis of distribution of residuals. With tools described in Sections 4.1 and 4.2, it is possible to answer questions like: "Does the model fit data?"," How similar models are?", and "Does the model miss any relevant information?". Examination of the behavior of residuals provides valuable information about a model. The scores related to this issue were, for example, Runs and Durbin-Watson scores. Presented plots were, among others, Residuals Plot, Residual Boxplot, and Residual Density.
Aside from describing existing methods and extending them to any model, we proposed new plots and scores. They are Model Ranking Plot, Two-Sided ECDF Plot and Half-Normal Score.
We implemented all introduced scores and plots in the auditor package for R. Included functions are based on a uniform grammar introduced in Figure 3. Documentation and examples are available at https: //mi2datalab.github.io/auditor/. The stable version of the package is on CRAN, the development version is on GitHub (https://github.com/MI2DataLab/auditor).
In Section 5 we showed the use-case of the model audit with the auditor. We have presented a broad exploration and comparisons of three models by analyzing their residuals.
There are many potential areas for future work that we would like to explore. This includes more extensions of model-specific diagnostic to model-agnostic methods and residual-based methods for investigating interactions. Another aim is to develop methods for local audit based on the diagnostic of a model around a single observation or a group of observations.