A dynamic rating for Brazilian football league 2011 season

##----------------------------------------------------------------------
## Required packages and session definition.

require(latticeExtra)
require(plyr)
require(reshape2)
require(gridExtra)

sessionInfo()
Sys.time()

##----------------------------------------------------------------------

rm(list=ls())

## Symbols.
## plot(1:30, pch=1:30)

## Color palette.
mycol <- brewer.pal(6, "Set1")

## Trellis graphical style.
ps <- list(
    box.rectangle=list(col=1, fill=c("gray70")),
    box.umbrella=list(col=1, lty=1),
    dot.symbol=list(col=1, pch=19),
    dot.line=list(col="gray50", lty=3),
    plot.symbol=list(col=1, cex=1.2),
    plot.line=list(col=1),
    plot.polygon=list(col="gray95"),
    superpose.line=list(col=mycol),
    superpose.symbol=list(col=mycol, pch=c(1:5)),
    superpose.polygon=list(col=mycol),
    strip.background=list(col=c("gray80","gray50"))
    )
trellis.par.set(ps)
## show.settings()

##----------------------------------------------------------------------
## Load the data of the season 2011.

## UNIX command to get the charset/encoding of a file.
system("file -bi brasileiro2004-2011.txt")

## iso-8859-1 aka latin1.

## Importing the dataset with seasons 2004 to 2011.
da <- read.table(
    file="brasileiro2004-2011.txt",
    header=TRUE,
    sep="\t",
    encoding="latin1",
    colClasses=c("factor","character","integer","NULL")[
        c(2,4,3,4,4,4,1,1,2)])

## Short names are more manageable.
## man: home team, away: away team.
names(da) <- substr(names(da), 1, 3)
names(da)[3:4] <- c("home", "away")
str(da)

## Same order in the factor levels.
l <- union(levels(da$home), levels(da$away))
da$home <- factor(da$home, levels=l)
da$away <- factor(da$away, levels=l)

## Converting string to data (calendar format).
da$dat <- as.POSIXct(x=da$dat, format="%d/%m/%Y")
da$yea <- as.integer(format(da$dat, format="%Y"))
    
## Get the goals of each team in a match.
da$h <- as.integer(gsub("^(\\d+)x\\d+$", "\\1", tolower(da$pla)))
da$a <- as.integer(gsub("^\\d+x(\\d+)$", "\\1", tolower(da$pla)))
da <- transform(da, d=abs(h-a), s=h+a)
str(da)

## Marginal statistics in goals scored for each season.
aggregate(cbind(home=h, away=a, absdiff=d, total=s)~yea,
          data=da, FUN=mean, na.rm=TRUE)

##----------------------------------------------------------------------

## 2010 season only.
db <- droplevels(subset(da, yea==2010,
                        select=c("rod","home","away","h","a")))
str(db)

subset(db, home=="Corinthians")

## Are the levels in the same order?
all(levels(db$home)==levels(db$away))

##----------------------------------------------------------------------

## Matrices to identify teams matches (design matrices).
Xh <- model.matrix(~-1+home, db) ## h: home.
Xa <- model.matrix(~-1+away, db) ## a: away.

## Check.
Xha <- Xh-Xa
db[1:5,c("home","away")]

print(as.table(t(Xha[1:5,])), zero.print=".")

##----------------------------------------------------------------------

## Place in 2010 season, from first to last.
cbind(scored=t(Xh)%*%db$h+t(Xa)%*%db$a,   ## goals scored
      conceded=t(Xa)%*%db$h+t(Xh)%*%db$a) ## goals conceded

## http://esportes.terra.com.br/futebol/brasileiro/2010/noticias/0,,OI4339585-EI15413,00-Classificacao+e+Jogos+Serie+A.html

aux <- with(db, {
    d <- h-a
    draw <- d==0
    winH <- d>0
    ## winA <- !draw & !winH
    winA <- !(draw | winH)
    x <- cbind(h=winH*3+draw, a=winA*3+draw)
    return(x)
})

## cbind(db, aux)[sample(1:nrow(db), size=10),]

pl10 <- t(Xh)%*%aux[,"h"]+t(Xa)%*%aux[,"a"]
pl10 <- cbind(team=levels(db$home), data.frame(pts=pl10))
## pl10$rk <- rank(pl10$pts, decreasing=TRUE)
pl10 <- arrange(pl10, -pts)
pl10$pos <- rank(-pl10$pts)
pl10

th16 <- pl10[16,]

rm(aux)

##----------------------------------------------------------------------
## Likelihood definition.

## log-likeligood funtion corresponding to a Poisson model for the goals
## scored as function of teams attacking streng (gamma), teams
## deffensive strength (delta), and home(away) field
## advantage(disadvatage) (omega,kappa), so theta=c(gamma, delta,
## omega, kappa).

## attack, defense and home field advantage.
ll1 <- function(theta, gh, ga, Xh, Xa){
    ## theta: parameter vector.
    ## gh, ga: goals for the home and away team.
    ## Xh, Xa: design matrices for home and away team.
    gamma <- 1:20; delta <- 21:40; omega <- 41
    eta.h <- (Xh%*%theta[gamma]-Xa%*%theta[delta])+theta[omega]
    eta.a <- (Xa%*%theta[gamma]-Xh%*%theta[delta])
    ll.h <- dpois(gh, lambda=exp(eta.h), log=TRUE)
    ll.a <- dpois(ga, lambda=exp(eta.a), log=TRUE)
    ll <- sum(ll.h+ll.a)
    return(ll)
}

## attack, defense, home field advantage and away field disadvantage.
ll2 <- function(theta, gh, ga, Xh, Xa){
    ## theta: parameter vector.
    ## gh, ga: goals for the home and away team.
    ## Xh, Xa: design matrices for home and away team.
    gamma <- 1:20; delta <- 21:40; omega <- 41; kappa <- 42
    eta.h <- (Xh%*%theta[gamma]-Xa%*%theta[delta])+theta[omega]
    eta.a <- (Xa%*%theta[gamma]-Xh%*%theta[delta])-theta[kappa]
    ll.h <- dpois(gh, lambda=exp(eta.h), log=TRUE)
    ll.a <- dpois(ga, lambda=exp(eta.a), log=TRUE)
    ll <- sum(ll.h+ll.a)
    return(ll)
}

##----------------------------------------------------------------------
## Optimization to get parameter estimates.

op1 <- optim(par=rep(1, 41),
             fn=ll1,
             gh=db$h, ga=db$a,
             Xh=Xh, Xa=Xa,
             control=list(fnscale=-1),
             hessian=TRUE)
str(op1)

op2 <- optim(par=rep(1, 42),
             fn=ll2,
             gh=db$h, ga=db$a,
             Xh=Xh, Xa=Xa,
             control=list(fnscale=-1),
             hessian=TRUE)
str(op2)

## Likelihoods.
l <- c(ll1=op1$value, ll2=op2$value)
l

pchisq(-2*diff(l), df=1, lower.tail=FALSE)

theta <- data.frame(th1=c(op1$par, NA), th2=op2$par)
theta

## plot(th2~th1, data=theta[1:40,])
## a <- with(theta, identify(th1, th2))
## dput(a)
a <- c(2L, 6L, 12L, 14L, 16L, 22L)
theta <- cbind(theta[1:40,],
               expand.grid(team=levels(db$home),
                           param=c("gamma","delta")))
theta[a,]

## Legend for gamma and delta.
key.gd <- list(columns=2,
               title="Parameter", cex.title=1.1,
               points=list(
                   col=ps$superpose.line$col[1:2],
                   pch=ps$superpose.symbol$pch[1:2]),
               text=list(c(expression(gamma), expression(delta))))

xyplot(th2~th1, data=theta, groups=param,
       key=key.gd, aspect="iso",
       xlab="Parameter estimates in the 1 field parameter model",
       ylab="Parameter estimates in the 2 field parameter model")+
    layer(panel.abline(a=0, b=1))+
             layer(panel.text(x=theta$th1[a], y=theta$th2[a], pos=4,
                              labels=as.character(theta$team[a])))

##----------------------------------------------------------------------

## Precision matrix.
## levelplot(-op$hessian[-41,-41])

## Variace-covariance matrix.
## levelplot(-solve(op$hessian)[-41,-41])

## Standard error of parameter estimates (so big).
## sqrt(diag(-solve(op$hessian)))

priors1 <- expand.grid(teams=levels(db$home),
                       par=c("gamma", "delta"),
                       stringsAsFactors=FALSE)
priors1 <- rbind(priors1, c("homefield", "omega"))

priors1$est <- sqrt(op1$par)        ## Estimate.
priors1$info <- diag(-op1$hessian)  ## Precision.
priors1$var <- 1/priors1$info       ## Inverse precision.
priors1$sd <- sqrt(priors1$var)

## Spliting the table of priors1.
L <- split(subset(priors1, select=c("teams","est","var")),
           f=priors1$par)

## Using better names.
names(L$gamma)[-1] <- c("mu","sigma2")
names(L$delta)[-1] <- c("theta","psi2")
names(L$omega)[-1] <- c("mu0","sigma20")
priors1 <- merge(L$gamma, L$delta)
priors1$teams <- factor(priors1$teams, levels=levels(db$home))
attr(priors1, "hf") <- list(omega=L$omega)

arrange(priors1, -mu)

attr(priors1, "hf")

##----------------------------------------------------------------------

priors2 <- expand.grid(teams=levels(db$home),
                       par=c("gamma", "delta"),
                       stringsAsFactors=FALSE)
priors2 <- rbind(priors2,
                 c("homefield", "omega"),
                 c("awayfield", "kappa"))

priors2$est <- sqrt(op2$par)        ## Estimate.
priors2$info <- diag(-op2$hessian)  ## Precision.
priors2$var <- 1/priors2$info       ## Inverse precision.
priors2$sd <- sqrt(priors2$var)

## Spliting the table of priors2.
M <- split(subset(priors2, select=c("teams","est","var")),
           f=priors2$par)

## Using better names.
names(M$gamma)[-1] <- c("mu","sigma2")
names(M$delta)[-1] <- c("theta","psi2")
names(M$omega)[-1] <- c("mu0","sigma20")
names(M$kappa)[-1] <- c("theta0","psi20")
priors2 <- merge(M$gamma, M$delta)
priors2$teams <- factor(priors2$teams, levels=levels(db$home))
attr(priors2, "hf") <- list(omega=M$omega, kappa=M$kappa)

arrange(priors2, -mu)

attr(priors2, "hf")

##----------------------------------------------------------------------

## Just the 2011 season.
dc <- droplevels(subset(da, yea==2011))
rownames(dc) <- NULL
str(dc)

##----------------------------------------------------------------------

## http://esportes.terra.com.br/futebol/brasileiro/2011/seriea/classificacao_jogos/

## Rating and position in 2011 season.
Xh <- model.matrix(~-1+home, dc) ## h: home.
Xa <- model.matrix(~-1+away, dc) ## a: away.

cbind(scored=t(Xh)%*%dc$h+t(Xa)%*%dc$a,   ## goals scored
      conceded=t(Xa)%*%dc$h+t(Xh)%*%dc$a) ## goals conceded

aux <- with(dc, {
    d <- h-a
    draw <- d==0
    winH <- d>0
    ## winA <- !draw & !winH
    winA <- !(draw | winH)
    x <- cbind(h=winH*3+draw, a=winA*3+draw)
    return(x)
})

pl11 <- t(Xh)%*%aux[,"h"]+t(Xa)%*%aux[,"a"]
pl11 <- cbind(team=levels(dc$home), data.frame(pts=pl11))
pl11 <- arrange(pl11, -pts)
pl11$pos <- rank(-pl11$pts)
## pl11

arrange(merge(pl10, pl11, by="team", all=TRUE), -pts.y)

##----------------------------------------------------------------------

## Distribution of goals.
histogram(~h+a+d+s, data=dc,
          breaks=with(dc, seq(min(d), max(s)+1, by=1)-0.5),
          xlab="Goals",
          ylab="Relative frequency (%)",
          strip=strip.custom(
              factor.levels=c("Goals of Home","Goals of Away",
                  "Absolute difference","Total of goals")))

## Classes of goals: 0, 1, ..., 3, >3.
cutin <- function(z){
    cut(z, breaks=c(0:4, Inf),
        right=FALSE, include.lowest=TRUE)
}

## Crossed frequency.
xt <- xtabs(~cutin(h)+cutin(a), data=dc)
addmargins(xt)

## Remove objects.
rm(xt, cutin)

##----------------------------------------------------------------------
## Set operation.

## Teams in 2010 that aren't in 2011 (the worst in seria A).
wA <- setdiff(levels(db$home), levels(dc$home))
wA

## The bests of serie B, came to 2011 season.
bB <- setdiff(levels(dc$home), levels(db$home))
bB

##----------------------------------------------------------------------
## Prior for the new teams, those in serie B that went to seria A.

priors1 <- droplevels(subset(priors1, !teams%in%wA))

## The team at 16th position.
t16 <- subset(priors1, teams==as.character(th16$team))

t16 <- t16[rep(1, 4),]
t16$teams <- bB

## Priors1 for all teams in the 2011 season.
priors1 <- rbind(priors1, t16)
rownames(priors1) <- NULL
priors1

## Prior mean and variance for the home-field advantage parameter.
attr(priors1, "hf") <- list(omega=L$omega)
attr(priors1, "hf")

## Same order in the factor levels.
priors1$teams <- factor(as.character(priors1$teams),
                        levels=levels(dc$home))

##----------------------------------------------------------------------

priors2 <- droplevels(subset(priors2, !teams%in%wA))

## Priors2 for all teams in the 2011 season.
priors2 <- rbind(priors2, t16)
rownames(priors2) <- NULL
priors2

## Prior mean and variance for the home-field advantage parameter.
attr(priors2, "hf") <- list(omega=M$omega, kappa=M$kappa)
attr(priors2, "hf")

## Same order in the factor levels.
priors2$teams <- factor(as.character(priors2$teams),
                        levels=levels(dc$home))

## Remove objects.
rm(wA, bB, t16, L, M)

## Check.
## all.equal(levels(db$home), levels(db$away))
## cbind(levels(dc$home), levels(dc$away), levels(priors1$teams))

##----------------------------------------------------------------------

priors12 <- rbind(cbind(priors1, pr=1),
                  cbind(priors2, pr=2))

key.gd <- list(columns=1,
               title="Model", cex.title=1.1,
               points=list(
                   col=ps$superpose.line$col[1:2],
                   pch=c(1,5)),
               text=list(c("1 field parameter",
                   "2 field parameter")))

xyplot(mu+sigma2+theta+psi2~teams, outer=TRUE,
       data=priors12, groups=pr, pch=c(1,5), key=key.gd,
       xlab="Teams", ylab="Prior", as.table=TRUE,
       scales=list(y="free", x=list(rot=90)),
       strip=strip.custom(
           factor.levels=expression(
               mu*": prior mean for"~gamma,
               sigma^2*": prior variance for"~gamma,
               theta*": prior mean for"~delta,
               psi^2*": prior variance for"~delta)))

##----------------------------------------------------------------------
## Prior densities.

## str(priors12)
## xtabs(~teams+pr, priors12)

L <- apply(priors12[,-1], MARGIN=1,
           FUN=function(pars){
               x <- seq(0.2, 1.8, length.out=50)
               ygamma <- dnorm(x, mean=pars[1], sd=sqrt(pars[2]))
               ydelta <- dnorm(x, mean=pars[3], sd=sqrt(pars[4]))
               data.frame(x=x, ygamma=ygamma, ydelta=ydelta, pr=pars[5])
           })

names(L) <- rep(levels(priors12$teams), times=2)
L <- ldply(L, .id="teams")
str(L)

## Legend for gamma and delta.
key.gd <- list(columns=1,
               title="Parameter", cex.title=1.1,
               lines=list(
                   col=ps$superpose.line$col[1:2]),
               text=list(c(expression(gamma), expression(delta))))

## Legend for 1 and 2 field parameter models.
key.12 <- list(columns=1,
               title="Model", cex.title=1.1,
               lines=list(
                   lty=1:2),
               text=list(c("1 field parameter",
                   "2 field parameter")))

## Create a grob by merging the two keys (latticeExtra).
keys <- mergedTrellisLegendGrob(
    a=list(fun=draw.key, args=list(key=key.gd)),
    b=list(fun=draw.key, args=list(key=key.12)),
    vertical=FALSE)

xyplot(ygamma+ydelta~x|teams, data=subset(L, pr==1),
       type="l", as.table=TRUE,
       xlab="Parameter space",
       ylab="Prior density",
       legend=list(top=list(fun=keys))
       )+
    as.layer(xyplot(ygamma+ydelta~x|teams, data=subset(L, pr==2),
                    type="l", lty=2))

## About multiple legends:
## https://procomun.wordpress.com/2012/02/18/maps_with_r_1/

##----------------------------------------------------------------------

## Number of teams and number of rounds.
nteams <- nlevels(db$home)
nrounds <- length(unique(db$rod))
c(nteams, nrounds)

## gamma: teams attacking strength.
## delta: teams deffensive strength.
## omega: home-field advantage parameter.
gamma <- matrix(numeric(), nrow=nteams, ncol=nrounds+1)
delta <- gamma
omega <- numeric(nrounds)

## To store values in each match inside a round.
lambdaH <- matrix(numeric(), nrow=nteams/2, ncol=nrounds)
lambdaA <- lambdaH

i <- 1
gamma[,i] <- priors1$mu
delta[,i] <- priors1$theta
omega[i] <- attr(priors1, "hf")$omega$mu0

## Fitting for each round.
for(i in 1:nrounds){
    ##-------------------------------------------
    ## Subset each round (10 matches).
    ## 
    dbi <- subset(dc, rod==i)
    ##
    ##-------------------------------------------
    ## Model matrices, which teams confronted each other.
    ## 
    Xh <- model.matrix(~-1+home, data=dbi)
    Xa <- model.matrix(~-1+away, data=dbi)
    ##
    ##-------------------------------------------
    ## Columns of Home and Away teams.
    ##
    posh <- apply(apply(Xh, MARGIN=1, FUN="==", 1), MARGIN=2, FUN=which)
    posa <- apply(apply(Xa, MARGIN=1, FUN="==", 1), MARGIN=2, FUN=which)
    ##
    ##-------------------------------------------
    ## Lambda values.
    ##
    lambdaH[,i] <- as.vector(
        exp(Xh%*%gamma[,i]-Xa%*%delta[,i]+omega[i])
        )
    lambdaA[,i] <- as.vector(
        exp(Xa%*%gamma[,i]-Xh%*%delta[,i])
        )
    ##
    ##-------------------------------------------
    ## Goals expected and goals scored.
    ## 
    ## cbind(goalRateHome=lambdaA[,i], goalRateAway=lambdaB[,i],
    ##       goalScorHome=dbi$x, goalScorAway=dbi$y)
    ##
    ##-------------------------------------------
    ## Updating priors for the next round.
    ##
    gamma[posh,i+1] <- 
        Xh%*%gamma[,i]+Xh%*%priors1$sigma2*(dbi$h-lambdaH[,i])
    gamma[posa,i+1] <- 
        Xa%*%gamma[,i]+Xa%*%priors1$sigma2*(dbi$a-lambdaA[,i])
    delta[posh,i+1] <- 
        Xh%*%delta[,i]-Xh%*%priors1$psi2*(dbi$h-lambdaA[,i])
    delta[posa,i+1] <- 
        Xa%*%delta[,i]-Xa%*%priors1$psi2*(dbi$h-lambdaH[,i])
    omega[i+1] <- 
        omega[i]+attr(priors1, "hf")$omega$sigma20*
            sum(dbi$h-lambdaH[,i])
}

gamma <- gamma[,-1]
delta <- delta[,-1]
omega <- omega[-1]

##----------------------------------------------------------------------

## Series for all teams parameters along the rounds.
dd <- melt(data=as.data.frame(delta), value.name="delta")
dd <- cbind(
    melt(data=cbind(team=levels(dc$home), as.data.frame(gamma)),
         id.vars="team", variable.name="rod", value.name ="gamma"),
    delta=dd$delta)
dd$rod <- as.integer(dd$rod)
str(dd)

sc <- list(x=list(at=c(1,10,20,30,38)))

## Legend for gamma and delta.
key.gd <- list(columns=2,
               title="Parameter", cex.title=1.1,
               type="o", divide=1,
               lines=list(
                   col=ps$superpose.line$col[1:2],
                   pch=ps$superpose.symbol$pch[1:2]),
               text=list(c(expression(gamma), expression(delta))))

## All teams.
xyplot(gamma+delta~rod|team, data=dd, type="l",
       as.table=TRUE,
       scales=sc,
       key=key.gd[c(-4,-5)],
       xlab="Championship round",
       ylab="Parameter value")

## Home field parameter.
xyplot(omega~seq_along(omega), type="b",
       scales=sc,
       xlab="Championship round",
       ylab="Home field advantage parameter")

##----------------------------------------------------------------------

## Select som teams to a depper view, 2 bests and 2 worsts.
sel <- levels(dd$team)[c(1,5,9,20)]
de <- droplevels(subset(dd, team%in%sel))

## Selected teams.
xyplot(gamma+delta~rod|team,
       data=de, type="o",
       as.table=TRUE,
       scales=sc,
       key=key.gd,
       xlab="Championship round",
       ylab="Parameter value")

key.tm <- list(columns=2,
               title="Teams", cex.title=1.1,
               type="o", divide=1,
               lines=list(
                   col=ps$superpose.line$col[1:4],
                   pch=ps$superpose.symbol$pch[1:4]),
               text=list(levels(de$team)))

xyplot(gamma+delta~rod, outer=TRUE, groups=team,
       data=de, type="o",
       as.table=TRUE,
       scales=sc,
       key=key.tm,
       strip=strip.custom(factor.levels=key.gd$text[[1]]),
       xlab="Championship round",
       ylab="Parameter value")

##----------------------------------------------------------------------
## Rating series.

## Rating along the matches.
dd$rating <- with(dd, gamma+delta)
de$rating <- with(de, gamma+delta)

xyplot(rating~rod|team,
       data=dd, type="o", cex=0.8,
       as.table=TRUE,
       scales=sc,
       xlab="Championship round",
       ylab=expression("Rating"~(gamma+delta)))

xyplot(rating~rod, groups=team,
       data=de, type="o", cex=0.8,
       as.table=TRUE,
       scales=sc,
       key=key.tm,
       strip=strip.custom(factor.levels=key.gd$text[[1]]),
       xlab="Championship round",
       ylab=expression("Rating"~(gamma+delta)))

##----------------------------------------------------------------------

## Average rating.
rtm <- arrange(aggregate(rating~team, data=dd, FUN=mean), -rating)
names(rtm)[2] <- "averating"

## Final rating.
frat <- subset(dd, rod==max(dd$rod))
names(frat)[5] <- "finrating"

## Final, average rating and rank.
fr <- arrange(merge(merge(pl11, frat, by="team"),
                    rtm, by="team"), -finrating)
fr

xyplot(finrating+averating~pts, data=fr,
       auto.key=TRUE, type=c("p","smooth"))

splom(fr[,c("pts","finrating","averating")],
      type=c("p","smooth"))

cor(fr[,c("pts","finrating","averating")],
    method="spearman")

## Como calcular as probabilidades de ganhar, empatar e perder?

##-----------------------------------------------------------------------------
## Who have more chance to gain, Palmeiras or Corinthians?

prob <- function(gh, ga, theta, hf=FALSE){
    ## 1: attack home
    ## 2: defense home
    ## 3: attack away
    ## 4: defense away
    ## 5: home field advantage
    if(hf){
        pa <- dpois(gh, lambda=exp(theta[5]+theta[1]-theta[4]))
    } else {
        pa <- dpois(gh, lambda=exp(theta[1]-theta[4]))
    }
    pb <- dpois(ga, lambda=exp(theta[3]-theta[2]))
    pa*pb
}
Prob <- Vectorize(prob, c("gh", "ga"))
rm("prob")

i <- c("Avaí", "Corinthians")
i <- rev(i)

sel <- subset(fr, team%in%i,
              select=c("team", "gamma","delta"))
sel <- split(sel, f=as.character(sel$team))
sel <- sel[i]
theta <- unlist(lapply(sel, "[", -1))
theta <- c(theta, omega=tail(omega, 1))
theta

g <- 0:10
## pla <- outer(g, g, function(x, y) Prob(gh=x, ga=y, theta))
pla <- outer(g, g, function(x, y) Prob(gh=x, ga=y, theta, hf=TRUE))
sum(pla) ## must be equal to 1.

res <- c(homeWin=sum(pla[lower.tri(pla)]),
         draw=sum(diag(pla)),
         awayWin=sum(pla[upper.tri(pla)]))

barchart(res~names(res), main=paste(i, collapse=" vs "),
         horizontal=FALSE, stack=FALSE, origin=0,
         xlab="Outcome of the match",
         ylab="Predicted probability")

## The outcome with highest probability.
m <- which(pla==max(pla), arr.ind=TRUE); m

##----------------------------------------------------------------------

## Combines team and round in a new variable.
dg <- dc
dg <- transform(dg,
                hr=paste(home, rod, sep="."),
                ar=paste(away, rod, sep="."))
dd$tr <- with(dd, paste(team, rod, sep="."))

dg <- merge(dg,
            subset(dd, select=c("gamma","delta","tr")),
            by.x="hr", by.y="tr")
names(dg)[ncol(dg)-1:0] <- c("hgamma", "hdelta")

dg <- merge(dg,
            subset(dd, select=c("gamma","delta","tr")),
            by.x="ar", by.y="tr")
names(dg)[ncol(dg)-1:0] <- c("agamma", "adelta")

omg <- data.frame(omega=omega, rod=seq_along(omega))
dg <- merge(dg, omg, by="rod")

dg <- arrange(dg, dat)
str(dg)

##----------------------------------------------------------------------

## HomeWin Draw AwayWin matrix.
hda <- with(dg, {
    d <- h-a
    draw <- d==0
    winH <- d>0
    ## winA <- !draw & !winH
    winA <- !(draw | winH)
    x <- cbind(winH, draw, winA)
    return(1*x)
})

## Estimated probability of HomeWin Draw AwayWin matrix.
g <- 0:10
ephda <- t(apply(dg[,13:17], MARGIN=1,
                 FUN=function(i){
                     pla <- outer(g, g,
                                  function(x, y){
                                      Prob(gh=x, ga=y, i, hf=TRUE)
                                  })
                     res <- c(homeWin=sum(pla[lower.tri(pla)]),
                              draw=sum(diag(pla)),
                              awayWin=sum(pla[upper.tri(pla)]))
                     return(res)
                 }))
## head(ephda)

## Per match DeFinetti Index.
dfi <- rowSums((ephda-hda)^2)

histogram(~dfi|factor(dg$rod),
          as.table=TRUE,
          xlab="DeFinetti Index")

xyplot(dfi~rod, data=dg,
       xlab="Round", ylab="DeFinetti Index")

## Average DeFinetti Index.
mean(dfi)

##----------------------------------------------------------------------
## Saving the results of MAP-1 approach.

map1 <- list(dd=dd, fr=fr, omega=omega, dg=dg, dfi=dfi)

##----------------------------------------------------------------------

## Number of teams and number of rounds.
nteams <- nlevels(db$home)
nrounds <- length(unique(db$rod))
c(nteams, nrounds)

## gamma: teams attacking strength.
## delta: teams deffensive strength.
## omega: home-field advantage parameter.
gamma <- matrix(numeric(), nrow=nteams, ncol=nrounds+1)
delta <- gamma
omega <- numeric(nrounds)
kappa <- omega

## To store values in each match inside a round.
lambdaH <- matrix(numeric(), nrow=nteams/2, ncol=nrounds)
lambdaA <- lambdaH

i <- 1
gamma[,i] <- priors2$mu
delta[,i] <- priors2$theta
omega[i] <- attr(priors2, "hf")$omega$mu0
kappa[i] <- attr(priors2, "hf")$kappa$theta0

## Fitting for each round.
for(i in 1:nrounds){
    ##-------------------------------------------
    ## Subset each round (10 matches).
    ## 
    dbi <- subset(dc, rod==i)
    ##
    ##-------------------------------------------
    ## Model matrices, which teams confronted each other.
    ## 
    Xh <- model.matrix(~-1+home, data=dbi)
    Xa <- model.matrix(~-1+away, data=dbi)
    ##
    ##-------------------------------------------
    ## Columns of Home and Away teams.
    ##
    posh <- apply(apply(Xh, MARGIN=1, FUN="==", 1), MARGIN=2, FUN=which)
    posa <- apply(apply(Xa, MARGIN=1, FUN="==", 1), MARGIN=2, FUN=which)
    ##
    ##-------------------------------------------
    ## Lambda values.
    ##
    lambdaH[,i] <- as.vector(
        exp(Xh%*%gamma[,i]-Xa%*%delta[,i]+omega[i])
        )
    lambdaA[,i] <- as.vector(
        exp(Xa%*%gamma[,i]-Xh%*%delta[,i])
        )
    ##
    ##-------------------------------------------
    ## Goals expected and goals scored.
    ## 
    ## cbind(goalRateHome=lambdaA[,i], goalRateAway=lambdaB[,i],
    ##       goalScorHome=dbi$x, goalScorAway=dbi$y)
    ##
    ##-------------------------------------------
    ## Updating priors for the next round.
    ##
    gamma[posh,i+1] <- 
        Xh%*%gamma[,i]+Xh%*%priors2$sigma2*(dbi$h-lambdaH[,i])
    gamma[posa,i+1] <- 
        Xa%*%gamma[,i]+Xa%*%priors2$sigma2*(dbi$a-lambdaA[,i])
    delta[posh,i+1] <- 
        Xh%*%delta[,i]-Xh%*%priors2$psi2*(dbi$h-lambdaA[,i])
    delta[posa,i+1] <- 
        Xa%*%delta[,i]-Xa%*%priors2$psi2*(dbi$h-lambdaH[,i])
    omega[i+1] <- 
        omega[i]+attr(priors2, "hf")$omega$sigma20*
            sum(dbi$h-lambdaH[,i])
    kappa[i+1] <- 
        kappa[i]-attr(priors2, "hf")$kappa$psi20*
            sum(dbi$a-lambdaA[,i])    
}

gamma <- gamma[,-1]
delta <- delta[,-1]
omega <- omega[-1]
kappa <- kappa[-1]

##----------------------------------------------------------------------

## Series for all teams parameters along the rounds.
dd <- melt(data=as.data.frame(delta), value.name="delta")
dd <- cbind(
    melt(data=cbind(team=levels(dc$home), as.data.frame(gamma)),
         id.vars="team", variable.name="rod", value.name ="gamma"),
    delta=dd$delta)
dd$rod <- as.integer(dd$rod)
str(dd)

sc <- list(x=list(at=c(1,10,20,30,38)))

## Legend for gamma and delta.
key.gd <- list(columns=2,
               title="Parameter", cex.title=1.1,
               type="o", divide=1,
               lines=list(
                   col=ps$superpose.line$col[1:2],
                   pch=ps$superpose.symbol$pch[1:2]),
               text=list(c(expression(gamma), expression(delta))))

## All teams.
xyplot(gamma+delta~rod|team, data=dd, type="l",
       as.table=TRUE,
       scales=sc,
       key=key.gd[c(-4,-5)],
       xlab="Championship round",
       ylab="Parameter value")

## Legend for omega and kappa.
key.ok <- list(columns=2,
               title="Parameter", cex.title=1.1,
               type="o", divide=1,
               lines=list(
                   col=ps$superpose.line$col[1:2],
                   pch=ps$superpose.symbol$pch[1:2]),
               text=list(c(expression(omega), expression(kappa))))

## Field parameters.
xyplot(omega+kappa~seq_along(omega), type="b",
       scales=sc, key=key.ok,
       xlab="Championship round",
       ylab="Field parameters")

##----------------------------------------------------------------------

## Select som teams to a depper view, 2 bests and 2 worsts.
sel <- levels(dd$team)[c(1,5,9,20)]
de <- droplevels(subset(dd, team%in%sel))

## Selected teams.
xyplot(gamma+delta~rod|team,
       data=de, type="o",
       as.table=TRUE,
       scales=sc,
       key=key.gd,
       xlab="Championship round",
       ylab="Parameter value")

key.tm <- list(columns=2,
               title="Teams", cex.title=1.1,
               type="o", divide=1,
               lines=list(
                   col=ps$superpose.line$col[1:4],
                   pch=ps$superpose.symbol$pch[1:4]),
               text=list(levels(de$team)))

xyplot(gamma+delta~rod, outer=TRUE, groups=team,
       data=de, type="o",
       as.table=TRUE,
       scales=sc,
       key=key.tm,
       strip=strip.custom(factor.levels=key.gd$text[[1]]),
       xlab="Championship round",
       ylab="Parameter value")

##----------------------------------------------------------------------
## Rating series.

## Rating along the matches.
dd$rating <- with(dd, gamma+delta)
de$rating <- with(de, gamma+delta)

xyplot(rating~rod|team,
       data=dd, type="o", cex=0.8,
       as.table=TRUE,
       scales=sc,
       xlab="Championship round",
       ylab=expression("Rating"~(gamma+delta)))

xyplot(rating~rod, groups=team,
       data=de, type="o", cex=0.8,
       as.table=TRUE,
       scales=sc,
       key=key.tm,
       strip=strip.custom(factor.levels=key.gd$text[[1]]),
       xlab="Championship round",
       ylab=expression("Rating"~(gamma+delta)))

##----------------------------------------------------------------------

## Average rating.
rtm <- arrange(aggregate(rating~team, data=dd, FUN=mean), -rating)
names(rtm)[2] <- "averating"

## Final rating.
frat <- subset(dd, rod==max(dd$rod))
names(frat)[5] <- "finrating"

## Final, average rating and rank.
fr <- arrange(merge(merge(pl11, frat, by="team"),
                    rtm, by="team"), -finrating)
fr

xyplot(finrating+averating~pts, data=fr,
       auto.key=TRUE, type=c("p","smooth"))

splom(fr[,c("pts","finrating","averating")],
      type=c("p","smooth"))

cor(fr[,c("pts","finrating","averating")],
    method="spearman")

## Como calcular as probabilidades de ganhar, empatar e perder?

##-----------------------------------------------------------------------------
## Who have more chance to gain, Palmeiras or Corinthians?

prob <- function(gh, ga, theta, hf=FALSE){
    ## 1: attack home
    ## 2: defense home
    ## 3: attack away
    ## 4: defense away
    ## 5: home field advantage
    ## 6: away field disadvantage
    if(hf){
        pa <- dpois(gh, lambda=exp(theta[1]-theta[4]+theta[5]))
        pb <- dpois(ga, lambda=exp(theta[3]-theta[2]-theta[6]))
    } else {
        pa <- dpois(gh, lambda=exp(theta[1]-theta[4]))
        pb <- dpois(ga, lambda=exp(theta[3]-theta[2]))
    }
    pa*pb
}
Prob <- Vectorize(prob, c("gh", "ga"))
rm("prob")

i <- c("Avaí", "Corinthians")
i <- rev(i)

sel <- subset(fr, team%in%i,
              select=c("team", "gamma","delta"))
sel <- split(sel, f=as.character(sel$team))
sel <- sel[i]
theta <- unlist(lapply(sel, "[", -1))
theta <- c(theta, omega=tail(omega, 1), kappa=tail(kappa, 1))
theta

g <- 0:10
## pla <- outer(g, g, function(x, y) Prob(gh=x, ga=y, theta))
pla <- outer(g, g, function(x, y) Prob(gh=x, ga=y, theta, hf=TRUE))
sum(pla) ## must be equal to 1.

res <- c(homeWin=sum(pla[lower.tri(pla)]),
         draw=sum(diag(pla)),
         awayWin=sum(pla[upper.tri(pla)]))

barchart(res~names(res), main=paste(i, collapse=" vs "),
         horizontal=FALSE, stack=FALSE, origin=0,
         xlab="Outcome of the match",
         ylab="Predicted probability")

## The outcome with highest probability.
m <- which(pla==max(pla), arr.ind=TRUE); m

##----------------------------------------------------------------------

## Combines team and round in a new variable.
dg <- dc
dg <- transform(dg,
                hr=paste(home, rod, sep="."),
                ar=paste(away, rod, sep="."))
dd$tr <- with(dd, paste(team, rod, sep="."))

dg <- merge(dg,
            subset(dd, select=c("gamma","delta","tr")),
            by.x="hr", by.y="tr")
names(dg)[ncol(dg)-1:0] <- c("hgamma", "hdelta")

dg <- merge(dg,
            subset(dd, select=c("gamma","delta","tr")),
            by.x="ar", by.y="tr")
names(dg)[ncol(dg)-1:0] <- c("agamma", "adelta")

omg <- data.frame(omega=omega, kappa=kappa, rod=seq_along(omega))
dg <- merge(dg, omg, by="rod")

dg <- arrange(dg, dat)
str(dg)

##----------------------------------------------------------------------

## HomeWin Draw AwayWin matrix.
hda <- with(dg, {
    d <- h-a
    draw <- d==0
    winH <- d>0
    ## winA <- !draw & !winH
    winA <- !(draw | winH)
    x <- cbind(winH, draw, winA)
    return(1*x)
})

## Estimated probability of HomeWin Draw AwayWin matrix.
g <- 0:10
ephda <- t(apply(dg[,13:18], MARGIN=1,
                 FUN=function(i){
                     pla <- outer(g, g,
                                  function(x, y){
                                      Prob(gh=x, ga=y, i, hf=TRUE)
                                  })
                     res <- c(homeWin=sum(pla[lower.tri(pla)]),
                              draw=sum(diag(pla)),
                              awayWin=sum(pla[upper.tri(pla)]))
                     return(res)
                 }))
## head(ephda)

## Per match DeFinetti Index.
dfi <- rowSums((ephda-hda)^2)

histogram(~dfi|factor(dg$rod),
          as.table=TRUE,
          xlab="DeFinetti Index")

xyplot(dfi~rod, data=dg,
       xlab="Round", ylab="DeFinetti Index")

## Average DeFinetti Index.
mean(dfi)

##----------------------------------------------------------------------
## Saving the results of MAP-2 approach.

map2 <- list(dd=dd, fr=fr, omega=omega, kappa=kappa, dg=dg, dfi=dfi)

##----------------------------------------------------------------------

## fp: number of field parameters.
map12 <- rbind(cbind(map1$dd, fp=1),
               cbind(map2$dd, fp=2))

## Select som teams to a depper view, 2 bests and 2 worsts.
sel <- levels(map12$team)[c(1,5,9,20)]
de <- droplevels(subset(map12, team%in%sel))

key.map <- list(columns=2,
               title="Model", cex.title=1.1,
               type="o", divide=1,
               lines=list(
                   col=ps$superpose.line$col[1:2],
                   pch=ps$superpose.symbol$pch[1:2]),
               text=list(c("MAP-1","MAP-2")))

sc <- list(x=list(at=c(1,10,20,30,38)),
           y=list(relation="free"))

useOuterStrips(combineLimits(
    xyplot(gamma+delta~rod|team, outer=TRUE, groups=fp,
           data=de, type="o",
           as.table=TRUE,
           scales=sc,
           key=key.map,
           xlab="Championship round",
           ylab="Parameter value")),
    strip.left=strip.custom(factor.levels=key.gd$text[[1]])
    )