Skip to main content
Mathematica

Return to Answer

Added an example using more factor levels as requested by the OP.
Source Link
JimB
JimB
  • 42.9k
  • 3
  • 51
  • 108

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

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
Source Link
JimB
JimB
  • 42.9k
  • 3
  • 51
  • 108

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} *)