# Is this an appropriate way of modeling the scores in a round-robin sports league?

I want to model the outcome of matches in a round-robin sports league based on which home team is playing which away team across several seasons.

Let's assume a league with four teams A, B, C, and D. In each league season, the four teams play each other two times with alternating home teams. So, for each season, there are 12 matches (home team A vs. away teams B, C, D; home team B vs. away teams A, C, D; home team C vs. away teams A, B, D; home team D vs. away teams A, B, C). Let's also assume that the outcome of each match is an integer score based on a normal distribution with m=0 and sd=5.

The following R code produces a random test data set for five seasons:

set.seed(123)
dat <- data.frame(
D=rep(1:12, 5),
Home=rep(c("A", "A", "A", "B", "B", "B", "C", "C", "C", "D", "D", "D"), 5),
Away=rep(c("B", "C", "D", "A", "C", "D", "A", "B", "D", "A", "B", "C"), 5),
Season=factor(rep(2014:2018, each=12)),
Score=round(rnorm(5*12, 0, 5)))


Now, I'm struggling with the appropriate model specification for fitting Score. The score will depend on which team plays which team, so Home and Away should go into the model. My first guess was just to fit a simple linear regression model, like so:

fit <- lm(Score ~ Home + Away + Season, data=dat)
summary(fit)

Call:
lm(formula = Score ~ Home + Away + Season, data = dat)

Residuals:
Min       1Q   Median       3Q      Max
-10.2750  -2.3604  -0.2375   2.9792  10.8583

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.9333     2.3096   0.404    0.688
HomeB         1.0000     1.8518   0.540    0.592
HomeC        -0.3750     1.8518  -0.203    0.840
HomeD         0.0750     1.8518   0.041    0.968
AwayB         1.0000     1.8518   0.540    0.592
AwayC        -1.5250     1.8518  -0.824    0.414
AwayD         0.4250     1.8518   0.230    0.819
Season2015   -2.0833     1.9520  -1.067    0.291
Season2016   -0.2500     1.9520  -0.128    0.899
Season2017   -1.4167     1.9520  -0.726    0.471
Season2018    0.2500     1.9520   0.128    0.899

Residual standard error: 4.781 on 49 degrees of freedom
Multiple R-squared:  0.08193,   Adjusted R-squared:  -0.1054
F-statistic: 0.4373 on 10 and 49 DF,  p-value: 0.9208


(never mind that the model is not performing at all – obviously in this league, every outcome is completely determined by accident.)

However, I see at at least two conceptual problems with the model. First, the intercept represents the average score for home team A playing away team A in season 2014, which is an event that is not in the observed data, and may never be observed. Second, the model doesn't seem to be aware really of the round-robin format. These thoughts make me suspicious that there is something fundamentally wrong with this approach, but I can't think of a more suitable model structure.

Is this an acceptable way to model a data set like this? Is there a better model specification that I should use instead that reflects the round-robin? Are there more appropriate models for this type of data?