Exact vs. estimated population variances. I suppose you do not know the exact values of $\sigma_P^2$ and $\sigma_T^2,$ and that you would use respective sample variances $S_P^2$ and $S_T^2$ as estimates. Then technically speaking, you have no choice but to do a t test.
(If you do a z test with estimated variances, you are really doing an approximate t test, depending that the large degrees of freedom with provide
a good approximation.)
Equal vs. unequal population variances. I see no basis for assuming
that $\sigma_P^2 = \sigma_T^2,$ so your test should not assume equality.
If you do a two-sample t test, then you should use the Welch (not pooled)
version. (In R, the Welch version is the default: you would have to use
parameter var.eq=T
in the t.test
procedure to select the pooled version.)
If you do an approximate z test, then be sure to supply proper estimates
for $\sigma_P^2$ and $\sigma_T^2$ in the denominator of the test statistic.
Checking for normal data. With several thousand in each group (P and T),
the t test should be sufficiently robust against mild non-normality that you
would have no qualms about doing the t test. If individual boxplots of P and T
groups separately show only a few moderate outliers, and show no evidence of extreme skewness, that would indicate that using the Welch t test is OK. (With data values between 1 and 14, I would not expect difficulties.)
As a technical matter it is the residuals that need to be normal for a precisely accurate t test. So if it is important to check for normality, you should look at P and T groups separately. (However, formal normality tests, such as Shapiro-Wilk, may be far too aggressive in detecting harmless departures from normality. See Comment of @SalMangiafico and my Addendum.)
Summary. Without knowing whether your data can be considered as a valid random sample from the two groups, without seeing actual data data values, and just guessing from academic
experience about PhD finishing times, my guess is that you should use a
Welch two-sample t test.
Addendum. (Prompted by various comments.)
Suppose both P and T samples have means somewhat above 4 years, as for
the fake data below. I am not suggesting that these data are similar to yours, but they do illustrate uses of various possibly relevant methods.
set.seed(2019)
prm = rbinom(3000, 9, .35) + 1
tmp = rbinom(2700, 9, .4) + 1
summary(prm); sd(prm)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 3.000 4.000 4.133 5.000 9.000
[1] 1.423371
summary(tmp); sd(tmp)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 4.000 5.000 4.566 6.000 9.000
[1] 1.480406
boxplot(prm,tmp, col="skyblue2", horizontal=T, notch=T)
The 'notches' in the sides of the boxes show nonparametric confidence
intervals calibrated so that non-overlapping CIs for two groups
suggests a significant difference in locations at the 5% level.
My (essentially binomial) data in Group P look nearly normal in a histogram (with 'best-fitting' normal density). However, the Shapiro-Wilk test strongly rejects their normality---possibly, mostly just because all values
are integers. See Note. [Similar results for the data in Group T are not shown.]
shapiro.test(prm)$p.val
[1] 1.315475e-29
Both the two-sample Welch t test and the (nonparametric) two-sample Wilcoxon
signed-rank test strongly reject the null hypothesis that the two populations
the same. With such large samples, you have good power to detect
even fairly small differences in location between time until PhD for P and T groups. In the figure, the intervals are far from overlapping.
t.test(prm, tmp)
Welch Two Sample t-test
data: prm and tmp
t = -11.236, df = 5581.4, p-value < 2.2e-16 # df=5698 for 'pooled'
alternative hypothesis:
true difference in means is not equal to 0
95 percent confidence interval:
-0.5088930 -0.3576996
sample estimates:
mean of x mean of y
4.133000 4.566296
wilcox.test(prm, tmp, conf.int=T)
Wilcoxon rank sum test with continuity correction
data: prm and tmp
W = 3371100, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
-9.999792e-01 -3.957901e-05
sample estimates:
difference in location
-2.474458e-05
Note: In the experiment below, the Shapiro-Wilk test does not recognize normal data rounded to
integers as being normal.
z = rnorm(3000, 5, 1); shapiro.test(z)$p.val
[1] 0.583755 # Accept null: data normal
rz = round(z); shapiro.test(rz)$p.val
[1] 5.377162e-37 # Reject null: rounded data not normal