Football predictions display

Having looked at the football data earlier, I wanted to look at predictions for new games. This consists of two parts, getting a predictive model, predicting and displaying the predictions. I decided to do this backwards, first to make the displays. This will make things easier when the time is there to compare models. To get the predictions I use a very simple model, which basically states, a club makes about x goals, irrespective of all other conditions. I don't believe this model, but it can give predictions.
model1 <- glm(Goals ~OffenseClub,data=StartData,family='poisson')
The consequence of this setup is that each game needs two predictions, one for the first club, one for the second. For clubs Vitesse and FC Groningen are used to make the predictions.

The prediction of the glm is a mean number of goals, which is still quite far from the reality of a number of goals. For this I use the Poisson distribution and treat the prediction as true. I do not include overdispersion nor standard error of parameters. The result shows FC Groningen has 30% of getting no goals, 35% chance of getting 1 goal, 22 % for two goals, after which the chances become quickly very low. 
top <- data.frame(OffenseClub=c('FC Groningen','Vitesse'))
prepred <- predict(model1,top,type='response')
(dp1 <- dpois(0:5,prepred[1]))[1:4]
[1] 0.29942768 0.36107456 0.21770672 0.08750956
Finally, the predictions need to be combined, to get a pair of goals. Not surprisingly, if the chance of a particular outcome, such as one goal, is 30%, then the chance of a pair of outcomes, such as 1-1 may be 30%*30%=9%. In this case it turns out to be slightly higher, 12%. This is the most probable outcome too. 
dp2 <- dpois(0:6,prepred[2])
oo <- outer(dp1,dp2)
rownames(oo) <- 0:6
colnames(oo) <- 0:6
round(oo,digits=3)

      0     1     2     3     4     5     6
0 0.073 0.103 0.073 0.034 0.012 0.003 0.001
1 0.088 0.124 0.088 0.041 0.015 0.004 0.001
2 0.053 0.075 0.053 0.025 0.009 0.002 0.001
3 0.021 0.030 0.021 0.010 0.004 0.001 0.000
4 0.006 0.009 0.006 0.003 0.001 0.000 0.000
5 0.002 0.002 0.002 0.001 0.000 0.000 0.000
6 0.000 0.000 0.000 0.000 0.000 0.000 0.000

It will be more useful to summarize the outcome as win-equal-lost. These are extracted as sums of probabilities. 

c(sum(oo[upper.tri(oo)]),sum(diag(oo)),sum(oo[lower.tri(oo)]))
[1] 0.4167411 0.2612236 0.3211237
It is practical to fit all this in a little function which creates these data in one go. The only new things are the introduction of a new class fboo which is used to direct the prediction to the appropriate accompanying print function and some attributes to administrate the clubs predicted.
fbpredict <- function(object,club1,club2) {
  top <- data.frame(OffenseClub=c(club1,club2),DefenseClub=c(club2,club1),OffThuis=c(1,0))
  prepred <- predict(object,top,type='response')
  dp1 <- dpois(0:9,prepred[1])
  dp2 <- dpois(0:9,prepred[2])
  oo <- outer(dp2,dp1)
  rownames(oo) <- 0:9
  colnames(oo) <- 0:9
  class(oo) <- c('fboo',class(oo))
  attr(oo,'row') <- club1
  attr(oo,'col') <- club2
  wel <- c(sum(oo[upper.tri(oo)]),sum(diag(oo)),sum(oo[lower.tri(oo)]))
  names(wel) <- c(club1,'equal',club2)
  return(list(details=oo,'summary chances'=wel))
}

print.fboo <- function(x,...) {
  cat(attr(x,'row'),'in rows against',attr(x,'col'),'in columns \n')
  class(x) <- class(x)[-1]
  attr(x,'row') <- NULL
  attr(x,'col') <- NULL
  oo <- formatC(x,format='f',width=4) # fixed format
  oo <- gsub('\\.0+$','       ',oo)   # replace trailing 0 by ' '
  oo <- substr(oo,1,6)                # and fix the width
  print(oo,quote=FALSE,justify='left')
}

fbpredict(model1,'FC Groningen','Vitesse')

$details

FC Groningen in rows against Vitesse in columns 
  0      1      2      3      4      5      6      7      8      9     
0 0.0730 0.0880 0.0531 0.0213 0.0064 0.0016 0.0003 0.0001 0      0     
1 0.1030 0.1242 0.0749 0.0301 0.0091 0.0022 0.0004 0.0001 0      0     
2 0.0727 0.0877 0.0529 0.0213 0.0064 0.0015 0.0003 0.0001 0      0     
3 0.0342 0.0413 0.0249 0.0100 0.0030 0.0007 0.0001 0      0      0     
4 0.0121 0.0146 0.0088 0.0035 0.0011 0.0003 0.0001 0      0      0     
5 0.0034 0.0041 0.0025 0.0010 0.0003 0.0001 0      0      0      0     
6 0.0008 0.0010 0.0006 0.0002 0.0001 0      0      0      0      0     
7 0.0002 0.0002 0.0001 0      0      0      0      0      0      0     
8 0      0      0      0      0      0      0      0      0      0     
9 0      0      0      0      0      0      0      0      0      0     

$`summary chances`
FC Groningen        equal      Vitesse 
   0.3213815    0.2612237    0.4173918 

Football (Eredivisie) goals

The football season has started in Netherlands, so I went and had a look at last year's scores. I did not find downloadable data, at http://www.eredivisiestats.nl/wedstrijden.php I could copy last season's data and paste into a spreadsheet. The default layout did not look good to me for statistical analysis (first lines shown below), so first part is reformatting to get how many goals a club did against another. This means that each game becomes two lines and an indicator variable is used for who was playing at home.

Datum	Thuisclub (home)	Uitclub (away)	Thuisscore	Uitscore
2011-08-05	Excelsior	Feyenoord	0	2
2011-08-06	RKC Waalwijk	Heracles Almelo	2	2
2011-08-06	Roda JC	FC Groningen	2	1

# read data
old <- read.csv('seizoen1112.csv',stringsAsFactors=FALSE)
# reformat lengthwise
StartData <- with(old,data.frame(
OffenseClub = c(Thuisclub,Uitclub),
DefenseClub = c(Uitclub,Thuisclub),
Goals = c(Thuisscore,Uitscore),
OffThuis = rep(c(1,0),each=length(Datum))
))
#remove space as first character in names
levels(StartData$OffenseClub) <- sub('^ ','',levels(StartData$OffenseClub))
levels(StartData$DefenseClub) <- sub('^ ','',levels(StartData$DefenseClub))
# check clubs have same name and ordening in both factors.
all(levels(StartData$OffenseClub)==levels(StartData$DefenseClub))

Distribution of goals

To get to know the data the distribution of goals is most interesting:

xt <- xtabs(~Goals ,data=StartData)

plot(xt)

The first question is; is this approximately Poisson distributed? MASS has some functions.

library(MASS)

(fd <- fitdistr(StartData$Goals, densfun="Poisson"))

     lambda  

  1.62908497 

 (0.05159364)

round(dpois(0:8,coef(fd))*sum(xt))

plot(xt,ylim=c(0,200),ylab='Frequency of goals')

points(x=0:8,y=dpois(0:8,coef(fd))*sum(xt),col='green')

It is also possible to do the same calculations with the fitdistrplus package and get some extra information. Two different algorithms give essentially the same information.

library(fitdistrplus)

(fd2 <- fitdist(StartData$Goals, distr='pois',method='mle'))

Fitting of the distribution ' pois ' by maximum likelihood 

Parameters:

       estimate Std. Error

lambda 1.629085 0.05159362

fitdist(StartData$Goals, distr='pois',method='mme')

Fitting of the distribution ' pois ' by matching moments 

Parameters:

       estimate

lambda 1.629085

plot(fd2)

This is nice. I get the same information as my own plot, but also a CDF and have to do less. If only I knew all packages in R, I would not have to program anything. But there are more extra's in fitdistrplus. Goodness of fit and bootstrap confidence interval.

gofstat(fd2,print.test=TRUE)

Chi-squared statistic:  18.46365 

Degree of freedom of the Chi-squared distribution:  4 

Chi-squared p-value:  0.001001433 

bd <- bootdist(fd2)

summary(bd)

Parametric bootstrap medians and 95% percentile CI 

  Median     2.5%    97.5% 

1.627451 1.527803 1.728736 

plot(density(bd$estim$lambda))

Finally, if desired the same can be done the Bayesian way via Jags

library(R2jags)

calcData <- with(StartData,list(N=length(Goals),Goals=Goals))

mijnmodel <- function() {

  for (i in 1:N) {

    Goals[i] ~ dpois(mu)

  }

  mu ~ dunif(0,3)

}

write.model(mijnmodel,'model.file')

jags.inits <- function() {

  list(

      mu = runif(1,0,3)

  )

}

jags.params=c('mu')

jags(data=calcData, inits=jags.inits, jags.params,model.file='model.file')

Inference for Bugs model at "model.file", fit using jags,

 3 chains, each with 2000 iterations (first 1000 discarded)

 n.sims = 3000 iterations saved

         mu.vect sd.vect     2.5%      25%      50%      75%    97.5%  Rhat n.eff

mu         1.633   0.051    1.533    1.598    1.632    1.667    1.732 1.001  2100

devianc 2048.287   1.338 2047.308 2047.403 2047.757 2048.671 2052.118 1.006  1500

For each parameter, n.eff is a crude measure of effective sample size,

and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)

pD = 0.9 and DIC = 2049.2

DIC is an estimate of expected predictive error (lower deviance is better).

This leads to essentially the same result as before.

Conclusion

Three ways to get parameters of a distribution are used. The most obvious way, via MASS gives parameters and standard deviation. A dedicated package gives all that plus significance and distribution. The Bayesian way gives the distribution too and uses most code by far. The estimates do not differ in any relevant way.

Random and fixed effects in sensory profiling

I am reading Introduction into mixed modelling by N.W. Galway. It is partly a repeat of things I know, but I expect to use mixed models quite a lot the coming time, so it is good to repeat these things.
My problem with this book is a sensory example in chapter 2. It is profiling data, but with some twists, as happens in real life. He also makes some odd choices.

Design

The test concerns four hot products (ravioli). Because the products are hot, it was found not feasible to use a proper balanced design within each day (e.g. Williams designs). Within a day each presentation (he used presentation for rounds) has one product, as shown in the table below. Which is fine in a way, I do believe in making designs which be executed in practice. However, I come from a facility where we would at least have tried to solve this with the 'au-bain-marie'.

Table 1, allocation of products A to D over presentations and days.

Day	Presentation 1	Presentation 2	Presentation 3	Presentation 4
1	B	A	C	D
2	C	D	B	A
3	A	C	B	D

Within each presentation the nine assessors get served portions. It was randomized though not registered who got what serving. Even so, a random variable was created to serve as proxy. I understand this does not make a difference if servings are nested in presentations and days. Still, it removes the option of actually looking at this effect except as nested within servings and days. I can actually imagine that the product served first is different, especially in case of sloppy sensory practices, so it would be nice to be able to check this.

Model

The following effects were used in the model(s)
Table 2. Allocation of effects

Effect	Allocation
Day	Random
Day.Presentation	Random
Day.Presentation.Serving	Random
Product	Fixed
Assessor	Fixed
Product.Assessor	Fixed

Here I got my ax to grind.

• Product.Assessor. To quote: '(assessor) ANA perceived Brand B as the least salty, whereas (assessor) GUI perceived Brand  A as the least salty. Such crossover effects may be important: they suggest an obstacle to designing a brand that will please all consumers'. In my (industry) experience the sensory panel may well be in a different country or even continent as the target market, so that is a bit too strong. Besides, nine assessors is a bit low to segment groups on. Typically segmenting would be done with a consumer group, say 120 persons, in the target market. In my view  sensory is intended to provide an objective measurement. The assessors are representing what a typical human may taste, and are hence random. With assessors random, product.assessor can only be found random too.
• Day.Presentation. To quote: 'Presentation 1 on Day 1 is not the the same as presentation 1 on day 2' and 'it would therefore not be meaningful to obtain the mean for presentation 1 over days'. Actually, presentation 1 is the same on day 1 as on day 2 . It is the tasting with a clean palette. While palette cleansing is supposed to make all presentations equal, that does not mean it really does. After all, this is on the trade-off between production (shorter cleaning time) and correctness (longer cleaning time). Looking at this effect is important. I would make it fixed. For round 1 specifically, it represents monadic tasting as in a consumer test. Surely a first impression is important to understand consumer liking, when correlating consumer data and sensory data I might look at solely presentation 1. Having said that, given the design in table 1, I would make presentation random. The information is not there to do anything else. 

Conclusion

I dislike the design, the sensory foundation and interpretation. However, if you ignore that sensory based dislike, the model actually makes sense and is a nice example. On top of that, the book has R code using lme (nlme package), with code on the download site http://www.wiley.com/legacy/wileychi/mixed-modelling/ has both nlme and lme4 examples, so is up to date. I am looking forward to reading the next chapters.

Trying Julia

In my previous post I tried building Williams designs in R. Since that code was running a bit slow, this was an ideal test for Julia. Big enough to be at least slightly realistic, small enough that it is doable.
I am very impressed. Almost twenty fold speed increase, even though this was the best I could do in R, the most naive way possible in Julia.
                                       R              Julia        Ratio
Double Williams design 4, 100 times    3.97 sec       0.212 sec    0.053   
Williams design 5, 10 times          289.5 sec       18.54  sec    0.064 


I'd surely love to use Julia more. For instance, if I could port some of the algorithm's of Professor Ng Machine Learning class to Julia? I don't think it is possible yet, but what is not, may come.

Julia code

macro timeit(ex,name,num)

    quote

        t0 = 0

        for i=1:$num

            t0 = t0 + @elapsed $ex

        end

        println("julia,", $name, ",", t0)

        #gc()

    end

end

function gendesign(ncol)

  nrow=ncol*2

  desmat = zeros(Uint8,nrow,ncol)

  desmat[1,1:ncol] = [1:ncol]

  for i = [1:ncol]

    desmat[2*i-1,1] = i

    desmat[2*i,1]=i

  end

  carover = zeros(Uint8,ncol,ncol)

  for i = 1:(ncol-1)

    carover[i,i+1] = 1

  end

  count   = 0

  addpoint( desmat , carover,count)

end

function numzero(matin) 

  length(matin) - nnz(matin) 

end   

 

function first0(matin)

   if numzero(matin)==0

      return -1

   end   

   nrow, ncol = size(matin)

   for row = 1:nrow

      for col = 1:ncol

         if matin[row,col] == 0

            return row, col

         end

      end

    end   

end   

function addpoint(desmat,carover,count)

  if nnz(desmat) == length(desmat)

     count +=1

     print("x")

     return count

  end

  row,col  = first0(desmat)

  for i = [1:size(desmat,2)]

     if numzero(desmat[row,:]-i) == 0

        if numzero(desmat[:,col]-i) < 2

           if carover[desmat[row,col-1],i] < 2 

              if (col !=2) | (desmat[row,1] != desmat[row-1,1]) | (desmat[row-1,col] < i)

                 desmat[row,col]=i

                 carover[desmat[row,col-1],i] +=1

                 count = addpoint(desmat,carover,count)

                 desmat[row,col]=0

                 carover[desmat[row,col-1],i] -=1

              end

            end 

         end 

      end 

   end 

   return count

end

@timeit gendesign(4) "design 4 " 100

@timeit gendesign(5) "design 5 " 10

Note

Both designs were created using the same script. If you ask an even number of design points using the odd algorithm, then this will work. You just get a design with double the rows, carry over balanced, every treatment equally often in each row and each column.

Williams designs with 5 products

In a previous post I created small Williams designs for an even number of products. This worked very well, also because the number of permutations could be restricted significantly due to symmetry. Unfortunately this does not work so well with an odd number of products. The symmetry is not so obvious. In fact, using brute force I can find non-symmetric solutions for five products. However, I am lacking force to handle more than five products in a reasonable time.

design setup

The design fixed parameters

To get the design I used a few assumptions. The number of rows is twice the number of columns. This is because I remember reading there are no solutions with the same number of rows as columns. The first row reads 1:n. The first column is all elements of 1:n twice.

Filling in the design

The approach is actually very simple.

• Look at the first not yet know point. 
• Examine which value it could have. 
• not yet used in current row
• at most once used in current column
• at maximum once used after the product in current row, one column to the left (carry over)
• if this is the second column and the first column of this row is the same as the first column of the previous row, then the current value must be higher than the value of one row up.
• Try the possible values in a loop
• Within that loop go back to step one

Implemented recursively. In all actually, it is quite brute force. I should be doing this in C or Julia rather than in R which is relatively slow. 

Tricks

Compared to the version with an even number of products I implemented some things a bit smarter. Removed some if within a for loop by selecting the objects over which the loop sequences. I also noticed some things which surprised me in making things faster. For instance, creating temporary variables (e.g. row <- desconst$todo[desobject$count,2L]) is faster than using desconst$todo[desobject$count,2L] to index matrices later on. Merging two simple lines of code into one complex line of code is slower. Apparently it is cheaper to store a small result than have the interpreter figure out a complex line of code. Hence I used quite a lot of local variables, certainly more than I find 'nice'. 

Designs

The result if 90 Williams designs for five products. Without doubt there is quite some duplication from permutations of the objects. The designs can be separated in three classes. Designs symmetric in columns (72 of those), designs with four rows symmetric in columns (five of those) and designs where are not symmetric in columns at all (13 of those).

design symmetric in columns

The first design has five unique rows each copied reversed to make the ten rows. For ease I colored the matching rows.

      [,1] [,2] [,3] [,4] [,5]

 [1,]    1    2    3    4    5

 [2,]    1    3    2    5    4

 [3,]    2    1    4    3    5

 [4,]    2    4    1    5    3

 [5,]    3    1    5    2    4

 [6,]    3    5    1    4    2

 [7,]    4    2    5    1    3

 [8,]    4    5    2    3    1

 [9,]    5    3    4    1    2

[10,]    5    4    3    2    1

designs four rows symmetric in columns

The colored rows have a reversed row in the same color. 

      [,1] [,2] [,3] [,4] [,5]
 [1,]    1    2    3    4    5
 [2,]    1    3    4    2    5
 [3,]    2    3    5    1    4
 [4,]    2    4    1    5    3
 [5,]    3    1    5    4    2
 [6,]    3    5    2    1    4
 [7,]    4    1    2    5    3
 [8,]    4    5    1    3    2
 [9,]    5    2    4    3    1
[10,]    5    4    3    2    1

design not symmetric in columns

      [,1] [,2] [,3] [,4] [,5]

 [1,]    1    2    3    4    5

 [2,]    1    3    5    2    4

 [3,]    2    4    1    5    3

 [4,]    2    5    3    4    1

 [5,]    3    1    4    2    5

 [6,]    3    2    1    5    4

 [7,]    4    3    5    1    2

 [8,]    4    5    2    1    3

 [9,]    5    1    4    3    2

[10,]    5    4    2    3    1

R code

generation

gendesign <- function(n=3) {

nc <- as.integer(n)

nr <- nc+nc

desmat <- matrix(NA,nrow=nr,ncol=nc)

desmat[1,] <- 1L:nc

desmat[,1] <- rep(1L:nc,each=2)

carover <- matrix(0L,nrow=nr,ncol=nc)

for (i in 1L:(nc-1L)) carover[i,i+1]  <- 1L

todo <- which(is.na(t(desmat)),arr.ind=TRUE)

desobject <- list(desmat=desmat,carover = carover,count=0)

desconst <- list(todo=todo,ntodo=nrow(todo),nc=nc,totest = 1L:nc)

desresult <- list()

addpoint(desobject,desresult,desconst)

}

addpoint <- function(desobject,desresult,desconst) {

desobject$count <- desobject$count+1L

if (desobject$count > desconst$ntodo) {

l <- length(desresult)

desresult[[l+1]] <- desobject$desmat

return(desresult)

}

row <- desconst$todo[desobject$count,2L]

col <- desconst$todo[desobject$count,1L]

dob <- desobject

currow <- desobject$desmat[row,]

currow <- currow[!is.na(currow)]

totest <- desconst$totest[-currow]

prev <- desobject$desmat[row,col-1L]

totest <- totest[desobject$carover[prev,totest]<2L]

counts <- rowSums(outer(totest,desobject$desmat[,col],function(x,y) as.integer(x==y)),na.rm=TRUE)

totest <- totest[counts<2L]

if (col==2 & desobject$desmat[row,1]==desobject$desmat[row-1L,1L]) {

totest <- totest[totest > (desobject$desmat[row-1L,2])]

}

for (i in totest) {

desobject$carover[prev,i] <- desobject$carover[prev,i] + 1L

desobject$desmat[row,col] <- i

desresult <- addpoint(desobject,desresult,desconst)

desobject <- dob

}

desresult

}

g5 <- gendesign(5)

checking

# check that each column has each object twice

sapply(g5,function(sub) all(sapply(1:5,function(x)   

      all(table(sub[,x])==rep(2,5)))))

# check that each product is followed by each other product twice

sub <- g5[[1]]

target <- Reduce('+',lapply(1:4,function(x) table(sub[,x],sub[,x+1])))

target  # shows a symmetric matrix with 0 on diagonal, 2 elsewhere.

sapply(g5,function(sub) identical(target,

      Reduce('+',lapply(1:4,function(x) table(sub[,x],sub[,x+1])))))

#check the number of reversed rows

sa <- sapply(g5,function(sub1) {

10-sum(1:10 %in% sapply(1:10,function(i) {

r1 <- rev(sub1[i,])

mat <- matrix(rep(r1,each=10),nrow=10,ncol=5)

which(rowSums(sub1==mat)==5)

}))

})

table(sa)

# number of reversed rows

# 0  4 10 

#13  5 72 

prints

g5[[1]]

g5[[4]]

g5[[8]]