1
$\begingroup$

Is there a way to impose a constraint on a generalized linear model fit in Mathematica? In R, when using the glm() function, you can set options(contrasts=c('YY.sum', 'ZZ.sum')). Is there something like that in Mathematica?

GeneralizedLinearModelFit[{X, response}, ExponentialFamily -> "Poisson", "INSERT SOME CONSTRAINT HERE"]

where X is the design matrix and response is the response matrix. Or possible manipulate the data to achieve this constraint? A "corner" constraint would be equivalent but would change the interpretation of the explanatory variables and I don't want to do this.

My regression model is not fully defined (rank of the design matrix is too low) unless a sum to zero constraint is implemented.

Thank you.

$\endgroup$
  • $\begingroup$ I am recreating D Karlis' paper,tolstoy.newcastle.edu.au/R/e8/help/att-6544/karlisntzuofras03.pdf. He mentions a sum to zero constraint. So I am modeling two years of France soccer data (league 1 and league 2) so there are 42 teams each with an attack and defense parameter. $\endgroup$ – Ray Troy Apr 13 '16 at 1:09
4
$\begingroup$

I'm not aware any coding options for such constraints but you can roll your own. Here's an example with using two "treatments" with different Poisson means.

(* Generate data from two Poisson distributions with very different sample sizes *)
n1 = 1000;
n2 = 10;
λ1 = 5;
λ2 = 15;
SeedRandom[1234];
v1 = RandomVariate[PoissonDistribution[λ1], n1];
v2 = RandomVariate[PoissonDistribution[λ2], n2];
v = Flatten[{v1, v2}];

(* So if we code things succesfully, we should get two parameter estimates
with approximately the following values:  Overall mean of the logs and
deviation of Log[λ1] from mean of logs *)
N[{(Log[λ1] + Log[λ2])/2, Log[λ1] - (Log[λ1] + Log[λ2])/2}]
(* {2.15874,-0.549306} *)

(* Design matrix for sum-to-zero constraints *)
m = Transpose[{Table[1, {i, n1 + n2}],
    Flatten[{Table[1, {i, n1}], Table[-1, {i, n2}]}]}];

(* Fit model *)
glm = GeneralizedLinearModelFit[{m, v}, ExponentialFamily -> "Poisson"];
estimates = glm["BestFitParameters"]
(* {2.19727,-0.562745} *)

Update

Suppose you have lots of factor levels and you want estimates of the mean of the factor effects and the individual deviations from that overall mean (i.e., the sum to zero constraint). You can just fit with a convenient set of dummy variables that fits the same overall model and then determine the sum-to-zero characterization afterwards. What I show below gets you the estimates but you should at some point want standard errors for the estimates. That's not difficult to get, too, but here I'm only showing how to get the estimates. (If one can be more specific about the models, then it's easier to write a more specific example.)

(* Generate some data from several Poisson distributions *)
n = {10, 5, 15, 12, 11}; (* Sample sizes *)
λ = {5, 15, 10, 9, 20};  (* Poisson means *)
SeedRandom[1234];
(* Observed counts *)
y = Flatten[
   Table[RandomVariate[PoissonDistribution[λ[[i]]], 
     n[[i]]], {i, Length[n]}]];
(* id is position of Poisson mean *)
id = Flatten[Table[Table[i, {j, n[[i]]}], {i, Length[n]}]];


(* Fit a "cell means" model where each unique value of the nominal
   variable results in estimate of the associated parameter *)
glm = GeneralizedLinearModelFit[Transpose[{id, y}], x, x, 
   NominalVariables -> x,
   IncludeConstantBasis -> False, ExponentialFamily -> "Poisson"];
mle = glm["BestFitParameters"]

(* Construct parameter estimates characterized by the sum-to-zero constraint *)
(* Overall mean *)
λ0 = Mean[mle]
(* Deviations from mean *)
λDeviations = mle - λ0
$\endgroup$
  • $\begingroup$ Sorry for the delay, was on holiday. Thank you. Is there a way to generalize it so if I have 90+ predictors it implements this constraint? $\endgroup$ – Ray Troy Apr 12 '16 at 0:43
  • $\begingroup$ 90+??? You're funny! I guess I would first try to talk you out of that as maybe there could be one of the categories that could be considered the "base level" ? Please! But, yes, a large number could be done. I'll put together a larger example (but not much increase in code). $\endgroup$ – JimB Apr 12 '16 at 1:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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