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 enter image description here

is the regression.

My data is here, https://drive.google.com/drive/u/0/folders/0B6WPzNPrW5SfMkxmbnRRX2FvY0k

It is two years of the France leagues.

rawData = Import[filePathImport];
header = First[rawData];
data = Riffle[Rest@rawData, Rest@rawData];
dataAssoc = AssociationThread[header, #] & /@ data;
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.

  • $\begingroup$ 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. $\endgroup$
    – MarcoB
    Apr 19, 2016 at 1:50
  • $\begingroup$ 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. $\endgroup$
    – Ray Troy
    Apr 19, 2016 at 2:16
  • $\begingroup$ 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. $\endgroup$
    – MarcoB
    Apr 25, 2016 at 3:45
  • $\begingroup$ @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. $\endgroup$
    – Ray Troy
    Apr 25, 2016 at 23:27
  • $\begingroup$ I'm glad to hear that. This is a potentially fascinating question; its formulation just needs some more help. $\endgroup$
    – MarcoB
    Apr 26, 2016 at 5:28

1 Answer 1


Using the data from the link,

This is only for model 2 of his paper.

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:

enter image description here

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?

  • $\begingroup$ @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. $\endgroup$
    – Ray Troy
    Apr 27, 2016 at 20:02
  • $\begingroup$ Hi Ray, have you found a solution for your issue? I would be really interested in it. I am trying to implement something similar too: basically, I applied simple poisson formula for home and for away teams for each match by considering the average of, respectively, home scored and away experienced goals (home) and home experienced and away scored goals (away). I would like to add other indicators like home factor, healt status and shoot attempts by implement numerical simulations but I cannot do it. Anyway, I think you'll find useful this: community.wolfram.com/groups/-/m/t/517631 $\endgroup$
    – Nate
    Feb 16, 2020 at 8:27
  • $\begingroup$ @Mick, I did find a solution but it wasn't using GeneralizedLinearModelFit. I submitted a feature request to Wolfram over 3 years ago and the sum to zero constraint is still not implemented. I just let my WolframOne subscription lapse as I am done with the Wolfram Language. If you direct message me (is that even possible?) I can send you over the notebook (s). I built many models and also tested those models vs the sportsbook odds using the Kelly Criterion. After 3 years on the project, gave it up as there was no edge to be had. I am would collaborate if you want to revive the project. $\endgroup$
    – Ray Troy
    Feb 20, 2020 at 2:13
  • $\begingroup$ Ray, I would like to direct message you but I think it is forbidden here. I would really appreciate to see your notebook because of my personal interest in this issue as well as to learn more on the wolfram programming environment (I am beginner). I see your website but it seems not working now, so tried to write you on facebook page (isn't Ray Troy Consulting - like the website - yours ?) $\endgroup$
    – Nate
    Feb 24, 2020 at 11:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.