Sum to Zero Constraint in GeneralizedLinearModelFit

Question

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.

Answer 1

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