# How to specify a hierarchical bayesian model with sum-to-zero constraints?

I'm working on the first model described in this paper ("Bayesian hierarchical model for the prediction of football [soccer] results"). The gist of the model is:

The model includes two sum-to-zero constraints, on the team-level attack and defense parameters. Below is the winBUGS code verbatim from the paper. Note that they enforce the sum-to-zero constraint by subtracting the mean.

model {
# LIKELIHOOD AND RANDOM EFFECT MODEL FOR THE SCORING PROPENSITY
for (g in 1:ngames) {
# Observed number of goals scored by each team
y1[g] ~ dpois(theta[g,1])
y2[g] ~ dpois(theta[g,2])
# Predictive distribution for the number of goals scored
ynew[g,1] ~ dpois(theta[g,1])
ynew[g,2] ~ dpois(theta[g,2])
# Average Scoring intensities (accounting for mixing components)
log(theta[g,1]) <- home + att[hometeam[g]] + def[awayteam[g]]
log(theta[g,2]) <- att[awayteam[g]] + def[hometeam[g]]
}

# 1. BASIC MODEL FOR THE HYPERPARAMETERS # prior on the home effect
home ~ dnorm(0,0.0001)

# Trick to code the ‘‘sum-to-zero’’ constraint for (t in 1:nteams) {
att.star[t] ~ dnorm(mu.att,tau.att)
def.star[t] ~ dnorm(mu.def,tau.def)
att[t] <- att.star[t] - mean(att.star[])
def[t] <- def.star[t] - mean(def.star[]) }

# priors on the random effects
mu.att ~ dnorm(0,0.0001)
mu.def ~ dnorm(0,0.0001)
tau.att ~ dgamma(.01,.01)
tau.def ~ dgamma(.01,.01)

}


My questions are:

1. Does it make sense to have mu.att be a hyper parameter, when the we're enforcing the sum-to-zero constraint by subtracting the mean anyway? When I run this model in python (code below), I'm able to reproduce their results, but sure enough, mu_att and mu_def don't converge.

2. Does the model need an intercept, representing the average goal scoring intensity? I've seen some other examples of models with sum-to-zero constraints, and they included an intercept.

3. If I were to switch sum-to-zero enforcement techniques, from subtract-the-mean to set-first-equal-to-negative-sum-of-remaining, what would that mean for questions 1. and 2.?

For python folks, here is my attempt at reproducing this in pymc. First, I extract arrays from the data:

observed_home_goals = [row['home_score'] for i, row in df.iterrows()]
observed_away_goals = [row['away_score'] for i, row in df.iterrows()]
who_played_whom = [(row['i_home'], row['i_away']) for i,row in df.iterrows()]
num_teams = len(df.i_home.unique())
num_games = len(who_played_whom)


and then here is my model:

home = pymc.Normal('home', 0, .0001, value=0)
mu_att = pymc.Normal('mu_att', 0, .0001, value=0)
mu_def = pymc.Normal('mu_def', 0, .0001, value=0)
tau_att = pymc.Gamma('tau_att', .1, .1)
tau_def = pymc.Gamma('tau_def', .1, .1)

#team-specific parameters
atts_star = pymc.Normal("atts_star", mu=mu_att, tau=tau_att, size=num_teams)
defs_star = pymc.Normal("defs_star", mu=mu_def, tau=tau_def, size=num_teams)

# trick to code the sum to zero contraint
@pymc.deterministic
def atts(atts_star=atts_star):
atts = atts_star.copy()
atts = atts - np.mean(atts_star)
return atts

@pymc.deterministic
def defs(defs_star=defs_star):
defs = defs_star.copy()
defs = defs - np.mean(defs_star)
return defs

@pymc.deterministic
def home_theta(who_played_whom=who_played_whom, home=home,
atts=atts, defs=defs):
home_attack = [atts[i[0]] for i in who_played_whom]
away_defense = [defs[i[1]] for i in who_played_whom]
return [math.exp(home + home_attack[i] + away_defense[i]) for i in range(num_games)]

@pymc.deterministic
def away_theta(who_played_whom=who_played_whom,
atts=atts, defs=defs):
away_attack = [atts[i[1]] for i in who_played_whom]
home_defense = [defs[i[0]] for i in who_played_whom]
return [math.exp(away_attack[i] + home_defense[i]) for i in range(num_games)]

home_goals = pymc.Poisson('home_goals', mu=home_theta, value=observed_home_goals, observed=True)
away_goals = pymc.Poisson('away_goals', mu=away_theta, value=observed_away_goals, observed=True)

mcmc = pymc.MCMC([home, mu_att, mu_def, tau_att, tau_def,
home_theta, away_theta, atts_star, defs_star, atts, defs, home_goals, away_goals])
map_ = pymc.MAP( mcmc )
map_.fit()
mcmc.sample(100000, 10000, 4)


2. Yes. It appears that they over-constrained the model. The symmetry in this model is that you can add $c$ to all $att_t$ and subtract $c$ from all $def_t$. This is one degree of freedom. But they imposed two constraints, removing two degrees of freedom. So you need to put one degree of freedom back. Alternatively, get rid of their two constraints and use the single constraint $\sum_t att_t = \sum_t def_t$.