Analyzing admissions metrics as predictors of graduate GPA and whether graduate GPA mediates PhD completion

Student level data for both this study and [1] was requested from physics departments that awarded more than 10 PhD’s per year for students who matriculated between 2000 and 2010, including information on undergraduate GPA (UGPA), GRE-Q, GRE-V, GRE-P, and graduate GPA (GGPA). Data collected also included the final disposition of students (PhD earned or not), start and finish years, and demographic information. GPA data is analyzed on a 4.0 scale while GRE scores are on the percentile scale.

We received data from 27 programs (approximately a 42% response rate), which spanned a broad range of National Research Council (NRC) rankings. The sample used in [1] consisted of all students in 21 programs for which start year was available. Given that the median time to degree across physics PhD programs is 6 years, some students who started before 2010 were still active at the time of data collection in 2016. The probability of not completing the physics PhD has an exponential time dependence with a time constant of 1.8 years. Thus, students who have been in their programs for three time constants have only a 5% chance of not completing. These students were thus categorized as completers in this study.

These data covered 3962 students. Of this subset, two programs did not report GGPA data for their students. Hence, the sample for this study excludes these students, thereby reducing the sample size to 3406 students across 19 programs. This corresponds to approximately 11% of matriculants to all U.S. physics PhD programs during the years studied.

Among the sample of US students, 16% are women ($N=317$). Although the authors generally advocate for a nuanced treatment of gender in physics education research and recognize the deficits associated with treating gender as a fixed binary variable [42], the present data set spans the years 2000 to 2010 during which the data collected by programs only allowed for the binary option of male/female. Hence, we must treat gender as a dichotomous variable in this analysis.

The racial composition of the dataset is 61.6% White, 1.3% Black, 2.1% Hispanic, 0.2% Native American, 3.5% Asian, 1.0% multiple or other races, and 30.2% undisclosed. Excluding the cases for which race was unavailable, the sample is thus roughly representative of annual PhD production in U.S. physics for gender, race, and citizenship [43]. We include race as a covariate in each analysis presented; however small $N$, particularly for Black, Hispanic, Native American, and Asian students, often precludes useful interpretation of the results pertaining to these subsets.

In order to focus on issues of diversity and inclusion associated most strongly with US applicants, we use only the subset of data from domestic graduate students. This decision is further motivated by research suggesting that it is difficult for admission committees to directly compare scores earned by US and international students, indicating that separate analyses are appropriate [9]. Using the subset of students who are from the US reduces the total sample size for the study to $N=1955$. A cursory visualization of the variables in the data set, as shown in Figure 1, shows that the distributions of scores for US and Non-US students are markedly different, which further justifies separate analyses of these two student populations.

Examining the distributions of scores in Figure 1, the presence of non-normality is evident in nearly all of the variables. Each of the continuous variables in the dataset fail the Shapiro-Wilk test of normality at the $\alpha$ level of .05. However, these tests are often of limited usefulness; in general distributions with skewness $\lvert\hat{\gamma_{1}}\rvert>3$ or kurtosis $\lvert\hat{\gamma_{2}}\rvert>10$ likely indicate that they violate any assumption of normality [44]. For this dataset, the GGPA distribution skewness $\hat{\gamma_{1}}=-3.31$ and kurtosis $\hat{\gamma_{2}}=20.43$, indicating severe non-normality. The GRE-Q distribution also falls into the problematic range ($\hat{\gamma_{1}}=-2.11$ and $\hat{\gamma_{2}=9.33}$). Ceiling effects are also present, since many students earned 4.0 grade point averages or earned the maximum score on the GRE examinations.

The data collection process was limited to gathering only cumulative graduate GPA rather than first-year graduate grades, which were not recorded by some programs. Thus, depending on whether a student persisted in a program, their graduate GPA may be based on many courses while others are based on only a few courses. The data set is also necessarily subject to range restriction, since data on student performance in graduate school is automatically limited to include only students who were accepted to undertake graduate study. We cannot know how students who were not accepted into graduate school would have performed had they been accepted. Range restriction may act to attenuate the strength of observed effects in subsequent analyses [45].

Although not used in a majority of this study, we briefly explore the role of the doctoral programs’ NRC ranking in PhD completion [46]. Since the NRC only gives confidence intervals for program rank, we created a ranking for this study by averaging the 5 and 95% confidence bounds for the NRC regression-based ranking (NRC-R) and rounded this up to the nearest five to protect the confidentiality of participating programs. This led to a ranking range of 5 to 105. We divided the programs into terciles of approximately equal number of records, and categorized as Tier 1 (highest ranked, NRC-R $\leq$ 20), Tier 2 (25 $\leq$ NRC-R $\leq$ 55), and Tier 3 (NRC-R $>$ 55).

Multiple imputation (predictive mean matching) is used to impute missing UGPA and GRE-P scores. Predictive mean matching is used due to the non-normality of the data. 160 students do not have data for either UGPA or GRE-P, while 400 are missing only UGPA and 263 are missing only GRE-P. All multiple imputation is conducted using the mice package in R [47]. 20 imputed data sets are used for each analysis. For consistency, incomplete variables are imputed using the same imputation model used in [1], in which the imputation model utilizes all other variables in the data set aside from graduate GPA and PhD completion. The model utilizes GRE-Q, GRE-V, program tier, gender, and race, as well as complete cases of UGPA and GRE-P. Although the imputation approach presented here is theoretically sound, we also present a comparison of several different models of data imputation in the Supplemental Materials.

The goal of this section is twofold. First, we seek to gain a cursory look at how graduate GPA is related to common admissions metrics. In doing this, we also wish to determine whether it is reasonable that admissions metrics could indirectly predict completion through graduate GPA. This section presents a series of analyses meant to elucidate the relationships between standard admissions metrics, students’ GGPA, and students’ final disposition. To make our analysis maximally accessible to readers of different statistical backgrounds, we describe here in detail the methods used in this section.

Bivariate correlation coefficients provide information about the level of association between two variables, and are therefore a useful starting point for analysis. We construct a correlation matrix (see Table 2) for all variables in the sample using Pearson correlation coefficients, which are equivalent to the standardized slope coefficients for a linear model predicting $y$ from $x$. These are given by

for any two continuous variables $x$ and $y$. Calculating $r_{xy}$ gives us a first glance at the relationships between the continuous variables UGPA, GGPA, GRE-P, GRE-Q, and GRE-V.

When the $x$ variable is treated as dichotomous (e.g., gender in this data set), Eq. (1) reduces to the point-biserial correlation coefficient $r_{pb}$,

where $\bar{y_{1}}$ and $\bar{y_{0}}$ are the means of the continuous $y$ variable for the two $x$ groups 1 and 0, $q$ and $p$ are the proportions of data belonging to these two groups, and $\sigma_{y}$ is the standard deviation for the $y$ variable. Like the Pearson coefficient, the quantity $r_{pb}$ ranges from -1 to 1 and indicates the strength of association between two variables. Conveniently, a significance test for the point-biserial correlation is identical to performing an independent t-test on the data [48]. Thus, the point-biserial correlation coefficient yields information about whether two group means are statistically different. For instance, the point-biserial correlation tests whether the GGPA of male students are statistically different from those of female students (we find that GGPAs are not significantly different by gender, see Table 2).

When $x$ and $y$ are both dichotomous, the Pearson coefficient reduces to the phi coefficient,

where $\chi^{2}$ is the chi-squared statistic for a 2x2 contingency table and $n$ is the total number of observations in the data. The phi coefficient also ranges from -1 to 1 and indicates the strength of association between two binary variables. This quantity allows us to examine whether final disposition is significantly associated with gender (we find that the association between gender and final disposition just meets the threshold for statistical significance ($\phi=-0.05\pm 0.04$, $p=0.04$), likely due to the large sample size of our data, but the very small phi coefficient indicates that the practical strength of this relationship is negligible [49]).

To characterize how GGPA and other numerical predictors vary across program tier, we conduct several one-way analysis of variance (ANOVA) tests using program tier as the independent variable. ANOVA tests allow us to determine whether there are significant differences between different groups, such as students in different program tiers. These tests produce an F-statistic, which is interpreted as the ratio of between-group variability to within-group variability. Thus, higher values of F indicate that between-group variability is large compared to within-group variability, which is unlikely if the group means all have a similar value.

Lastly we present the results of a multiple regression analysis in which we regress GGPA on common admissions metrics and demographic factors. Regression allows us to examine the unique predictive effects of these predictors.

The classical linear regression model is written mathematically for an outcome variable $Y$ as

where $i=1,...,n$, the number of observations in the data, and $k$ represents the number of predictors in the model. Error terms $\epsilon_{i}$ are assumed to be independent and normally distributed with mean 0 and standard deviation $\sigma$. $\hat{\alpha}$ is the vector of regression coefficients that minimizes the sum of squared errors

for the given data. The regression coefficients can be interpreted as the difference in the outcome variable $Y$, on average, when comparing two groups of units that differ by 1 in one predictor $X$ while keeping all the other predictors the same.

We report both unstandardized and standardized versions of the regression coefficients. Unstandardized coefficients are the result of regression analysis using the original, unscaled variables. Thus, the unstandardized regression coefficients represent the predicted average change in the outcome variable $Y$ when the corresponding predictor $X$ is changed by one unit. This allows for a straightforward interpretation since the variables are not scaled, but does not yield insight into the relative predictive strengths of the independent variables since they are scaled differently. Standardized regression coefficients result from regression analyses using continuous variables that have been mean-centered and divided by their standard deviation, resulting in variables with variances equal to 1. Thus standardized regression coefficients represent the average number of standard deviations changed in the outcome variable when a predictor variable is increased by one standard deviation. By calculating the standardized coefficients, we exchange a simple interpretation of score change for an interpretation of which variables have the greatest effect on the dependent variable.

Using mediation analysis we seek to answer the question of whether graduate GPA mediates the predictive ability of common admissions metrics on PhD completion. Whereas analyses such as logistic regression [1] yield information about whether independent variables such as UGPA and GRE-P affect final disposition of a graduate student, they do not offer insight into the explanation of why and how UGPA, GRE-P, and other admissions metrics affect completion. Mediation analysis is one technique that allows us to probe the underlying process by which some variables influence others [19, 20, 21, 22, 23].

Figure 2 graphically depicts a prototypical mediation model, where $X$, $Y$ and $M$ represent the model’s independent, dependent, and mediating variables. $C$ represents a covariate. As a hypothetical example, let’s say that previous research has shown a positive relationship between the use of active learning activities in physics class and student exam grades. Researchers might hypothesize that this relationship is actually due to a third mediating variable, student engagement. Using active learning activities in class may cause students to become more engaged with the material, making their subsequent exam grades increase. Engagement is a mediating variable in this case. Meanwhile, since students who are already interested in physics could be predisposed to being more engaged and performing better on exams, the researcher might take students’ prior physics interest into account by including it as a covariate in their analyses.

In this section, we wish to discern whether graduate GPA mediates the relationship between common admissions metrics and students’ likelihood of completing graduate school. In practice, mediation analysis is done by simultaneously estimating a set of regression equations [50]. The goal is to partition the total effect of the independent variable $X$ on the dependent variable $Y$ into two parts: the direct effect of $X$ on $Y$ and the indirect effect of $X$ on $Y$ through the mediating variable $M$. For the simple example given above, the set of regression equations to be solved are:

Traditional mediation literature [51, 52] defines the direct effect of $X$ on $Y$ as the coefficient $\beta_{2}$ and the indirect effect of $X$ on $Y$ as the product of coefficients $\gamma_{1}\beta_{1}$, corresponding to the products of the path coefficients along the mediated path shown in Figure 2. Statistically significant values of $\gamma_{1}\beta_{1}$ indicate that the relationship between $X$ and $Y$ is mediated by $M$. This method demonstrates the general intuitive ideas underlying mediation analysis, but is subject to several important limitations. Foremost among the limitations associated with traditional mediation analysis is that its applicability to model categorical variables (e.g., binary outcomes) and nonlinearities is not well defined, as these situations preclude the use of sums and products of coefficients [53, 19, 54]. The difficulties associated with binary outcomes are therefore problematic for a model predicting final disposition, a binary outcome. Furthermore, traditional mediation models leave the causal interpretation of their results ambiguous [55].

Recent work in the field of causal inference [16, 17, 18] has formalized and generalized mediation analysis to resolve these limitations, allowing for categorical outcomes while also clarifying that under certain conditions the results may be interpreted causally. In this framework, often called the “potential outcomes” framework, the traditional product-of-coefficients mediation analysis is a special case for which the mediator and outcome variables are both continuous, while the functional forms of the direct and indirect effects for other situations become more complicated [50].

For the primary analysis of this paper, we calculate the direct and indirect effects defined by the potential outcomes framework for the case of a continuous independent variable (UGPA, GRE-P, and GRE-Q), a continuous mediator (GGPA), and a dichotomous outcome (final disposition). The simultaneous regression equations to be calculated are still Eqs. (6) and (7), except the binary $Y$ is replaced with $Y^{*}$, a continuous unobserved latent variable which represents the observed binary variable. Once estimated, the direct and indirect effects reduce to simple differences in probability of completing a PhD between students across different values of the independent variables. Mathematically, the effects for a change in the independent variable from a value $x_{0}$ to $x_{1}$ at a particular value of the control $c$ are given by

where $\Phi$ represents the normal cumulative distribution function and $\text{probit}(x_{a},x_{b})$ is given by

and $v(x_{a})$ is

Note that these expressions are all still simply combinations of the coefficients from the regression equations (6) and (7). Thus, the potential outcomes framework allows us to calculate the total predicted change in probability of completing a PhD due to an independent variable and decompose it into that variable’s indirect effect on final disposition through GGPA as well as its direct effect (see Figure 4).

Using this mediation framework can help to give powerful insights into nuanced relationships between the variables in an observational study [56]. However, giving a truly causal interpretation to the results of this mediation analysis requires a set of strong assumptions to be met, and in practice it can be difficult for any observational study fully meet these conditions [57]. Hence, we try to avoid making explicitly causal claims in our discussion of the results. We discuss the assumptions needed for a causal interpretation as well as the robustness of the current study to violations of those assumptions in Section V.2.

Mediation analyses were conducted using the mediation package in R [58]. Checks for the consistency of results across different computational approaches were done by performing duplicate mediation analyses in the R package medflex [59] as well as the statistical software Mplus [60].

Correlations An initial question related to predicting a student’s final disposition is whether UGPA, GRE scores, and GGPA are reliably correlated with a student’s final outcome. We are also interested in the strength of association between GGPA, UGPA, and GRE scores, as this information yields insight into whether GGPA could serve as a mediating variable in predicting final disposition. Table 2 contains the bivariate correlations (Pearson, point-biserial, and phi) between each pair of measures for the sample in the lower diagonal. The 95% confidence intervals are reported in the upper diagonal. Confidence intervals that do not include a value of zero indicate that the correlation is statistically significant. The means and standard deviations of the continuous variables are presented in the table’s first column (GPA data is analyzed on a 4.0 scale while GRE scores are on the percentile scale).

Inspection of Table 2 reveals that GGPA is the predictor most strongly correlated with final disposition ($r_{pb}=0.43$). This value is statistically significant ($p<.001$) and positive, meaning that students with higher GGPA are more likely to finish their PhD program successfully. This trend is visually apparent in Figure 3(a), which shows boxplots of GGPA grouped by PhD completers and non-completers. We also observe that UGPA ($r_{xy}=0.29$, $p<.001$), GRE-Q ($r_{xy}=0.18$, $p<.001$), and GRE-P ($r_{xy}=0.22$, $p<.001$) are all positively correlated with GGPA, albeit weakly, meaning that students with higher scores in these metrics tend to earn higher GGPAs. Taken together, the observation that higher UGPA and GRE scores positively correlate to GGPA, which in turn correlate with a student’s likelihood of completion, implies that GGPA might play an important role in mediating the influence of these admissions metrics on PhD completion.

The lack of a statistically significant correlation between gender and GGPA indicates that the disparity in scores on the GRE-P between males and females does not manifest itself in subsequent GGPA performance. Indeed, there is no statistical difference between average GGPA for males and females ($r_{pb}=-0.01$, $p=0.57$, equivalent to a non-significant independent $t$-test), as demonstrated in Figure 3(b). Yet there exists a statistically significant difference between males and females in GRE-P performance in our data ($r_{pb}=-0.30$, $p<.001$, equivalent to a statistically significant independent $t$-test). Thus GGPA does not differ between genders despite the known performance gap between males and females on the GRE-P exam. Furthermore, the phi coefficient measuring the association between gender and PhD completion is negligible despite barely meeting the threshold of statistical significance ($\phi=-0.05$, $p=0.04$).

Still, the bivariate correlations shown in Table 2 do not control for possible relationships between the variables of interest. For instance, the moderate correlation between GRE-Q and GRE-P ($r_{xy}=0.51$, $p<.001$) indicates that there may be a spurious relationship between one of these variables and GGPA. In addition, we observe low but statistically significant correlations between UGPA and GRE-P ($r_{xy}=0.31$, $p<.001$), GRE-Q ($r_{xy}=0.30$, $p<.001$), and GRE-V ($r_{xy}=0.17$, $p<.001$). This is expected, as UGPA likely contains some information regarding the specific aspects of students’ aptitudes tested by GRE exams. These results motivate the use of multiple regression analysis later in this section in order to disentangle the unique effects of each independent variable on GGPA. That analysis reveals that when we isolate the unique predictive effects of each variable in the regression model, UGPA and GRE-P remain significant but weak predictors of GGPA, while GRE-Q does not retain statistical significance.

Similarly, although UGPA ($r_{pb}=0.15$, $p<.001$), GRE-P ($r_{pb}=0.14$, $p<.001$) and GRE-Q ($r_{pb}=0.15$, $p<.001$) are positively correlated with PhD completion, the magnitudes of these correlations are very weak and do not account for other parameters that may be associated with completion. Multivariate approaches allows us to isolate how individual metrics relate to PhD completion, which we explore in the mediation analysis presented in Section IV.2. Consistent with previous studies of PhD completion using multivariate approaches [1], we find that when accounting for other parameters, only UGPA remains a statistically significant predictor of completion.

Results of one-way independent ANOVA tests show that the main effect of program tier on GGPA is significant, $F(2,1949)=26.31$, $p<.001$, which reflects the upward trend in GGPA from Tier 3 to Tier 1 and 2 programs. A Tukey post hoc test reveals that the GGPA was significantly higher for students at Tier 2 ($M=3.59$, $SD=0.40$, $p<.001$) and Tier 1 ($M=3.60$, $SD=0.42$, $p<.001$) institutions than those at Tier 3 institutions ($M=3.41$, $SD=0.59$). There was no statistically significant difference between the Tier 1 and Tier 2 groups ($p=0.97$).

Multiple Regression To disentangle the unique effects of each predictor on GGPA we conduct a multiple linear regression analysis. Multiple linear regression allows us to simultaneously fit many independent variables to measure each of their relative effects on a single dependent variable, GGPA. Analyzing the raw coefficients fitted by the regression analysis yields insight into the predicted change in GGPA due to changes in one variable while holding all others constant. Standardized coefficients allow for a comparison of the relative effect sizes of the independent variables.

Our model includes all available GRE scores (GRE-P, GRE-Q, and GRE-V) as well as UGPA in order to examine the unique predictive effects of each measure. We considered the possibility that including both GRE-P and GRE-Q in the same model would raise collinearity concerns, but find these concerns unfounded. The bivariate correlation ($r_{xy}=0.51$) between the two is not high enough to warrant genuine concern [61, 62]; furthermore, the variance inflation factor (VIF) for every imputed dataset’s regression model was below 1.75, well below the commonly cited threshold of 10. Hence, we deem the model posed in the study as most appropriate to answer the research questions raised in this study (further discussion regarding collinearity concerns in the data is available in the Supplemental Material).