# Modeling of soccer results using GeneralizedLinearModelFit

I am trying to recreate the regression in "Analysis of sports data by using bivariate Poisson models," which can be found here, http://tolstoy.newcastle.edu.au/R/e8/help/att-6544/karlisntzuofras03.pdf.

It states

is the regression.

It is two years of the France leagues.

rawData = Import[filePathImport];
data = Riffle[Rest@rawData, Rest@rawData];
formatDataAssoc = Quiet@MapAt[DateList, #, Key["date"]] & /@ dataAssoc;
formatMatrix[row_, number_] := Module[{rowNumber},rowNumber =number[[1]]; If[EvenQ[rowNumber],rowFormatHome[row],rowFormatAway[row]]]
rowFormatAway[dataRow_] := {toExpression[dataRow[["away"]], "A"], toExpression[dataRow[["home"]], "D"]
rowFormatHome[dataRow_] := {toExpression[dataRow[["home"]], "A"], toExpression[dataRow[["away"]], "D"], home}
toExpression[one_, two_] := ToExpression[one <> two]
design = MapIndexed[formatMatrix, formatDataAssoc]
teams = Union[Flatten[design]]
DesignMatrix[design, teams, teams, NominalVariables -> All]


I then planned to use GeneralizedLinearModelFit[{designmatrix,reponse},teams, teams]

where the response is

response = Flatten[formatDataAssoc[[All, {"aScore", "hScore"}]] // Values]


each row has the predictors and am not sure why it is not working. I have very little experience with regressions.

In the paper it also states "To achieve identifiability of the above model parameters, we may use any standard set of constraints. Here we propose to use either sum-to-zero or corner constraints, depending on the interpretation that we prefer."

I do not know how to include this constraint in my code.

It might not be best to use the DesignMatrix function to create the data to use in the GeneralizedLinearModelFit.

Thank you.

• Ray, could you explain the terms used in the model? What are lambda, mu, att, etc? What should your model be able to predict? Could you explain a bit more what you intended to do with your code? I am not sure that I understand, for instance, what data = Riffle[Rest@rawData, Rest@rawData]; should do, since it seems to me that it simply doubles the raw data points. – MarcoB Apr 19 '16 at 1:50
• Some of these terms I use might not be correct as regression analysis is new to me. A quick read of section 4 of the paper might help in explanation. @MarcoB the lambda is the correlation coefficient, att is the attack coefficients of each team. These are also the explanatory variables (X1,X2...). def is the defense of each team again, explanatory variables. The riffle is to turn the data from two rows containing data into a single vector. So there are two rows for every game. With one row I use that for lamda1 and the next row for lamda2. This might not be the correct approach. – Ray Troy Apr 19 '16 at 2:16
• Ray, as you can see the question has not attracted much attention so far. This may be because you expect potential responders to read the whole paper and understand its notation, correlate it to your data, then provide a solution. This is unlikely. Perhaps some efforts on your part to clarify the question would help: e.g pointing specifically at the relevant section in the paper, explaining exactly what you want to model and predict, explaining the meaning of the terms in the model, and possibly sharing any more code you have produced in the meantime. Otherwise this is work for hire. – MarcoB Apr 25 '16 at 3:45
• @MarcoB Thank you for your advice. I did not mean to mean it to be "work for hire." I will allow this to expire and then post a better formulated question. – Ray Troy Apr 25 '16 at 23:27
• I'm glad to hear that. This is a potentially fascinating question; its formulation just needs some more help. – MarcoB Apr 26 '16 at 5:28

Using the data from the link,

This is only for model 2 of his paper.

{FileNameSetter[Dynamic[traingingSetDirectory]],Dynamic[traingingSetDirectory]}
parsedMatchesTraining = Import[traingingSetDirectory]
teamToVarFunc[team_] := <|ToString@team -> <|"attack" -> ToExpression[team <> "A"], "defense" -> ToExpression[team <> "D"], "betaA" -> ToExpression[team <> "BA"], "betaH" -> ToExpression[team <> "BH"]|>|>


I couldn't get this to format correctly so I took an image:

setup[var_, teamesToVars_, mu_, hfieldConst_, lamdaConst_, indicator_] := Module[{awayTeamName, homeTeamName, homeGoals, awayGoals, homeAttack, awayAttack, homeDefense, awayDefense, awayParm, homeParm, awayAttackAlpa, homeAttackAlpha, awayDefenseAlpha, homeDefenseAlpha, awayLamda, homeLamda, lamdaFinal}, awayTeamName = Keys[var][[1]]; homeTeamName = Keys[var][[2]]; awayGoals = var[[1, "goalF"]]; homeGoals = var[[2, "goalF"]]; awayAttack = teamesToVars[[awayTeamName, "attack"]]; homeAttack = teamesToVars[[homeTeamName, "attack"]]; awayDefense = teamesToVars[[awayTeamName, "defense"]]; homeDefense = teamesToVars[[homeTeamName, "defense"]]; awayParm = Exp[mu + awayAttack + homeDefense]; homeParm = Exp[mu + homeAttack + awayDefense + hfieldConst]; awayLamda = teamesToVars[[awayTeamName, "betaA"]]; homeLamda = teamesToVars[[homeTeamName, "betaH"]]; lamdaFinal = Exp[Which[indicator == 2, lamdaConst, indicator == 3, lamdaConst + awayLamda, indicator == 4, lamdaConst + homeLamda, indicator == 5, lamdaConst + homeLamda + awayLamda]]; -bivariatePoissonLogLike[awayGoals, homeGoals, awayParm, homeParm, lamdaFinal]]

hfieldConst = hField; lamdaConst = lamda; muConst = mu; teamesToVars = <|teamToVarFunc /@ teamNames|>;
varsToSolveFor = Flatten[{teamesToVars[[All, {"attack", "defense"}]] // Values // Values, muConst, hfieldConst, lamdaConst}];

attackConstriant = Total[teamesToVars[[All, "attack"]]]; defenseConstraint = Total[teamesToVars[[All, "defense"]]];
summed = Total[setup[#, teamesToVars, muConst, hfieldConst, lamdaConst, 2] & /@ parsedMatchesTraining];
solvedParamsLamda = NMinimize[{summed, attackConstriant == 0, defenseConstraint == 0}, varsToSolveFor, MaxIterations -> 10^100]


It is not a direct answer to the question but it answers the larger problem of recreating the paper in Mathematica. In his work he used a Generalized Linear Fit to get initial values to feed into his expectation maximization algorithm. I wanted to recreate his methodology exactly but I do not think it is possible to use GeneralizedLinearModelFit for this problem, @Jimbaldwin can explain more.

One may notice this code is not as one would expect. You could solve the problem using a matrix focused approach and @JimBaldwin has that solution.

With some datasets, you get funny numbers and I am not sure why.

Formatting this is difficult. Why can you not just copy past from Mathematica into Stack Exchange?

• @MarcoB here is the solution. Not pretty but it works. I would love to figure out how to use GenealizedLinearModelFit in this problem and develop my own EM algo just as the author did. If that is posted, it would really be something. My next stop is state space time series modeling of soccer, by S.J. Koopman. – Ray Troy Apr 27 '16 at 20:02