What is the factor structure of the Problem Gambling Severity Index (PGSI; Ferris & Wynne, 2001)? Since its creation, it has been assumed that there is only a single factor underlying this gold-standard index of disordered gambling symptomatology. However, recently, Tseng, Flack, Caudwell, and Stevens (2023) tested and indicated support for their hypothesis that the nine PGSI items reflect two factors—problem gambling behaviours (four items) and negative consequences of problem gambling behaviours (five items). Importantly, they used confirmatory factor analysis (CFA; Kline, 2016) for their analyses. A two-factor model provided a better fit to the data relative to a single-factor model. However, neither model provided a good fit to the data via the chi-square test of model fit (χ2) or the approximate fit indices, including the comparative fit index (CFI) and root mean square error of approximation (RMSEA). Fit statistics are reported in Table 1. Fit is considered excellent when the χ2 test is not statistically significant, CFI ≥ 0.95, the RMSEA ≤ 0.05 and its the 90 % confidence interval contains zero (Kline, 2016). Also, in the two-factor model, the two factors were highly correlated (r = 0.92), which indicates near-complete overlap between the scores on both factors. Taken together, these observations raise concerns about the internal validity of the results from the CFAs. Lastly, other plausible alternative models were not tested.

Herein, we argue that a hierarchical model may best fit the factor structure of the PGSI. A hierarchical model of the PGSI would include one global factor underlying all nine items that co-exists with the two sub factors (i.e., problem behaviours and negative consequences of problem gambling behaviours). In the hierarchical model, the global factor is not correlated with the sub factors. If support is found for a hierarchical model of the PGSI, then the two sub factors observed by Tseng et al. account for unexplained variance in the items that is not accounted for by the global factor. Put differently, the PGSI may be, as Ferris and Wynne (2001) observed, a single factor. However, there may also be two sub factors that account for method or residual variance due to item domain.

We contacted Tseng et al. to obtain a copy of the PGSI data used in their published paper. They promptly provided the polychoric correlation matrix they used as input in their CFAs. They also provided the raw data. For better accuracy, we used the raw data in our analyses. The fit of a single-factor model and a two-factor model were tested, then the hierarchical model was tested. The nine PGSI items were modelled as ordered categorical because of their severe positive skew (range: 5.597 to 17.954) and kurtosis (range: 30.013 to 389.234) using the robust weighted least squares estimator in Mplus version 8.2 (Muthén, & Muthén, 1998–2017). Models were compared using a chi-square difference test (Δχ2). The w effect size was used to calculate the degree of increase in relative fit between the competing models (Newsom, 2015). A value of 0.10 is small, 0.30 is medium, and 0.50 or more is large.

Fit statistics are reported in Table 1. The χ2 values for the single-factor and two-factor models were statistically significant indicating poor fit, but the fit indices suggested approximate fit. There were no residuals > |0.10| in both models. In the two-factor model, the two factors were strongly correlated, r = 0.97, p < 0.0001. Although the two-factor model provided a better fit to the data compared to the single-factor model, Δχ2(1) = 6.29, p = 0.01, the relative increase in fit was very small, w = 0.04.

As expected, the hierarchical model provided an excellent fit to the data via χ2 and the approximate fit indices (see Table 1). There were no residuals > |0.10| and the two sub factors were correlated, r = 0.69, p < 0.001. Critically, the hierarchical model provided a better fit to the data relative to the two-factor model, Δχ2(9) = 42.79, p < 0.001. The magnitude of the increase in fit due to adding a global factor was medium in size, w = 0.32. Also, the hierarchical model provided a better fit to the data relative to the single-factor model, Δχ2(10) = 48.84p < 0.001. The magnitude of the increase in fit due to adding the two sub factors was medium in size, w = 0.36. Together, these observations indicate that a hierarchical model provided the best fit to the data.

The standardized factor loadings from the hierarchical model are reported in Table 2. Of note, all items loaded strongly on the global factor. Internal consistency reliability analyses using coefficient omega for ordered-categorical items (see Green & Yang, 2009) showed that the PGSI total-score had high reliability, 0.87. Thus, the global factor was empirically well-defined.

As for the two sub factors, two of the four items on the behaviours factor and two of the five items on the consequences factor had loadings that were not statistically significant. Also, among loadings that were statistically significant, all were < 0.50 on the behaviours factor and only one was > 0.50 on the consequences factor. Lastly, the reliable portion of the variance in the behaviours and consequences subscales scores that is independent from the global factor were 0.04 and 0.07, respectively. Thus, although the inclusion of the two sub factors in the hierarchical model enhanced model fit relative to the single-factor model, the two subfactors were not empirically well-defined devoid of the global factor. The unique variance in each sub factor independent of the global factor may reflect method variance due to item domain.

All input and output files of the statistical analyses reported herein are publicly available via the Open Science Framework: https://osf.io/sa9mk/.