Just so we’re clear, this:
x \sim T(\mu, \sigma, \nu)
Is equivalent to this:
x \sim \mu + T(0, \sigma, \nu)
Including if \mu \sim T(\mu_0, \sigma_0, \nu_0)
My point is that when you add together year_mean and runner_mean, you end up with a unique mean for each runner-year combination. It is exactly equivalent to what you originally wrote, meanRunner = pm.StudenT('meanRunner', mu=meanYear, ...). The number of combinations of runner + year is governed by the indices – each unique addition of two means leads to a unique mean.
If there are multiple races within each year, this is no problem. You just need to introduce another index, race_id. The data becomes something like this:
times
runner year race_id
A 2018-01-01 0 -0.050662
1 0.921232
2 -0.665283
2019-01-01 0 -1.239528
1 -0.884322
2 0.423450
3 -1.070264
2020-01-01 0 2.510278
1 1.525603
2 1.790918
B 2018-01-01 0 1.141352
1 1.347654
2 4.210434
3 1.485402
2019-01-01 0 1.636708
You don’t need to factorize race_id though, just factorize runner and year. All the races of racer A in 2018 will share the same mean, all the races of racer A in 2019 will share a mean, and so on.