Identifying the Big Shots—A Quantile-Matching Way in the Big Data Context

2.1 The P-value and Big Data

The p-value was first introduced by the UK statistician Ronald Fisher [17] in the 1920s to determine whether the observed data complied with the proposed hypotheses by calculating the difference between the predicted and the observed data series. Fisher regarded the p-value as an informal way of checking whether the evidence of a treatment effect was worth a second look. However, in the late 1920s, during the movement to make evidence-based decision-making more rigorous and objective, many non-statistician authors combined Fisher's easy-to-calculate p-value with Neyman and Pearson's reassuringly rigorous rule-based approach to create a hybrid system, and thus a p-value of 0.05 became enshrined as “statistically significant” [43].

The abuse and misinterpretation of the p-value were first noted in the late 20th century. Researchers were mainly concerned about two aspects. First, the p-value does not establish the probability as to whether the investigator's hypothesis is correct; it only represents the false-positive rate given the observations [15, 21, 22]. Second, using the p-value to divide the results into “significant” and “insignificant” categories is arbitrary [16]. More extensive concerns about the p-value have recently been voiced and are presented in Table 1. Among these concerns, the deflated p-value problem has received considerable attention. This problem is not a new phenomenon. In 1993, Kass and Raftery [30] pointed out “frequently tests tend to reject null hypotheses almost systematically in very large samples”, an observation that has been further elaborated by Lin et al. [39], Kim and Ji [33], and Kim et al. [32]. Because the p-value measures the distance between the data and the null hypothesis united by standard errors, as n approaches infinity, the standard error will be close to 0. This means even a tiny deviation from a given value can be detected.

The influence of the deflated p-value problem has grown as large datasets have become extensively applied in research. Lin et al. [39] surveyed articles in MIS Quarterly and Information Systems Research (ISR) from 2004 to 2010 along with abstracts from the Workshop on IS and Economics and symposia on Statistical Challenges in Electronic Commerce Research. They found over 50% of the papers used extremely large samples (>10,000) and solely relied on low p-values and the sign of the coefficient. Kim and Ji [33] also reported 42% of the articles published during 2012 in four finance journals exploited extremely large samples and used the p-value null hypothesis testing method only. Similarly, Kim et al. [32] found 39% of the studies published during 2014 in eight accounting journals used extremely large samples; each study used the p-value and only one reported the confidence interval. Additionally, we have reviewed empirical papers published in MIS Quarterly and ISR from 2014 to 2018 and found that although IS researchers were more cautious with a large sample—among papers exploiting extremely large samples—58% also reported the effect size (\( \theta \)) or confidence interval along with the p-value, and a few of them reported the predictive power [3, 49]. Furthermore, there is a slight decrease in the percentage of papers (48%) dealing with extremely large samples, but the average sample size increased to 4.8 million. These findings highlight the urgent need to develop another test statistic to overcome the deflated p-value problem in the large sample context.

2.2 Effect Size

Several methods have been proposed to deal with the deflated p-value problem, such as changing the significance threshold [9, 33], reporting the confidence interval of the estimated coefficient [4, 32, 36, 39], and reporting the effect size [32, 39]. Lowering the p-value threshold is easy to implement, but it cannot fully address the deflated p-value concern as the sample size can be easily increased, particularly if p-hacking behavior is considered [14]. A confidence interval describes the range that the true value of the unknown parameter would fall into with a certain probability, and the width of the interval decreases with the sample size, however, it can only be applied at the variable level. Effect size focuses on the size of an observed effect rather than whether the effect is there, providing a more quantitative description of the observed effect and a different statistical inference. Because variables will become significant unless they are strictly equal to zero in very large data samples, it is natural to use the effect size to draw statistical inferences in large sample research.

The existing effect size measures can be divided into two types: those specific to comparing two conditions (Cohen's d, Hedges’ g, Glass's d, etc.) and those describing the proportion of variability explained (\( {\eta ^2} \), \( \eta _p^2 \), R2, adjusted R2, etc). We found that the most reported effect size measure is the regression coefficient (\( \theta \))—among 58% of the IS papers published between 2014 and 2018 that reported effect size, nearly all of them used the regression coefficient (\( \theta \)). The regression coefficient (\( \theta \)) describes the sensitivity of the dependent variable to a unit change in the independent variable; in that sense, it belongs to the first type of effect size measure. Although the use of the regression coefficient (\( \theta \)) alleviates the deflated p-value concern to some extent, this method has limitations. First, the regression coefficient (\( \theta \)) is an unstandardized effect size measure whose value is influenced by the choice of scale. This hinders the comparison of effect size across different variables and different studies (standardized coefficient can help with this issue but it still faces the second limitation as pointed below). Second, the regression coefficient (\( \theta \)) method can only be applied at the variable level, not the factor level. In many cases, we are interested in the implication of the underlying factors, instead of variables. For instance, if we want to study the relationship between people's personalities and their job performance, where the personality is measured from five dimensions with many items, we are interested in the effect size of the Big Five personality factors in explaining peoples’ job performance, rather than the effect size of each item. In this case, the traditional effect size (\( \theta \)) won't be able to provide many insights. In addition, our e-mail marketing case in Section 4 indicates the effect size for the Seasonal/Festive variable increases the e-mail opening odds by 27%. But the effect sizes for Membership Tier and Message Type cannot be obtained because they are categorical factors and are operationalized as several variables.

In the current study, we have used the LRT statistic to enable the application of effect size in multivariate situations. This test statistic directly measures the joint effect size of all the variables corresponding to the focal factor by calculating the unique amount of explanatory power provided by the focal factor to the dependent variable. Furthermore, we have developed a quantile matching method and an absolute multivariate effect size test statistic to enable the comparison of effect size across different factors and different studies to overcome the scaling effect.

3.1 Log-likelihood Ratio Test Statistic and Multivariate Effect Size

When comparing the relative tenability of two competing nested models, the most common method is the log likelihood ratio test, where \( \begin{equation*} {\rm{\Lambda }} = \frac{{{\rm{max}}[{L_0}(Null{\rm{\ }}Model|Data)]}}{{{\rm{max}}[{L_1}(Alternative{\rm{\ }}Model|Data)]}}. \end{equation*} \)Then \( - 2log{\rm{\Lambda }} \) has an asymptotic \( {\chi ^2} \) distribution with q degrees of freedom, where q is the difference in the number of free parameters between the general and restricted hypotheses [29]. By testing whether \( -2log{\rm{\Lambda }} \) is significantly different from 0, we can tell whether the alternative model is preferred or not. Similarly, we apply this method to identify the important factors by constructing the LRT statistic \( LR{T_i} \): (1) \( \begin{equation} LR{T_i} = - 2( {{\cal L}( {{\boldsymbol{\tilde{\theta }}}}) - {\cal L}( {{\boldsymbol{\hat{\theta }}}})}), \end{equation} \)where \( {\cal L}( {\boldsymbol{\theta }} ) \) is the log likelihood of \( {\boldsymbol{\theta }} \), \( {\boldsymbol{\hat{\theta }}} = ( {{{\hat{\theta }}_1},{{\hat{\theta }}_2}, \ldots {{\hat{\theta }}_v},{{\hat{\theta }}_{v + 1}}, \ldots {{\hat{\theta }}_w}} ) \) is the maximum likelihood (ML) estimator for \( {\boldsymbol{\theta }} \) over the full parameter set, and \( {\boldsymbol{\tilde{\theta }}} = ( {{{\tilde{\theta }}_1},{{\tilde{\theta }}_2}, \ldots {{\tilde{\theta }}_v},0, \ldots 0} ) \) is the ML estimator under the null hypothesis \( {H_0} \): the exclusion of \( {X_i} \) has no impact on the model's goodness of fit, i.e., \( \ {\theta _{v + 1}} = \cdots = {\theta _w} = 0 \). Then \( LR{T_i} \) measures the impact of the exclusion of \( {X_i} \) from the model on the model's goodness of fit, which can be used as an indicator of the effect size of \( {X_i} \).

Based on Equation (1), the LRT statistic can be rewritten as: \( \begin{equation*} LR{T_i} = - 2( {{\cal L}( {{\boldsymbol{\tilde{\theta }}}} ) - {\cal L}( {{\boldsymbol{\hat{\theta }}}})}) \approx - \mathop \sum \limits_{i = v + 1}^w \frac{{{\partial ^2}{\cal L}}}{{\partial {{\hat{\theta }}_i}\partial {{\hat{\theta }}_i}}}\hat{\theta }_i^2 = - \mathop \sum \limits_{i = v + 1}^w {\lambda _i}p_{ii}^2\hat{\theta }_i^2. \end{equation*} \)(Proof is presented in Appendix A.) Therefore, the LRT statistic measures the overall effect size for the focal factor by summing the square of the independent coefficient \( \theta \) for each variable weighted by the second derivative of the log likelihood on \( \theta \). This can be interpreted as the amount of information the variable carries about \( \theta \) or the credibility of \( \theta \). Note that the second derivative also equals \( {\lambda _i}p_{ii}^2 \), in which \( {\lambda _i} \) is the eigenvalue and \( {p_{ii}} \) is the value from the eigenvector matrix of the corresponding Hessian matrix, as eigenvectors give the directions in which \( \theta \) increases or decreases the most and eigenvalues give the magnitudes of those changes in \( \theta \), the LRT statistic shows a joint effect from individual \( \hat{\theta }_i^2 \). Through the eigenvalue–eigenvector pairs, the LRT statistic helps address the issue of aggregating multiple univariate effect sizes. Therefore, we have used the LRT statistic as a basis to develop our test statistic. When \( w - v = 1 \) (the factor is operationalized as one variable), the LRT statistic gives us the univariate effect size by calculating the square of its corresponding coefficient weighted by its credibility. That is, the LRT statistic is closely related to the traditional effect size (\( \theta \)) in the univariate setting. But it can be easily extended into a multivariate setting. We have also provided numerical evidence for the close relationship between the LRT statistic and traditional effect size (\( \theta \)) in the univariate setting in Table 9, Section 4— the correlation between the LRT statistic and the square of coefficient (\( \theta \)) is extremely close to 1. By applying the LRT method in our email marketing case, we can easily obtain the multivariate effect size for Membership Tier and Message Type—44,733.9 and 39,527.7, respectively. This can be considered as a joint effect from all the variables corresponding to these two factors.

Although the LRT statistic can help quantify the joint effect from multiple variables, its expected value and uncertainty are determined by its degrees of freedom and may not enable the comparison of the two factors unless both factors have the same degrees of freedom. In the next section, we will propose the following quantile-matching method to facilitate the comparison of factors with different degrees of freedom and provide theoretical justifications for quantile matching.

3.2 Quantile-Matching Transformation Method for Multivariate Effect Size

As discussed, we cannot compare the effect size of factors directly from the LRT statistic because more variables in a factor will lead to a higher expected explanatory power and greater variance. To resolve this problem, we must adjust the effect of the number of variables underlying the factor of interest. Alternative methods to achieve this have been proposed: (1) using cumulative probability, which provides a concise and universal summary of different statistical tests that no other statistic can achieve [37]; and (2) matching the representative characteristics of the distribution. Four characteristics are typically used to describe a distribution: location, spread, skewness, and peakedness. The location can be measured by the mean or median and the SD or interquartile ranges can be used to describe the spread. The classic method of matching both mean and SD (location and spread) of the distributions of the LRT statistic can be achieved by: \( \begin{equation*} {Z_i} = \frac{{{{\widehat {LRT}}_i} - E\left( {{{\widehat {LRT}}_i}} \right)}}{{\sqrt {Var\left( {{{\widehat {LRT}}_i}} \right)} }} = \frac{{{{\widehat {LRT}}_i} - d}}{{\sqrt {2d} }}, \end{equation*} \)where the mean of \( {Z_i} \) is zero and its variance is 1.

However, both methods above have major limitations, particularly in the large sample context. Although cumulative probability is a powerful tool for comparisons across contexts, the prevalence of a large sample significantly limits its applicability. The pchisq() function in RStudio (1.1.383) can display a maximum cumulative probability of 1–10⁻³²³. However, the cumulative probability for a \( {\chi ^2} \) test statistic larger than 1,590 and with degrees of freedom of less than 20 will be calculated as 1. Such a large test statistic is quite common with a large sample: in the example in Section 4, with n = 1,565,720, the LRT statistics for all factors exceeded 1,590. Matching several representative characteristics of the distribution method only uses the aggregated information of the distribution (e.g., location and spread), leaving other information about the distribution unused. This leads to comparison inconsistency when using the cumulative probability method. In the example presented in Section 1, for the two test statistics \( LR{T_1} = 21.11\sim{\chi ^2}( 3 ) \) and \( LR{T_2} = 25.75\sim{\chi ^2}( 5 ) \), the corresponding cumulative probabilities for \( LR{T_1} \) and \( LR{T_2} \) are both 0.9999, whereas the transformed mean/SD statistics \( {Z_1} \) and \( {Z_2} \) are 7.39 and 6.56, respectively, leading us to wrongly claim the first factor \( {F_1} \) has a larger effect size.

In this study, we will propose a method of adjusting the impact of the number of variables underlying the factor of interest based on matching across quantiles of distributions. This quantile matching transformation method can achieve the three Fs: (1) Feasibility—it has a one-to-one corresponding relationship with cumulative probabilities under a perfect matching situation, but it can also be applied to scenarios where the cumulative probability method fails due to the use of extremely large samples; (2) Flexibility—the quantile-matching transformation may be done on the full distribution or a specific interval, providing tremendous flexibility in real applications; and (3) Fast calculation—it is easy to implement and calculate.

The quantile-matching method is as follows: let M be a random variable and \( {\boldsymbol{N}} = ( {{N_1}, \ldots ,{N_p}} ) \) be a collection of p random variables. The goal for quantile-matching is to find a linear combination \( {\boldsymbol{\beta 'N}} \) where \( {\boldsymbol{\beta }} \) minimizes the integrated squared difference between the quantile functions of M and \( {\boldsymbol{\beta 'N}} \) across [0, 1] [50]. The quantile function is the inverse function of the cumulative distribution function, suppose \( {Q_\xi }( \alpha ) \) denotes the \( \alpha \)th quantile of the random variable \( \xi \), then \( {\rm{P}}\{ {\xi \le {Q_\xi }( \alpha )} \} = \alpha ,{\rm{\ for\ }}\alpha \in [ {0,1} ] \). In addition, if only a certain quantile interval \( [ {{\alpha _1},{\alpha _2}} ] \) is of interest, the integration interval can be changed to \( [ {{\alpha _1},{\alpha _2}} ] \) instead of [0, 1].

Note that the two LRT statistics \( {\widehat {LRT}_1} \) and \( {\widehat {LRT}_2} \) can be rewritten as the quantile functions \( {Q_{{k_1}}}( {{q_1}} ) \) and \( {Q_{{k_2}}}( {{q_2}} ) \), where \( {Q_k}( \alpha ) \) denotes the \( \alpha \)th quantile of the random variable, \( {q_1} \) and \( {q_2} \) are the corresponding cumulative probabilities for \( {\widehat {LRT}_1} \) and \( {\widehat{LRT}_2} \), and \( {k_1} \) and \( {k_2} \) are the corresponding degrees of freedom of \( {\widehat {LRT}_1} \) and \( {\widehat {LRT}_2} \). In this article, we propose the transformed LRT statistics, \( Z_1^q \) and \( Z_2^q \), based on the linear quantile-matching transformation to adjust the effect of the number of variables or the effect of the degrees of freedom in different factors. We define \( Z_1^q \) and \( Z_2^q \) as the functions of the cumulative probability, q, as \( \frac{{{Q_{{k_1}}}( q ) - {a_{{k_1}}}}}{{{b_{{k_1}}}}} \) and \( \frac{{{Q_{{k_2}}}( q ) - {a_{{k_2}}}}}{{{b_{{k_2}}}}} \), where in practice, \( {a_{{k_1}}} \), \( {b_{{k_1}}} \), \( {a_{{k_2}}} \), \( {b_{{k_2}}} \) are the estimated values obtained by minimizing the integrated squared difference between the quantile functions. In an ideal situation where we can conduct perfect matching at \( {q_1} \) and \( {q_2} \) for two LRT statistics—\( {\widehat {LRT}_1} \) and \( {\widehat {LRT}_2} \)—with degrees of freedom \( {k_1} \) and \( {k_2} \), we have an exact match in \( Z_1^q \) and \( Z_2^q \) at the two cumulative probabilities \( {q_1} \) and \( {q_2} \): (2) \( \begin{equation} \frac{{{Q_{{k_1}}}\left( {{q_1}} \right) - {a_{{k_1}}}}}{{{b_{{k_1}}}}} = \frac{{{Q_{{k_2}}}\left( {{q_1}} \right) - {a_{{k_2}}}}}{{{b_{{k_2}}}}}, \end{equation} \)and: (3) \( \begin{equation} \frac{{{Q_{{k_1}}}\left( {{q_2}} \right) - {a_{{k_1}}}}}{{{b_{{k_1}}}}} = \frac{{{Q_{{k_2}}}\left( {{q_2}} \right) - {a_{{k_2}}}}}{{{b_{{k_2}}}}}. \end{equation} \)

We will demonstrate in the following that under an ideal situation that leads to Equations (2) and (3), the transformed test statistic \( Z_i^q \) will have a one-to-one corresponding relationship with the cumulative probability. Specifically, the following relationship will hold: \( \begin{equation*} If{\rm{\ }}{q_1} > {q_2},{\rm{\ }}then{\rm{\ }}Z_1^q > Z_2^q,{\rm{\ }}and{\rm{\ }}vice{\rm{\ }}versa. \end{equation*} \)

If \( {q_1} > {q_2} \), then \( {Q_{{k_1}}}( {{q_1}} ) > {Q_{{k_1}}}( {{q_2}} ) \), and according to Equation (3) we have: \( \begin{equation*} \frac{{{Q_{{k_1}}}\left( {{q_1}} \right) - {a_{{k_1}}}}}{{{b_{{k_1}}}}} > \ \ \frac{{{Q_{{k_1}}}\left( {{q_2}} \right) - {a_{{k_1}}}}}{{{b_{{k_1}}}}} = \frac{{{Q_{{k_2}}}\left( {{q_2}} \right) - {a_{{k_2}}}}}{{{b_{{k_2}}}}} \end{equation*} \)that is, \( Z_1^q > Z_2^q \).

If \( Z_1^q > Z_2^q \), then: \( \begin{equation*} \frac{{{Q_{{k_1}}}\left( {{q_1}} \right) - {a_{{k_1}}}}}{{{b_{{k_1}}}}} > \ \ \frac{{{Q_{{k_2}}}\left( {{q_2}} \right) - {a_{{k_2}}}}}{{{b_{{k_2}}}}} \end{equation*} \)and according to Equation (3), we have: \( \begin{equation*} \frac{{{Q_{{k_1}}}\left( {{q_1}} \right) - {a_{{k_1}}}}}{{{b_{{k_1}}}}} > \ \ \frac{{{Q_{{k_1}}}\left( {{q_2}} \right) - {a_{{k_1}}}}}{{{b_{{k_1}}}}} \end{equation*} \)therefore, \( {q_1} > {q_2} \). Therefore, the comparison based on \( Z_1^q \) and \( Z_2^q \) will be identical to the comparison based on \( {q_1} \) and \( {q_2} \). We recognize that in a real application, Equations (2) and (3) will not always hold perfectly since (1) we have matched the quantiles across distributions numerically (shown in Section 3.3), so the quantile matching results may be less accurate when a large step size is used; and (2) the quantile cannot be calculated explicitly when the sample size is very large, and therefore an explicit match on those points cannot be achieved. Although the one-to-one corresponding relationship between the cumulative probability and the quantile-matching-based transformed statistic can be violated to a small extent, the relationship roughly holds.

3.2 The Implementation of Quantile-Matching Transformation for LRT Statistics

The above quantile-matching method can be applied to any type of distribution. In our study, we are specifically interested in \( {\chi ^2} \) distributions because our LRT statistic \( {\widehat {LRT}_i} \) follows an asymptotic \( {\chi ^2} \) distribution. We have attempted to find a linear relationship between \( {\chi ^2} \) distributions, with degrees of freedom from 2 to 20 with \( \chi _{( 1 )}^2 \) according to the following matching formula: (4) \( \begin{equation} \chi _{\left( n \right)}^2 \approx a + b\chi _{\left( 1 \right)}^2 \end{equation} \)by minimizing the following integrated squared difference of the two quantile functions: \( \begin{equation*} \mathop \int \limits_0^1 {\left\{ {{Q_{\chi _{\left( n \right)}^2}}\left( \alpha \right) - {Q_{a + b\chi _{\left( 1 \right)}^2}}\left( \alpha \right)} \right\}^2}d\alpha \end{equation*} \)

We have chosen a linear model specification due to its simplicity and our numerical estimation have shown it can also achieve a reasonable model fit (the adjusted R² is above 0.8 for all analyses). We have used the numerical method proposed by Sgouropoulos et al. [50] to estimate \( a \) and \( b \), which provides a better fitting at the tails of the distributions and works as follows: \( \begin{equation*} ( {\hat{a},\hat{b}} ) = \arg \mathop {\min }\limits_{a,b} \mathop \sum \limits_{j = v}^{mv} {\left\{ {{Q_{\chi _{\left( n \right)}^2}}\left( j \right) - {Q_{a + b\chi _{\left( 1 \right)}^2}}\left( j \right)} \right\}^2}, \end{equation*} \)where \( v \) denotes the stepwise increment used to divide the quantile interval.

In our study, we have chosen the quantile ranging from [0.000001,0.999999] with a step size of 1E-6. We have also excluded 0 and 1 from the range because \( {Q_{\chi _{( n )}^2}}( 0 ) = 0 \) and \( {Q_{\chi _{( n )}^2}}( 1 ) = Inf \) for all n, as the inclusion of infinity would cause problems in the integration.

Then, if \( n \ge 2 \), \( {\widehat {LRT}_i} \) can be transformed in the following way: (5) \( \begin{equation} Z_i^q = \frac{{{{\widehat {LRT}}_i} - {{\hat{a}}_n}}}{{{{\hat{b}}_n}}} \end{equation} \)

The estimated \( {\hat{a}_n} \) and \( {\hat{b}_n} \) for each n are summarized in Table 2:

We have compared the transformed test statistic using the mean/SD method and the QM transformation method at the 0.95 quantiles of \( {\chi ^2} \) distribution with degrees of freedom from 1 to 10. The results in Table 3 show that the QM transformation method provides more consistent results with the cumulative probability than the mean/SD method—for the QM transformation, the difference from \( Z_1^q \) remains below 0.5%, whereas the difference from \( Z_1^{mean/SD} \) quickly surpasses 7% for the mean/SD method.

To address the concern where the p-value deflates quickly to zero in the large sample context and achieve a high matching accuracy, we have repeated our analysis in the following two intervals according to Equation (4): [0.99,0.9999999] with step size 1E-7 and the extreme scenario [1–10⁻³²⁰,1–10⁻³²²] with step size 10⁻³²⁰. The estimated \( {\hat{a}_n} \) and \( {\hat{b}_n} \) for each n based on [0.99,0.9999999] are summarized in Table 4. For the 0.9999 quantile of the \( {\chi ^2} \) distribution with degrees of freedom from 1 to 10, we have compared the transformed test statistic using the mean/SD method, the QM transformation method (full distribution), and the QM transformation method ([0.99,0.9999999]). Results are presented in Table 5. We can see that the QM transformation method using the [0.99,0.9999999] interval produces the most consistent results with the cumulative probability method. The reason is that the \( {\chi ^2} \) test statistic changes much more quickly at the right tail with the same unit change in the cumulative probability; full distribution matching would underestimate this change and achieve lower accuracy, as other parts of the distribution are also taken into consideration, whereas interval matching does not have this problem. This demonstrates the advantage of quantile-matching-based transformation—\( {\hat{a}_n} \) and \( {\hat{b}_n} \) can be calculated according to different integration intervals, which provides more flexibility for different scenarios.

Then, we have tried repeating our analysis according to Equation (4) in the extreme intervals, e.g., [1–10⁻³²⁰,1–10⁻³²²]. However, as the corresponding \( {\chi ^2} \) statistics are not calculatable for quantiles larger than 1-10-320 in RStudio, direct matching is not possible. To address this issue, we have used QM and natural spline interpolation methods to obtain the estimated \( {\chi ^2} \) statistics for those quantiles by the following steps [45]. Details are presented in Appendix B:

• We have transformed the probability distribution function (PDF) of the \( {\chi ^2} \) distribution to obtain a second-order ordinary differential Equation (ODE) by using the quantile mechanic approach.

• Power series approach is used to obtain the estimates of the first few items in ODE, which are served as the initial approximates.

• Natural spline interpolation is used to force the initial approximates to converge to a solution close to the R values.

• The obtained spline function is then used to estimate the \( {\chi ^2} \) statistics in the extreme interval.

Then, the estimated \( {\chi ^2} \) statistics are used for quantile matching by Equation (4) and the estimated \( {\hat{a}_n} \) and \( {\hat{b}_n} \) are presented in Table 6. In practice, we suggest IS researchers to use the following criteria to decide how the transformation should be done:

• If the LRT statistics exceed 1,500, Table 6 should be used.

• If the corresponding quantile of the LRT statistic is larger than 0.99, Table 4 should be used; otherwise, Table 2 is sufficient.

The idea of the transformation method is essential. It ensures comparison accuracy and flexibility for a multivariate effect size across different factors with different numbers of operationalized variables. As shown in Table 18, in the context of the factors that influence accident severity, the LRT statistics for Location and Time are 24,522.13 and 8,699.50; because they are operationalized as 11 and one variables, theoretically their LRT statistics cannot be directly used in comparison and their cumulative probabilities are computed as 1 in RStudio. If we use the standardized LRT statistics for comparison, the transformed statistic for Location is \( \frac{{24,522.13 - 11}}{{\sqrt {22} }} = 5,\!225.79 \), while the transformed statistic for Time is \( \frac{{8,699.50 - 1}}{{\sqrt 2 }} = 6,\!150.77 \). We will conclude Time is more influential in explaining accident severity than Location. However, if we apply the quantile-matching transformation method, then the transformed statistic for Location is (24,522.13-51.6892)/1.0064 = 24314.83, which is much larger than 8,699.50. That is, in fact, Location is more influential in explaining accident severity than Time.

3.4 The Development of Absolute Multivariate Effect Size

The method outlined above enables us to identify the influential factors within one model. However, sometimes it would be more meaningful to conclude whether the factor is influential based on an absolute threshold. To do that, we have used the proportion of explanatory power provided by the focal factor similar to R2, (6) \( \begin{equation} {p_i} = \frac{{LR{T_i}}}{{LR{T_{\boldsymbol{X}}}}}, \end{equation} \)where \( LR{T_i} \) is the explanatory power provided by \( {X_i} \), and \( LR{T_{\boldsymbol{X}}} \) is the total explanatory power provided by all Xs.

However, because one factor can be operationalized as several variables and a greater number of variables provides more freedom for the model to estimate the variance in the dependent variable [20, 44], a factor operationalized as more variables tend to achieve a higher absolute multivariate effect size using Equation (6). To adjust the bias caused by the number of operationalized variables, we have used the quantile-matching transformed statistic to calculate the absolute multivariate effect size, which is similar to the adjusted R²,

(7)

Moreover, we can provide the thresholds for researchers to refer to by building connections with Cohen's \( {f^2} \). Cohen's \( {f^2} \) is defined in the following way: \( \begin{equation*} {f^2} = \frac{{{R^2}}}{{1 - {R^2}}}, \end{equation*} \)which is the ratio between the variance explained by the IVs and the variance explained by residuals. According to Cohen's [13] guidelines, \( {f^2} \ge 0.02,\ {f^2} \ge 0.15 \) and \( {f^2} \ge 0.35 \) represent small, medium, and large effect sizes, respectively, and the corresponding \( {R^2} \) is 0.196, 0.36, and 0.51. Because \( {R^2} \) is similar to \( {p_i} \) (as one measures the proportion of variance explained, and the other measures the proportion of explanatory power provided by IVs), we have adopted the same thresholds for \( {p_i} \). Furthermore, the \( Adjusted\ {p_i} \) can be rewritten as \( \begin{equation*} Adjusted\ {p_i} = \frac{{Z_i^q}}{{Z_{\boldsymbol{X}}^q}} = \left\{ \begin{array}{@{}*{1}{c}@{}} {{{\hat{b}}_N} \cdot \frac{{{{\widehat {LRT}}_i}}}{{{{\widehat {LRT}}_{\boldsymbol{X}}} - {{\hat{a}}_N}}},\ \ n = 1}\\[6pt] {{{\hat{b}}_N} \cdot \frac{{\frac{{{{\widehat {LRT}}_i} - {{\hat{a}}_n}}}{{{{\hat{b}}_n}}}}}{{{{\widehat {LRT}}_{\boldsymbol{X}}} - {{\hat{a}}_N}}},\ \ n > 1} \end{array}\right., \end{equation*} \)where n is the number of variables by which the focal factor \( {X_i} \) is operationalized and N is the number of IVs in the full model. Because \( {\hat{a}_N} \) is relatively small compared to \( {\widehat {LRT}_{\boldsymbol{X}}} \), \( Adjusted\ {p_i} \) is approximately equal to \( {\hat{b}_N} \cdot {p_i} \) when \( n = 1 \), and \( {\hat{b}_N} \cdot {\tilde{p}_i} \) when \( n > 1 \), where \( {\tilde{p}_i} = \) \( \frac{{\frac{{{{\widehat {LRT}}_i} - {{\hat{a}}_n}}}{{{{\hat{b}}_n}}}}}{{{{\widehat {LRT}}_{\boldsymbol{X}}}}} \) is the proportion of explanatory power provided by the focal factor after adjusting its number of operationalized variables. Then, the new thresholds for \( Adjusted\ {p_i} \) can be obtained by multiplying the corresponding \( {\hat{b}_N} \) in Tables 2, 4, or 6 with the degrees of freedom equal to the number of IVs in the full model.

4.1 Comparison with Predictively Influential Factors

Although explanation is central to theory development, methods for assessing the predictive claims of a model have received increasing attention in recent years. As a result, we have also compared the influential factors identified by our method with those identified by the predictive measurement of the area under the ROC curve (AUROC).

First, we will demonstrate below how we have identified influential factors based on their predictive accuracy by applying ROC analysis. After fitting a logistic regression model: (8) \( \begin{equation} \log \left( {\frac{p}{{1 - p}}} \right) = {\theta _0} + {\theta _1}{X_1} + {\theta _2}{X_2} + \cdots {\theta _k}{X_k}, \end{equation} \)where p = Pr(Y = 1), we can estimate the value of p, denoted by \( \hat{p} \), representing a score for prediction. Predicted values of Y can be derived by setting a “threshold value” \( v \) . An observation is predicted as “Y = 1” if \( \hat{p} \ge v \), where \( \hat{p} \) is calculated by substituting \( {\theta _i} \) with the ML estimate \( {\hat{\theta }_i} \).

A ROC curve is commonly used to assess the predictive performance of a binary regression, which is a plot of sensitivity (true positive rate) versus specificity (1 – false positive rate). By varying the threshold value \( v \) and calculating the AUROC, we can measure a model's predictive performance and the probability of making a correct binary classification [8, 26]. In practice, the calculation of AUROC is based on sample misclassification errors. To obtain a reliable estimate of AUROC, we have used N-fold cross-validation [27]. For example, with N = 10, we have randomly separated the data set into 10 non-overlapping blocks. To predict the category for observations in block j, j = 1, …, 10, we have obtained \( {\hat{\theta }_i} \) to determine the formula for the score \( \hat{p} \) using the other nine blocks and compute the AUROC(j). The AUROC is then estimated from the sample mean \( \mathop \sum \nolimits_{j = 1}^{10} AUROC( j )/10 \).

Using cross-validation, we have obtained \( AURO{C_{full}} \) for the full model. To assess the effect size of \( {X_i} \), we have refitted the logistic regression model in (8) without \( {X_i} \) in each block j and recalculate \( AURO{C_i} \). The level of influence is the difference in AUROC: \( \begin{equation*} DAURO{C_i} = AURO{C_{full}} - AURO{C_i}. \end{equation*} \)

In our example, we have found Previous E-mail Opened is the most influential factor—even from the perspective of prediction—followed by Message Type and Membership Tier. Less predictive power results came from Seasonal/Festive and Company Name. We have also replicated the ROC analysis with various sample sizes and have reported the results in Tables 12 and 13; the main results remain unchanged. However, we have noted that the influence ranking of factors is not stable when the sample size is too small (below 10%).

Next, we will compare the results obtained from the LRT and the AUROC analyses and have found both sets of results are highly correlated. First, we have calculated the Pearson correlation of the test statistic value produced by these two methods (presented in Table 14). We have discovered the Pearson correlation is above 0.999 across all sample sizes. We then calculate the Spearman correlation of the two sets of rankings. The correlation decreases, as the AUROC analysis identifies Message Type as the second most influential factor, whereas the LRT analysis indicates that Membership Tier is the second most influential. However, the correlation is still above 0.8 when the AUROC analysis results become stable (as shown in Table 15).

Initially, our finding may seem surprising as there should be a disparity between the ability to explain phenomena at the conceptual level and the ability to generate predictions at the measurable level. This is because explaining and predicting are different as the type of uncertainty associated with explanation is of a different nature than that associated with prediction [28]. According to Shmueli [51], measurable data are not accurate representations of their underlying constructs, therefore, operationalization of theories and constructs into statistical models and measurable data creates a disparity between the ability to explain phenomena at the conceptual level and the ability to generate predictions at the measurable level. That is, the results of the explanation are usually different from those of the prediction. However, Konishi and Kitagawa [35] have pointed out there may be no significant difference between inferring the true structure and making a prediction if an infinitely large quantity of data is available or if the data are noiseless. Our results have provided empirical support for this view and this demonstration is particularly relevant to big data research where predictive performance is the main objective when building statistical models. Furthermore, LRT analysis substantially outperforms ROC analysis in terms of computation time (shown in Table 16). On average, ROC analysis is around 100 times more time-consuming than our method.

4.2 Replication on the US Accidents Dataset

We have replicated our analysis on a public dataset—the US accidents dataset. This is a countrywide car accident dataset, which covers 49 states of the USA and was collected by Moosavi et al. [41]. The accident data were collected from February 2016 to Dec 2020 and include 47 variables. After excluding the records that contain missing values, there are 1,464,180 records remaining.

Traffic accidents are a major public health problem. Helping to predict their causes or occurrences have been a topic of interest in the field of machine learning. Several studies have used Moosavi's accident dataset to study this issue: Moosavi et al. [40] uses a deep-neural network model to predict real-time traffic accidents based on traffic events, weather data, time, and points-of-interest. They evaluated their model with the F1-score defined as \( \frac{{2{\rm{*}}precision{\rm{*}}recall}}{{precision + recall}} \), and demonstrated a significant improvement for class of accidents. Kebede [31] used the image data in this dataset to predict locations where car accidents are likely to happen with transfer learning with the CNN method. Parra et al. [46] evaluated different explainable machine learning models (e.g., random forests and decision trees) on this dataset in predicting road traffic crashes. Brodeur et al. [10] examined the impact of COVID-19 safer-at-home policies on car crashes with this dataset and found a 20% reduction in vehicular collisions.

In our example, we have tried understanding how the severity of an accident is influenced by the following three factors: weather, location, and time. The original accident severity in the dataset has four levels. For demonstration purposes, we have transformed it to a binary variable by considering levels 1 and 2 as less severe and levels 3 and 4 as very severe, such that we could apply logit model in this situation. Weather is described by temperature, wind_chill, humidity, pressure, pressure, visibility, wind_speed, and precipitation, Location is described by bump, crossing, give_way, junction, no_exit, railway, roundabout, station, stop, traffic calming, traffic signal, and turning loop, Time is indicated by the variable sunrise_sunset (whether it is day or night time). A simple logit model is built, and the results are presented in Table 17. Again, clear evidence of the deflated p-value problem is observed as the p-values for most variables are \( < 2E - 16 \).

First, we have used the LRT method to obtain the raw multivariate effect size at the factor level, then applied quantile-matching transformation using Table 6. The results are presented in Table 18. We have found that the most influential factor in explaining accident severity is Location, followed by Weather. We have also calculated the absolute multivariate effect size for each factor according to Equation (7) and have presented the results in Table 18, we have found Location provides about 47% explanatory power of the whole model and Weather provides about 31% explanatory power. As the thresholds for small, medium, and large effect sizes are 0.198, 0.364, and 0.515, we can conclude Location has a medium absolute effect size, Weather has a small effect size, and Time has a lower than small effect size.

Furthermore, we have performed the ROC analysis with 10-folder cross-validation and calculated the DAUROC for each factor—the results are presented in Table 18. We have found Weather and Location are more influential in predicting accident severity, which is consistent with the results of the multivariate effect size method, the Person correlation is 0.9401, and the Spearman correlation is 1.

Notice that if we use the standardized LRT statistics for comparison of the influence level of factors, we will conclude that Time is the most influential factor in explaining accident severity, followed by Location and Weather. However, the AUROC method suggests Location is the most influential factor in predicting accident severity, followed by Weather and Time, which is consistent with results generated by the transformed multivariate effect size method. This highlights the importance of the quantile-matching transformation technique in the large sample context. In addition, we have observed the direct use of raw LRT statistics can generate similar results as quantile-matching transformation. However, the direct use of raw LRT statistics lacks theoretical grounding. Therefore, we recommend using the quantile-matching transformation method in general for easy and fast calculation.

4.3 Replication of the Airbnb Listing Dataset

We have used the logit regression model in the previous two examples. In the following example, we have replicated our analysis on the Airbnb listing dataset¹ with a linear regression model. This dataset contains 250,000+ listings in ten major cities, including information about hosts, pricing, location, room type, and review scores.

Price determinants of Airbnb listings have always been of interest to researchers, given the uniqueness of Airbnb's listings and the heterogeneity of the landlord [56]. Wang and Nicolau [53] explored the impacts of the five factors—host attributes, site and property attributes, amenities and services, rental rules, and online review ratings—on prices using 180,533 accommodation rental offers in 33 cities. Cai et al. [11] examined the impacts of five groups of explanatory variables on Airbnb prices in Hong Kong: listing attributes, host attributes, rental policies, listing reputation, and listing location. Wu and Qiu [2019] focused on the influence of nice factors (external factors, landlord characteristics, location characteristics, listing characteristics, room facilities, rental rules, trust, sociality, and tenant characteristics) on listings prices using 51,874 listings in 36 cities in China. Similarly, in this example, we have investigated the impact of host attributes, property attributes, location, and online review ratings on listing prices (the prices are converted to Euros). Host attributes include the duration the host has been on Airbnb, whether the host's identity is verified, whether the host is a “Superhost”, and the total number of listings the host has on Airbnb. The property attribute describes the property as an entire place, a private room, or a shared room. Location refers to the city the property is in. Online review ratings include the overall rating and the respective rating of cleanness, check-in experience, communication experience with the host, location within the city, and listing's value relative to its price. There are 161,447 listings remaining after cleaning the data. A simple linear model is built, and the results are presented in Table 19. Our observation of the p-values for quite a few variables is <2e-16, indicating the deflated p-value problem.

The raw multivariate effect size of each factor is first obtained through the LRT method, then quantile-matching transformation is applied using Tables 4 and 6 (Table 4 is used for Host Attributes and Online Review Ratings, and Table 6 is used for Property Attributes and Location), and the results are presented in Table 20. We find that the most influential factors in explaining Airbnb listing prices are Property Attributes and Location. The absolute multivariate effect size for each factor is also calculated according to Equation (7) and is presented in Table 20. Property Attributes and Location have medium absolute effect sizes (both can provide about 42% explanatory power of the whole model), while the effect sizes of Host Attributes and Online Review Ratings is lower than small.

Furthermore, we have performed the Root-Mean-Square Error (RMSE) analysis with 10-folder cross-validation and calculated the Difference in RMSE for each factor. The results are presented in Table 20. RMSE is a measure of the differences between values predicted by a model and the values observed. It is commonly used to measure the prediction performance of linear regression models. RMSE is calculated as follows: \( \begin{equation*} RMSE = {\rm{\ }}\sqrt {\mathop \sum \limits_{i = 1}^n \frac{{{{\left( {{{\hat{y}}_i} - {y_i}} \right)}^2}}}{n}.} \end{equation*} \)

We have found Property Attributes and Location are more influential in predicting listing prices. This is consistent with the results of the multivariate effect size method, and the Person correlation is 0.9999.

Again, if we use the standardized LRT statistics to compare the influence level of factors, we will conclude Property Attributes is the most influential factor in explaining listing prices. However, shown by the RMSE method, Property Attributes and Location are similarly influential in predicting listing prices, which is consistent with results generated by the transformed multivariate effect size method. This underlines the usefulness of the quantile-matching transformation technique in the large sample context.