Categorical data are prevalent in all areas, including economics, marketing, finance, psychology, and clinical studies. To analyze categorical data, the probit or logit models are often used to make inferences. To perform model assessment and comparison, researchers often rely on goodness-of-fit measures, such as R 2 (also known as the coefficient of determination). For example, the ordinary least square (OLS) R 2 is one of the most extensively used goodness-of-fit measures for linear models in continuous data analysis. For categorical data analysis, however, there is no such universally adopted R 2 measure [19, 39]. Continuous efforts have been made to develop sensible R 2 measures for probit/logistic models, and more generally, generalized linear models [29, 30, 14, 11, 22, 42, 27, 21, 26]. Among the existing R 2 measures, McKelvey-Zavoina’s R M Z 2 [30] and McFadden’s R 2 [29] are probably the most well-known and widely used in domain research [19, 39]. But as demonstrated in [26], Mckelvey-Zavoina’s R M Z 2 does not hold monotonicity, which means a larger model may have a smaller R M Z 2 . This serious defection of R M Z 2 may be misleading in practice and misguide the model-building process. On the other hand, McFadden’s R 2 relies on the ratio of likelihoods, and it does not preserve the interpretation of explained variance. Neither of these two R 2 measures meets all of the three criteria considered in [26]:

[26] proposed a so-called surrogate R 2 that satisfies all three criteria for probit models. This surrogate R 2 uses the notion of surrogacy, namely, generating a continuous response S and using it as a surrogate for the original categorical response Y [24, 25, 8, 18, 23]. In the context of probit analysis, [26] used the truncated distributions induced by the latent variable structure to generate a surrogate response S. This surrogate response S is then regressed on explanatory variables through a linear model. The OLS R 2 of this linear model is used as a surrogate R 2 for the original probit model. This surrogate R 2 meets all three criteria (C1)–(C3).

The goals of this paper are (i) developing an R package to implement [26]’s method; (ii) demonstrating how this new package can be used jointly with other existing R packages for variable selection and model diagnostics in the model building process; and (iii) illustrating how this package can be used to compare different empirical models trained from two different samples (a.k.a. comparability) in real data analysis.

Specifically, we first develop an R package to implement the surrogate R 2 method for probit/logistic regression models. This package contains the R functions for generating the point and interval estimates of the surrogate R 2 measure. The point and interval estimates allow researchers and practitioners to evaluate the model’s overall goodness of fit and understand its uncertainty. In addition, we develop an R function that calculates the percentage contribution of each variable to the overall surrogate R 2 . This percentage reflects each variable’s contribution to the model’s overall explanatory power. Based on the contribution’s relative size, our R function provides an “importance” ranking of all the explanatory variables.

Second, to provide practical guidance for categorical data modeling, we use the developed R package to demonstrate how it can be used jointly with other R packages developed for variable screening/selection and model diagnostics ( leaps [28], step() function from the R core, glmnet [17], ordinalNet [40], ncvreg [5], grpreg [6], SIS [34], sure [18], PAsso [44]). In particular, we recommend a workflow that consists of three steps: variable screening/selection, model diagnostics, and goodness-of-fit analysis. The workflow is illustrated in the analysis of wine-tasting preference datasets.

Third, the comparability of the surrogate R 2 across different samples and/or models allows us to compare goodness-of-fit analysis from similar studies. The comparison can lead to additional scientific and business insights which may be useful for decision making. To illustrate this, we conduct goodness-of-fit analysis separately for the red wine and white wine samples to demonstrate the comparability of the surrogate R 2 . Our analysis result reveals that (i) the same set of explanatory variables has different explanatory power for red wine and white wine (43.8% versus 31.0%), and (ii) the importance ranking of the explanatory variable (in terms of their contribution to the surrogate R 2 ) is different between red wine and white wine.

Our SurrogateRsq package has broad applicability. It is compatible with the following R functions that can fit probit/logistic regression models for a binary or ordinal response: glm() in the R core, polr() in the MASS package [33], clm() in the ordinal package [9], and vglm() in the VGAM package [41], although we only demonstrate the functions of SurrogateRsq package through ordinal probit regression models using plor() in this paper. More examples and details can be found on the website: https://xiaorui.site/SurrogateRsq/.

We briefly review the surrogate R 2 measure in the study of [26]. For the model setting, we consider a probit/logit model with a set of explanatory variables. The categorical response is either a binary or ordinal variable Y that has J categories { 1 , 2 , … , J }, with the order 1 < 2 < ⋯ < J, (2.1) Pr { Y ≤ j } = G { α j − ( β 1 X 1 + ⋯ + β l X p ) } , j = 1 , … , J , where − ∞ < α 1 < ⋯ < α J < + ∞. The link function G ( · ) can be a probit ( G ( · ) = Φ ( · )) link or a logit ( G ( η ) = 1 / ( 1 + e − η )). Each generic symbol of { X 1 , … X p } in Model (2.1) can represent a single variable of interest, a high-order term (e.g., X 2 ), or an interaction term between X and another variable. It is well-known that an equivalent way to express Model (2.1) is through a latent variable. For example, if the link is probit, the latent variable has the following form with a normally distributed ϵ: Z = α 1 + β 1 X 1 + ⋯ + β p X p + ϵ , ϵ ∼ N ( 0 , 1 ) .

The categorical response Y can be viewed as generated from censoring the continuous latent variable Z in the following way: Y = 1 if − ∞ < Z ≤ α 1 + α 1 , 2 if α 1 + α 1 < Z ≤ α 2 + α 1 , ⋯ J if α J − 1 + α 1 < Z < + ∞ .

To construct a goodness-of-fit R 2 , [26] adopted the surrogate approach proposed by [24]. The idea of the surrogate approach is to simulate a continuous variable and use it as a surrogate for the original categorical variable in the analysis [24, 25, 8]. In the context of probit models, [26] proposed to generate a surrogate response variable using the following truncated conditional distribution: S ∼ Z ∣ − ∞ < Z ≤ α 1 + α 1 if Y = 1 , Z ∣ α 1 + α 1 < Z ≤ α 2 + α 1 if Y = 2 , ⋯ Z ∣ α J − 1 + α 1 < Z < + ∞ if Y = J .

[26] proposed to regress the surrogate response S on { X 1 , … , X p } using a linear model below: (2.2) S = α 1 + β 1 X 1 + ⋯ + β p X p + ϵ , ϵ ∼ N ( 0 , 1 ) . Their approach used the OLS R 2 measure of this linear model as a surrogate R 2 for Model (2.1): R ( S ) 2 { X 1 , … , X p } = the OLS R 2 of the linear model (2.2) .

[26] showed that the surrogate R ( S ) 2 measure has three desirable properties. First, it approximates the OLS R 2 calculated using the latent continuous outcome Z. This property enables us to compare surrogate R 2 ’s and OLS R 2 ’s across different models and samples that address the same scientific question. Second, as it is the OLS R 2 calculated using the continuous surrogate response S, the surrogate R ( S ) 2 has the interpretation of the explained proportion of variance. It measures the explained proportion of the variance of the surrogate response S through the linear model. This explained proportion of variance implies the explanatory power of all the features in the fitted model. Third, the surrogate R ( S ) 2 maintains monotonicity between nested models, which makes it suitable for comparing the relative explanatory power of different models. In contrast, the well-known McFadden’s R 2 does not preserve the first two properties of the surrogate R ( S ) 2 . McFadden’s R 2 relies on the ratio of likelihoods, so it neither approximates the OLS R 2 nor preserves the interpretation of explained variance. On the other hand, [26] showed that McKelvey-Zavoina’s R M Z 2 did not necessarily maintain monotonicity between nested models. This serious issue may make McKelvey-Zavoina’s R M Z 2 an unsuitable tool for measuring the goodness of fit.

To make inferences for the surrogate R ( S ) 2 , [26] provided procedures to produce point and interval estimates. Since the surrogate response S is obtained through simulation, [26] used a multiple-sampling scheme to “stabilize” the point estimate. They also provided an implementation to produce an interval estimate with a 95 % confidence level. This confidence interval is constructed through a bootstrap-based pseudo algorithm. When the sample size is large (e.g., n = 2000), [26]’s numerical studies show that the interval measure of the surrogate R ( S ) 2 can approximate the nominal coverage probability.

It is also worth noting that [26]’s method requires a full model. This paper will illustrate how to use existing tools, such as variable selection and model diagnostics, to initiate a full model. The full model is used to generate a common surrogate response S, which is then used to calculate the surrogate R ( S ) 2 ’s of whatever reduced models. We will demonstrate how to carry it out in a real data analysis presented in Section 5.

We develop an R package SurrogateRsq for goodness-of-fit analysis of probit models [43]. This package contains functions to provide (i) a point estimate of the surrogate R 2 ; (ii) an interval estimate of the surrogate R 2 ; (iii) an importance ranking of explanatory variables based on their contributions to the total surrogate R 2 of the full model; and (iv) other existing R 2 measures in the literature. In this section, we explicitly explain the inputs and outputs of these functions. In the next two sections, we will demonstrate the use of these functions through a recommended workflow and real data examples.

In empirical studies, goodness-of-fit analysis should be used jointly with other statistical tools, such as variable screening/selection and model diagnostics, in the model-building and refining process. In this section, we discuss how to follow the workflow in Figure 1 to carry out statistical modeling for categorical data. We also discuss how to use the SurrogateRsq package with other existing R packages to implement this workflow. As [26]’s method requires a full model, researchers and practitioners can also follow the process in Figure 1 to initiate a full model to facilitate goodness-of-fit analysis. 1.

In Step-0, we can use the AIC/BIC/LASSO or any other variable selection methods deemed appropriate to trim or prune the set of explanatory variables to a “manageable” size (e.g., less than 20). The goal is to eliminate irrelevant variables so researchers can better investigate the model structure and assessment. The variable selection techniques have been studied extensively in the literature. Specifically, one can implement (i) the best subset selection using the function regsubsets() in the leaps package; (ii) the forward/backward/stepwise selection using the function step() in the R core; (iii) the shrinkage methods including the (adaptive) LASSO in the glmnet package; (iv) the regularized ordinal regression model with an elastic net penalty in the ordinalNet package; and (v) the penalized regression models with minimax concave penalty (MCP) or smoothly clipped absolute deviation (SCAD) penalty in the ncvreg package [36, 46, 45, 35, 40]. When the dimension is ultrahigh, the sure independence screening method can be applied through the SIS package [16]. When the variables are grouped, one can apply the group selection methods including the group lasso, group MCP, and group SCAD through the grpreg package [6]. In some cases, Step-0 may be skipped if the experiment only involves a (small) set of controlled variables. In these cases, the controlled variables should be modeled regardless of statistical significance or predictive power. We limit our discussion here because our focus is on goodness-of-fit analysis.

In this section, we demonstrate how to use our SurrogateRsq package, coupled with R packages for model selection and diagnostics, to carry out statistical analysis of the wine rating data. A critical problem in wine analysis is to understand how the physicochemical properties of wines may influence human tasting preferences [10]. For this purpose, [10] collected a dataset that contains wine ratings for 1599 red wine samples and 4898 white wine samples. The response variable, wine ratings, is measured on an ordinal scale ranging from 0 (very bad) to 10 (excellent). The explanatory variables are 11 physicochemical features, including alcohol, sulphates, acidity, dioxide, pH, and others.

Our analysis of the wine rating data follows the workflow discussed in Section 4. Specifically, in Section 5.1, we initiate a full model using several R packages for variable selection and model diagnostics. In Section 5.2, we use our SurrogateRsq package to evaluate (i) the goodness-of-fit of the full model and several reduced models; (ii) the contribution of each individual variable to the overall R 2 ; and (iii) the difference between the red wine and white wine in terms of how physicochemical features may influence human tasting differently.

In this paper, we have developed the R package SurrogateRsq for categorical data goodness-of-fit analysis using the surrogate R 2 . The package applies to probit/logistic regression models, and it is compatible with commonly used R packages for binary and ordinal data analysis. With SurrogateRsq, we are able to obtain the point estimate and the interval estimates of the surrogate R 2 . An importance ranking table for all explanatory variables can be produced as well. These new features can be used in conjunction with other R packages developed for variable selection and model diagnostics. This “whole-analysis” is summarized in a workflow diagram, which can be followed in practice for categorical data analysis. To examine the utility of this package in real data analysis, we have used a wine rating dataset as an example and provided the sample codes. In addition, we have used the package SurrogateRsq to demonstrate that the surrogate R 2 allows us to compare different models trained from the red wine sample and white wine sample. The comparison has led to new findings and insights that deepen our understanding of how physicochemical features influence wine quality. The result suggests that our package can be used in a similar way to analyze multiple studies (and/or models) that address the same or similar scientific or business question.

We use the red wine data to examine the computational time of the functions in our package. Table 3 presents the comparison, where the column (n = 1597) corresponds to the real data, and the other columns (n = 3000, 6000, 12000) are based on pseudo-real data sets generated by randomly sampling more rows from the real data. The numbers are the average running time in seconds over 10 times of repetition on an Apple Macbook Pro Max with the M1 Max CPU. The upper panel of Table 3 shows if only a point estimate of the surrogate R 2 is needed, our surr_rsq() function takes almost no time to provide the result (e.g., merely 0.213 seconds when n = 12000). However, the bottom panel of Table 3 shows if confidence intervals are wanted, our surr_rsq_ci() function takes longer time with one core of CPU (e.g., 1421 seconds = 23.6 minutes when n = 12000). Given that a CPU with 6 10 cores is quite common nowadays, we recommend using parallel computing aforementioned. It reduces the computing time to 211.93 seconds = 3.5 minutes when n = 12000.

If software developers want to build or modify this package for their specific scientific inquiries, they can modify one or all of the three components of our package. First, what we really need as an input for the functions in our SurrogateRsq package is the fitted model from another model training package (e.g., glm(), polr()). Software developers can replace the object with the model of their interest. For example, this surrogate R 2 approach may still work for the discrete choice models studied by [8]. Second, depending on the form of the model, software developers can choose or modify what distribution to use for simulating the surrogate response. Third, once the surrogate responses are available, one can follow the inference procedures discussed in our paper and tailor them to meet specific needs.