[This article was first published on Wiekvoet, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
In this post the football model is programmed into JAGS. There are all the reasons to do so. Jags 3.3 is recently released, I was stimulated by Gianluca‘s post . Obviously I could copy the model in his paper, but that would be too easy and not sure about copyright either. So, in this post, variations of the model of an earlier post is recreated.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Data structure
The starting data looks like this:
head(Jold)
Seizoen Datum Thuisclub Uitclub Thuisscore Uitscore
1 2011-2012 2011-08-05 Excelsior Feyenoord 0 2
2 2011-2012 2011-08-06 RKC Waalwijk Heracles Almelo 2 2
3 2011-2012 2011-08-06 Roda JC FC Groningen 2 1
4 2011-2012 2011-08-06 SC Heerenveen NEC 2 2
5 2011-2012 2011-08-06 VVV-Venlo FC Utrecht 0 0
6 2011-2012 2011-08-07 ADO Den Haag Vitesse 0 0
The data has to be made ready to be placed into JAGS. All factors have to be made numbers. Note that each game is still one row. The two outcomes will be resolved in JAGS
JAGSData <- with(Jold,list(
N=length(Datum),
nclub = nlevels(Thuisclub),
HomeClub=(1:nlevels(Thuisclub))[Thuisclub],
OutClub=(1:nlevels(Uitclub))[Uitclub],
HomeMadeGoals=Thuisscore,
OutMadeGoals=Uitscore))
str(JAGSData)
List of 6
$ N : int 306
$ nclub : int 18
$ HomeClub : int [1:306] 5 14 15 16 18 1 3 4 11 7 …
$ OutClub : int [1:306] 9 10 6 12 8 17 13 2 7 3 …
$ HomeMadeGoals: int [1:306] 0 2 2 2 0 0 3 1 0 2 …
$ OutMadeGoals : int [1:306] 2 2 1 2 0 0 1 4 1 0 …
Model
The model in JAGS has to contain all the details and all distributions. That makes it much longer, but also gives the opportunity to fiddle around with details. For now it is reasonably simple. Each club has two properties, an attack strength (Astr) and a defense strength (Dstr). Based on the combination of these strengths is the number of goals, again Poisson distributed. The interesting part is in the way attack and defense strength are constructed. These are created from an overall strength (TopStr) and an attack/defense balance (AD).fbmodel1 <- function() {
for (i in 1:N) {
HomeMadeGoals[i] ~ dpois(HS[i])
OutMadeGoals[i] ~ dpois(OS[i])
log(HS[i]) <- Home + Astr[HomeClub[i]] + Dstr[OutClub[i]]
log(OS[i]) <- Dstr[HomeClub[i]] + Astr[OutClub[i]]
}
for (i in 1:nclub) {
Astr[i] <- TopStr[i]+AD[i]
Dstr[i] <- TopStr[i]-AD[i]
TopStr[i] ~ dnorm(0,tauTop)
AD[i] ~ dnorm(0,tauAD)
}
tauAD <- pow(sigmaAD,-2)
sigmaAD ~ dunif(0,100)
tauTop <- pow(sigmaTop,-2)
sigmaTop ~ dunif(0,100)
Home ~ dnorm(0,.0001)
}
To make this model work it needs a series on incantations in R. To keep things short only a basic output plot is given
params <- c(“TopStr”,”AD”,”sigmaAD”,”sigmaTop”,”Home”)
inits <- function(){
list(TopStr=rnorm(JAGSData$nclub),
AD=rnorm(JAGSData$nclub),
sigmaAD=runif(1),
sigmaTop=runif(0,.5),
Home=rnorm(1)
)
}
jagsfit <- jags(JAGSData, model=fbmodel1, inits=inits,
parameters=params,progress.bar=”gui”,
n.iter=4000)
plot(jagsfit)
inits <- function(){
list(TopStr=rnorm(JAGSData$nclub),
AD=rnorm(JAGSData$nclub),
sigmaAD=runif(1),
sigmaTop=runif(0,.5),
Home=rnorm(1)
)
}
jagsfit <- jags(JAGSData, model=fbmodel1, inits=inits,
parameters=params,progress.bar=”gui”,
n.iter=4000)
plot(jagsfit)
Parameter extraction
As the interest is in the teams, their properties are extracted. As Gianluca wrote we considered a slightly more complex structure in which we included information on each team’s propensity to be “good”, “average”, or “poor”. This helped avoid overshrinkage in the estimations To examine this I looked at the teams relative scores.The plots look like they have a shoulder, so there is a point there.
jags.smat <- jagsfit$BUGSoutput$summary
for (i in 1:nlevels(Jold$Thuisclub))
rownames(jags.smat) <- sub(paste(‘[‘,i,’]’,sep=”),levels(Jold$Thuisclub)[i],
rownames(jags.smat),fixed=TRUE)
plot(density(jags.smat[grep(‘^Top’,rownames(jags.smat)),1])
,main=’Distribution of Strength’)
points(x=jags.smat[grep(‘^Top’,rownames(jags.smat)),1],
y=rep(0,nlevels(Jold$Thuisclub)))
plot(density(jags.smat[grep(‘^AD’,rownames(jags.smat)),1]),
main=’Distribution of Attack / Defense’)
points(x=jags.smat[ grep(‘^AD’,rownames(jags.smat)),1],
y=rep(0,nlevels(Jold$Thuisclub)))
model variation
There are a number of ways to create distributions for attack and defense strength (Astr and Dstr). One variation is to make them independent, but also to give them fatter tails, which is handled by a t distribution rather than a normal distribution. This fatter tail could be accommodate teams much better or worse than the majority. Unfortunately this did not work, the number of degrees of freedom for the t distribution varied a lot and was too high to make this an alternative model.
# section replaced in fbmodel1
for (i in 1:nclub) {
Astr[i] ~ dt(0,tauStr,nuStr)
Dstr[i] ~ dt(0,tauStr,nuStr)
}
tauStr <- pow(sigmaStr,-2)
sigmaStr ~ dunif(0,100)
nuStr <- 1/InuStr
InuStr ~ dunif(0,.5)
To leave a comment for the author, please follow the link and comment on their blog: Wiekvoet.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.