Yes, this line is correct:
Y_obs = pm.Normal(“Runner”, mu=mu,sigma=sig,observed=data.times)
To recover specific combinations of runner + times, I think the easiest thing will be to make the full matrix of runner-year pairs and save it as a deterministic, and give it named dims, like this:
all_mu = pm.Deterministic('all_mu', runner_mean[:, None] + year_mean[None], dims=['runner', 'year'])
I am broadcasting the two effects together by making runner_mean into a column vector, so the runners are on the rows, and year_mean into a row-vector, so years are on the columns. Thus, all_mu[i, j] will be the mean for the i-th runner in the j-th year. To plot the effects, you can then use the named dimensions. This plots the KDE for the posterior mean of runner A in year 2018:
az.plot_posterior(trace, var_names=['all_mu'], coords={'runner':['A'], 'year':[2018]})
This also gives you counter-factual means for the runner-year pairs that you didn’t observe in the data, which might or might not be interesting.