How to implement Bayesian hierarchical model for the prediction of football results in R?

Question

I'm trying to implement this idea in R and I use the English Premier League as my data set https://discovery.ucl.ac.uk/id/eprint/16040/1/16040.pdf

In page 4, it is stated that

So I tried to set the attack and defence ability of the first team to zero but that is not meaningful and I want to add the sum-to-zero constraint. How can I do that in R?

Below is the entire code

# Load the data
my_df = read.csv('https://www.football-data.co.uk/mmz4281/2223/E0.csv')
my_df = my_df[, c('HomeTeam', 'AwayTeam', 'FTHG', 'FTAG')]

# Data preparation
my_df$HomeTeam=as.character(my_df$HomeTeam)
my_df$AwayTeam=as.character(my_df$AwayTeam)

teams = unique(c(my_df$HomeTeam, my_df$AwayTeam))

# JAGS list with the data from my_df in integer strings
DATAlist = list(FTHG = my_df$FTHG, FTAG = my_df$FTAG, 
                HomeTeam = as.numeric(factor(my_df$HomeTeam, levels=teams)),
                AwayTeam = as.numeric(factor(my_df$AwayTeam, levels=teams)),
                n_teams = length(teams), n_games = nrow(my_df))

# Generator for the type of column names Jags outputs.
column_name = function(name, ...) {
  paste0(name, "[", paste(..., sep=",") , "]")
}

# Model string
model_string = "model {
for(i in 1:n_games) {
  FTHG[i] ~ dpois(lambda_home[HomeTeam[i],AwayTeam[i]])
  FTAG[i] ~ dpois(lambda_away[HomeTeam[i],AwayTeam[i]])
}

for(HomeTeam_i in 1:n_teams) {
  for(AwayTeam_i in 1:n_teams) {
    lambda_home[HomeTeam_i, AwayTeam_i] = exp(int + hfa + 
                                    attack[HomeTeam_i] - defence[AwayTeam_i])
    lambda_away[HomeTeam_i, AwayTeam_i] = exp(int + 
                                    attack[AwayTeam_i] - defence[HomeTeam_i])
  }
}

attack[1] = 0
defence[1] = 0

for(j in 2:n_teams) {
  attack[j] ~ dnorm(group_skill, group_tau)
  defence[j] ~ dnorm(group_skill, group_tau)
}  

group_skill ~ dnorm(0, 0.0625)
group_tau = 1 / pow(group_sigma, 2)
group_sigma ~ dunif(0, 3)
int ~ dnorm(0, 0.0625)
hfa ~ dnorm(0, 0.0625)
}"

bayes_model = jags.model(textConnection(model_string), 
                    data=DATAlist, n.chains=3, n.adapt=1000)
update(bayes_model, 1000)
params = c("int",'hfa',  "attack", 'defence', "group_skill", "group_sigma")
bayes_mod_sim = coda.samples(bayes_model, variable.names = params, n.iter = 1e4)
bayes_mod_sim = as.matrix(bayes_mod_sim)

Check out bayes_mod_sim, you will see that the attack and defence ability of the first team is zero.