2.1 The effect of unbalanced sample sizes on r_pb

The derivation and limitations of r_pb are reviewed by McGrath and Meyer [3]. Two sets, $$ {\mathbf{y}}_{\mathrm{A}}^{}\in {\mathbb{R}}^{N_{\mathrm{A}}}$$ and $$ {\mathbf{y}}_{\mathrm{B}}^{}\in {\mathbb{R}}^{N_{\mathrm{B}}}$$, are combined to form a set of paired values, $$ \left\{(c_i,y_i)\right|c_i=\text{'}A\text{'}\vee \text{'}B\text{'},y_i\in {\mathbb{R}}^1,1\le i\le N,N=N_{\mathrm{A}}^{}+N_{\mathrm{B}}^{}\}$$, where c_i is a group membership label, and the {(c_i, y_i)} data correspond to the vectors, (c, y). The standard practice is to invoke a numeric {0, 1} representation for c to obtain an indicator vector, $$ {\mathbf{I}}_{\mathrm{c}}\in {\mathbb{R}}^N$$. Then, application of the Pearson product-moment formula produces the point-biserial correlation coefficient [3]

Fig 1

Quadratic dependence of the point-biserial correlation coefficient, r_pb.

For the fixed value r_pb = 0.2, there is a range for Cohen’s d and the sample size proportion, p_A. This ambiguity complicates the interpretation of r_pb as an effect size measure.

Fig 2

Nonoverlap proportion and point-biserial correlation.

Theoretical curves and estimated values for point-biserial correlation, r_pb, nonoverlap proportion, ρ_pb, and sample size adjusted correlation, r_pbd, for simulated data with unequal sample sizes (N_A : N_B = 15000 : 500) and the difference between mean values, $$ {\overline{y}}_{\text{A}}-{\overline{y}}_{\text{B}}$$. Compared to r_pbd, r_pb is attenuated due to the confounding effect of the binomial sampling factor. A: Uniform unit width $$ (\sigma =1/\sqrt{12})$$ distributions. B: Standard normal (σ = 1) distributions.

2.2 Statistical parameters for point-biserial variation

In this section, we consider the question of how to generate a set of parameters for statistical variation in point-biserial data. The fact that r_pb is subject to confounding effects suggests that replacing categorical labels with {0, 1} numeric values is an improper procedure, because the labels acquire arithmetic properties in an ad-hoc way. Instead, we propose a new framework where sort is used as an intrinsic property of both numbers and labels. Suppose there is a machine which generates numbers with labels, (c_i, y_i), in no particular order, placing them in a data table to produce a point-biserial data set. Then, the table can be sorted using either c or y, to obtain orderings denoted as y_c and c_y, respectively. As we discuss next, these orderings are associated with statistical parameters, v_c and v_y, respectively. However, there is no rule that specifies which parameterization, v_c or v_y, might be preferred. Therefore, we make the following proposition,

Proposition 1. Point-biserial variation is parameterized by the Cartesian product of statistical parameters for the y_c and c_y orderings,

2.3 Coordinates for a two-component system of distributed effects

In this section, we discuss the fact that d and ρ_pb are only two elements of a minimal set of parameters for representing point-biserial variation. The one-to-one correspondence, d ↔ r_pbd, will be discussed in section 2.4. Algebraically, v_c corresponds to the Cartesian product of statistical parameters for two sets of $$ {\mathbb{R}}^1$$ data, $$ {\mathbf{\text{v}}}_{\mathbf{\text{c}}}=({\overline{y}}_{\text{A}},S_{\text{A}}^2,N_{\text{A}})\times ({\overline{y}}_{\text{B}},S_{\text{B}}^2,N_{\text{B}})$$. Introducing the center of mass parameter, $$ \mu =({\overline{y}}_{\mathrm{A}}^{}+{\overline{y}}_{\mathrm{B}}^{})/2$$, the mean values vector is expressed as

The term ‘substantive significance’ has been used to refer to the magnitude of an effect that would be regarded as practically important in a given application [6]. Suppose functional or engineering requirements are expressed in terms of a vector, h, of system parameters. Then, the utility of an effect would be specified as a mapping, $$ u:\mathbf{h}\mapsto {\mathbb{R}}^1$$. The specification of u(h) would account for differences in cost-benefit trade-offs for variation in the {h_i} components. The substantive significance for the effect size would be determined by the mapping, u(h) → u(v_pb). Without this information, it is difficult to reach a consensus on the merits of an effect size. This explains the criticism of Cohen’s thresholds for small, medium, and large effects as “somewhat arbitrary” [16] and suggestions that the significance of the magnitude of an effect size depends on the research question [3, 17, 18].

A fundamental limitation arises from the fact that the (δ, μ) center of mass decomposition does not extend to higher dimensions in a straightforward way. Consider the group means vector for three sets, i.e., $$ ({\overline{y}}_{\mathrm{A}}^{},{\overline{y}}_{\mathrm{B}}^{},{\overline{y}}_{\mathrm{C}}^{})$$. The default center of mass parameter is defined as $$ \mu =({\overline{y}}_{\mathrm{A}}^{}+{\overline{y}}_{\mathrm{B}}^{}+{\overline{y}}_{\mathrm{C}}^{})/3$$. However, there is no standard procedure for choosing the two additional deviation parameters needed to specify a complete basis. Consequently, the formulation of an effect size measure for multiple group variation is not a well-posed problem, i.e., there is no unique solution [19]. This explains why Cohen’s d does not generalize to schemes involving more than two groups [20] and provides support for previous recommendations to break down ‘complicated hypotheses’, p. 526 [21], and ‘reduce any multiple-level or multiple-variable relationship’ into a set of two-variable effect size relationships [17]. This provides the raison d'être for the development of exploratory methodologies such as CART in high-dimensional data analytics [22, 23].

2.4 Homogeneous coordinates for Pearson correlation

In the effect size literature, it is accepted practice to distinguish three different types of effect size measure, ‘relationship’, ‘group difference’, and ‘group overlap’ [3, 7]. In this section, we discuss the fact that this classification is misleading. We have already discussed the fact that Cohen’s d, r_pbd and ρ_pb all serve as measures of nonoverlap (section 2.2). Now, we point out that r_pbd and Cohen’s d are two sides of the same coin because relationship and group difference correspond to different coordinate systems for representing fractional variation. Such correspondences are quite useful in exploring statistical dependence in high-dimensional data. Consider a vector $$ (a,b)\in {\mathbb{R}}^2$$. Division by the y-component produces the ratio vector, $$ \{\alpha =(\alpha ,1)\in {\mathbb{P}}^1|\alpha =a/b,b\ne 0\}$$. Ratios can be distinguished by their representations as points in the projective line, $$ {\mathbb{P}}^1$$. However, normalization of a ratio vector by the Euclidean length, $$ \Vert \alpha \Vert =\sqrt{{\alpha }^2+1}$$, produces the unit vector $$ \widehat{\alpha }$$, which is a point in the positive half-circle $$ {\mathbb{S}}_{+}^1$$. Thus, a fractional quantity can be represented as a point in either $$ {\mathbb{P}}^1$$ or $$ {\mathbb{S}}_{+}^1$$. Algebraically, the $$ {\mathbb{P}}^1$$ and $$ {\mathbb{S}}_{+}^1$$ representations are related by linear fractional transformations. In the terminology of projective geometry, a ratio corresponds to a perspective function, P(u, t) = u/t, for vector u [15]. The scaling invariance property of α is represented by the equivalence relation

Fig 3

Projective spaces for the representation of point-biserial correlation.

The point-biserial correlation coefficient, r_pb, corresponds to the point $$ (r_{\text{pb}}^{},\sqrt{1-r_{\text{pb}}^2})$$ on the positive half-circle, $$ {\mathbb{S}}_{+}^1$$, and the point $$ (\sqrt{p_{\text{A}}{p}_{\text{B}}}d,1)$$ on the projective line, $$ {\mathbb{P}}^1$$. The homogeneous coordinates $$ (\sqrt{p_{\text{A}}{p}_{\text{B}}}d,1)t$$ for $$ t\in {\mathbb{R}}^1$$ correspond to points on the line through the origin. {p_A, p_B}: sample size proportions, d: Cohen’s d.

Table 1

Homogeneous coordinates for Pearson correlation.

Effect size	$$	$$	$$
Pearson correlation	$$	$$	$$
Point-biserial correlation	$$	$$	$$
rpbd	$$	$$	$$

The representations for Pearson product-moment correlation as homogeneous coordinates in $$ {\mathbb{R}}^2$$, the vector $$ (r,\sqrt{1-r^2})\in {\mathbb{S}}_{+}^1$$, and the vector $$ (r/\sqrt{1-r^2},1)\in {\mathbb{P}}^1$$. Corresponding representations for the point-biserial correlation, r_pb, and sample size adjusted correlation, r_pbd, are also listed. Cohen’s $$ d=({\overline{y}}_{\mathrm{A}}-{\overline{y}}_{\mathrm{B}})/S_p$$, $$ \{{\overline{y}}_{\mathrm{A}}^{},{\overline{y}}_{\mathrm{B}}^{}\}$$: mean values, S_p: pooled variance, {p_A, p_B}: sample size proportions for ‘A’ and ‘B’ data, $$ t\in {\mathbb{R}}^1$$.

2.5 Point-biserial variation in regression tree analysis

The CART association graph was introduced in Paper1 as a new method for analyzing statistical association in point-biserial data. In this section, we investigate the role of point-biserial variation in rCART, particularly the connection between IG_MSE and r_pb, and introduce the rCART graph as a new method for analyzing association for (x, y) data. The CART decision tree algorithm creates a decision tree by recursive partitioning of the association between response and independent variables [2, 14]. Each node of the tree corresponds to a binary partition of the range of an independent variable. In standard implementations, the partition parameters for a node are determined by maximizing the information gain (IG) for the response variable in an exhaustive search of associations over all independent variables. The rCART implementation is of particular interest because it involves the analysis of point-biserial variation. In each iteration, the set of statistics obtained for partitions of an independent variable constitutes a CART association graph [2]. For the partition value $$ x_j\in {\mathbb{R}}^1$$, the data for a node (V) are divided into two subsets, i.e., V_A = {(x_i, y_i)|x_i ≤ x_j} and V_B = {(x_i, y_i)|x_i > x_j}, from which data vectors {y_A, y_B} are obtained. Alternatively, if x_j is categorical, the subsets are specified using matching criteria V_A = {(x_i, y_i)|x_i = x_j} and V_B = {(x_i, y_i)|x_i ≠ x_j}. The standard rCART impurity measure is the mean square error for the response, $$ \text{MSE}\left(\mathbf{y}\right)={\sum }_i{(\overline{y}-y_i)}^2/N_V^{}$$, where N_V is the sample size and $$ \overline{y}$$ is the mean [14]. Then, IG is defined as the parent node impurity minus the weighted impurities for the subsets

3.1 rCART association graphs for NHC quality measures

NHC data of the fourth quarter of 2018 were retrieved for 20 quality measures (Q_i) for 15341 nursing homes; detailed descriptions of these continuous variables can be found on the NHC website [26]. A histogram of the nursing home occupancy is shown in Fig 4A. Since performance estimates for nursing homes with low occupancy would be less reliable, a minimum occupancy criterion of at least 50 ‘Average number of residents per day’ was applied to obtain a restricted data set of 11053 nursing homes for further analysis [27]. Pearson correlation coefficients, r(Q_i, Q_j), and association graphs were calculated for all pairs of quality measures, {(Q_i, Q_j)|i ≠ j}. On average, the information gain for the rCART partition is larger when the (Q_i, Q_j) variables are highly correlated (Fig 5A); the r(Q_i, Q_j) correlations are distributed with 95% less than 0.16 and a maximum of 0.65. The distribution for ‘Number of outpatient emergency department visits per 1000 long-stay resident days’ (‘Emergency visits’) versus ‘Number of hospitalizations per 1000 long-stay resident days’ (‘Hospitalizations’) with correlation r = 0.37 is skewed, with a long tail towards larger values (Fig 4B). rCART association graphs are shown for the ‘Hospitalizations’ response and ‘Emergency visits’ partition variables (Fig 6A and 6B), and for the reverse, i.e., ‘Emergency visits’ response and ‘Hospitalizations’ partition variables (Fig 6C and 6D). The high correlation between r_pb and $$ \sqrt{p_{\text{A}}{p}_{\text{B}}}$$ (r = 0.99) is typical and indicates that variation in the binomial sampling factor overrides the smaller variation in Cohen’s d (Eq 2). We also note that the graphs for r_pb and IG_MSE (not shown) are superimposable, as expected from Eq 20 and because the variation in S_p is small. Thus, r_pb and IG_MSE mainly correspond to the variation in sample size proportion. In general, we observe that the association curves for r_pbd and ρ_pb can be categorized as monotonically increasing or decreasing, or even U-shaped (concave up), depending on how the (Q_i, Q_j) data are distributed. Here, the U-shaped dependence of r_pbd correlates well with δ (r = 0.999) and contrasts sharply with the concave down variation for r_pb. Consequently, r_pb and r_pbd produce very different rCART partitions (Table 2). In Fig 6A, the r_pb partition for the split value, x_j = 0.8, produces subnodes with comparable sample sizes, N_A = 5742 and N_B = 4890 (Table 2). It is useful to view this partition from a statistical perspective. As a first approximation, we expect that the majority of nursing homes belong to a broad distribution for average performance. Then, the r_pb partition with a split value close to the median, 0.85, is analogous to splitting a normal distribution nearly in half, producing subsets with different mean ‘Emergency visits’ values {0.5, 1.4} that nevertheless correspond to entities with average performance. Thus, r_pb and IG_MSE produce rCART subsets that are not well distinguished from a functional perspective. In comparison, for r_pbd, there are two possible rCART partitions at either low (x_j = 0.3) or high (x_j = 2.5) split values. Each partition produces a large subset corresponding to a broad distribution for average performance and a much smaller subset for either above- or below-average performance. Thus, r_pbd produces more functionally relevant classifications.

Fig 4

Skewed distributions for NHC quality measures.

A. Histogram of ‘Average number of residents per day’ for 15341 nursing homes. B. Two-dimensional Gaussian kernel density estimate of the distribution of ‘Number of outpatient emergency department visits per 1000 long-stay resident days’ (‘Emergency visits’) versus ‘Number of hospitalizations per 1000 long-stay resident days’ (‘Hospitalizations’), with correlation r = 0.37.

Fig 5

The relation between r_pbd and ρ_pb in rCART.

These graphs display data obtained from association graphs for 380 pairs of quality measures, {(Q_i, Q_j)|i ≠ j}. A. r_pbd effect size for rCART split versus correlation r(Q_i, Q_j). On average, the largest information gain is obtained when the response and partition variables are highly correlated. B. Correlation r(r_pbd, ρ_pb) between effect size and r(Q_i, Q_j) for association graphs. There is good correlation between r_pbd and ρ_pb in many cases, but there are exceptions.

Fig 6

rCART association graphs for effect size.

A,B: ‘Hospitalizations’ response versus ‘Emergency visits’ partition variables, with correlation r(r_pbd, ρ_pb) = 0.93. C,D: ‘Emergency visits’ response versus ‘Hospitalizations’ partition variables, with correlation r(r_pbd, ρ_pb) = 0.49. Bar plot histograms are shown for ‘Emergency visits’ (B inset) and ‘Hospitalizations’ (D inset). r_pb: point-biserial correlation coefficient, {p_A, p_B}: sample size proportions, r_pbd: sample size corrected correlation coefficient, ρ_pb: nonoverlap proportion, (δ, μ): center of mass parameters $$ ({\overline{y}}_{\mathrm{A}}-{\overline{y}}_{\mathrm{B}},({\overline{y}}_{\mathrm{A}}+{\overline{y}}_{\mathrm{B}})/2)$$.

Table 2

rCART subnode parameters.

Response variable	Partition variable	Split value	Subnode A	Subnode B
Hospitalizations	Emergency visits	rpbd: 0.3	1.8, 0.7, 9909	1.2, 0.7, 723
”	”	rpb: 0.8	2.0, 0.7, 5742	1.5, 0.7, 4890
”	”	rpbd: 2.5	2.5, 0.8, 320	1.7, 0.7, 10312
Emergency visits	Hospitalizations	rpbd: 0.7	1.0, 0.6, 10137	0.5, 0.4, 495
”	”	rpb: 1.7	1.2, 0.7, 5318	0.8, 0.5, 5314
”	”	rpbd: 3.3	1.5, 1.0, 330	1.0, 0.6, 10302

Summary of the rCART partition values and subnode statistics $$ (\overline{y},\sigma ,N)$$ for the association graphs in Fig 6. r_pb: point-biserial correlation coefficient, r_pbd: sample size corrected correlation, $$ (\overline{y},\sigma ,N)$$: mean value, standard deviation, sample size.

The importance of accounting for variation in both degrees of freedom, (r_pbd, μ), is illustrated in Fig 6B and 6D. Here, μ is monotonically increasing, and one of the r_pbd partitions might be preferred depending on μ. However, this requires an assessment of the cost-benefit trade-offs for (r_pbd, μ) variation, which will depend on the particular application. A close correspondence between r_pbd and ρ_pb is observed in many cases, with r(r_pbd, ρ_pb) ≥ 0.8 in 68% of the association graphs (Fig 5B), but there are many cases where they differ depending on how the (Q_i, Q_j) data are skewed. Fig 6C shows an example of the difference between the ρ_pb and r_pbd curves with r(r_pbd, ρ_pb) = 0.49. The r_pbd partition for the lower split value might be preferred because it is associated with higher ρ_pb, depending on how the cost-benefit trade-off is assessed for (r_pbd, ρ_pb) variation. Consequently, three coordinates (r_pbd, μ, ρ_pb) are needed to distinguish different forms of point-biserial variation. These observations provide support for previous remarks stating that interpreting the magnitude of an effect size as a measure of substantive significance depends on the particular application [6, 18]. A more precise approach would take into account the multidimensional nature of point-biserial variation and involve the specification of functional or engineering requirements for a relevant vector basis. Then, an analysis of the effect size for the system response could involve separate thresholds for each coordinate. The ability to account for all relevant degrees of freedom is also important in assessing reproducibility. A one-parameter representation using an effect size such as r_pbd or Cohen’s d gives an incomplete picture and leads to ambiguous results because of the loss of information.

3.2 Distributed effects in point-biserial variation

The reproducibility of nursing home performance data depends on stochastic effects in the measurement of patient outcome. Then, the observed data are associated with a distribution of data sets, $$ \mathcal{P}\left(\mathbf{y}\right)$$, and corresponding distributions of the statistical parameters $$ \mathcal{P}\left({\mathbf{v}}_{\text{pb}}^{}\right)$$ and effect size. The specification of $$ \mathcal{P}\left(\mathbf{y}\right)$$ must be based on a realistic assessment of all sources of error and uncertainty to form an error model for the data, $$ \mathcal{E}\left(\mathbf{y}\right)$$. Then, the determination of the distribution for the effect size requires propagation of the error in $$ \mathcal{P}\left(\mathbf{y}\right)$$. For fractional quantities such as Cohen’s d and r_pbd, it is necessary to account for stochastic effects in both the numerator and denominator. However, analytical methods for estimating distributions for ratios [28, 29], proportions [30, 31], and correlation coefficients [32] are complicated by fractional transformation, a bounded range, and discreteness. Thus, iterative procedures are needed for the analysis of noncentral effect size distributions and estimating confidence intervals for deviations above and below the effect size estimate [5, 18]. Alternatively, Monte Carlo (MC) methods [2, 33, 34] provide a more practical approach to estimating the distribution for the effect size. In an MC simulation, $$ \mathcal{E}\left(\mathbf{y}\right)$$ specifies error parameters for each observed value in the original data. Then, a point-biserial MC data set is obtained by random sampling to produce MC instances for y_A and y_B. The MC sampling process is repeated many times to obtain a collection of MC data sets to form an estimate, $$ {\mathcal{P}}_i\left(\mathbf{y}\right)$$. Statistical parameters are calculated for the data sets in $$ {\mathcal{P}}_i\left(\mathbf{y}\right)$$ to obtain estimates of distributions and histograms for point-biserial effects. Many MC runs are performed to obtain a set, $$ \left\{{\mathcal{P}}_i\left(\mathbf{y}\right)\right|1\le i\le N_{\text{MCruns}}^{}\}$$, which allows the determination of the degree of convergence for the MC simulation. However, the information needed to construct an error model is not included in the NHC quality measures data. For this demonstration, we provided a rudimentary ‘Emergency visits’ error model, where σ_i = y_i/5. MC simulations for (r_pbd, μ) and (r_pbd, ρ_pb) for ‘Emergency visits’ response with ‘Hospitalizations’ rCART split value, 3.3 (Table 2), are shown in Fig 7. The discrete structure of the ρ_pb distribution is due to stochastic effects in the c_y sorting. The separate confidence intervals in Fig 6 for positive and negative deviation from the observed effect size estimate were estimated from the MC distributions. In practical applications, the advantage of the MC method is that it allows detailed simulation of the data acquisition process, including heterogeneity within groups, and specifications for $$ \mathcal{E}\left(\mathbf{y}\right)$$ can include heteroscedasticity, measurement error, and misclassification [17, 35, 36].

Fig 7

Monte Carlo simulation of the distribution of stochastic effects for point-biserial variation.

2D histograms of MC distributions for (r_pbd, μ) (A) and (r_pbd, ρ_pb) (B) for ‘Emergency visits’ response with ‘Hospitalizations’ rCART split value, 3.3 (Table 2). The 1σ error bars for the r_pbd histogram (A inset) serve as an indication of convergence for the simulation; the mean for the normal curve corresponds to the observed r_pbd value, 0.398. r_pbd: sample size corrected correlation, ρ_pb: nonoverlap proportion, μ: center of mass parameter $$ ({\overline{y}}_{\mathrm{A}}+{\overline{y}}_{\mathrm{A}})/2$$, number of MC runs: 25, samples per MC run: 4000.