Re: [R] Constrain coefs. in linear model to sum to 0

From: Martin Maechler <maechler_at_stat.math.ethz.ch>
Date: Tue 08 Aug 2006 - 18:23:48 EST

>>>>> "BillV" == <Bill.Venables@csiro.au> >>>>> on Tue, 8 Aug 2006 15:52:20 +1000 writes:

    BillV> Gorjanc Gregor asks:

>> Hello!
>>
>> I would like to use constrain to sum coeficients of a
>> factor to 0 instead of classical corner contraint i.e. I
>> would like to fit a model like
>>
>> lm(y ~ 1 + effectA + effectB)
>>
>> and say get parameters
>>
>> intercept
>> effectA_1
>> effectA_2
>> effectB_1
>> effectB_2
>> effectB_3
>>
>> where effectA_1 represents deviation of level A_1 from intercept and
>> sum(effectA_1, effectA_2) = 0 and the same for factor B.
>>
>> Is this possible to do?

    BillV> Yes, but not quite as simply as you would like. If you set

    BillV> options(contrasts = c("contr.sum", "contr.poly"))

    BillV> for example, then factor models are parametrised as
    BillV> you wish above, BUT you don't get all the effects
    BillV> directly

    BillV> In your case above, for example, if fm is the fitted     BillV> model object, then

    BillV> coef(fm)

    BillV> Will give you intercept, effectA_2, effectB_2,
    BillV> effectB_3.  The remaining effects*_1 you will need to
    BillV> calculate as the negative of the sum of all the
    BillV> others.

    BillV> This gets a bit more complicated when you have     BillV> crossed terms, a*b, but the same principle applies.

Further, note that there have been functions

    dummy.coef()
and

    model.tables()

that have been intended to do exactly this. Unfortunately, these two functions only work correctly in "simpler" cases (e.g., complete balanced designs, IIRC)

E.g. help(model.tables) has

 >> Warning:
 >> 
 >>      The implementation is incomplete, and only the simpler cases have
 >>      been tested thoroughly.

and help(dummy.coef)

  >> Details:
  >> 
  >>   [........]
  >> 
  >>      The method used has some limitations, and will give incomplete
  >>      results for terms such as 'poly(x, 2))'.  However, it is adequate
  >>      for its main purpose, 'aov' models.
  >> 
  >>   [........]
  >> 
  >> Warning:
  >> 
  >>      This function is intended for human inspection of the output: it
  >>      should not be used for calculations.  Use coded variables for all
  >>      calculations.
  >> 
  >>      The results differ from S for singular values, where S can be
  >>      incorrect.

When we last discussed this, IIRC,
we (some within R-core) were waiting for interested (and knowledgable) users to provide improved code for these functions. Hey, if you want to become immortal in the R annals, now is your chance! ;-)

    BillV> Bill Venables.

>>
>> Lep pozdrav / With regards,
>> Gregor Gorjanc
> [.........]

>> ----------------------------------------------------------------------
>> "One must learn by doing the thing; for though you think you know it,
>> you have no certainty until you try." Sophocles ~ 450 B.C.

    BillV> Well, now's your chance!

indeed!

Martin Maechler, ETH Zurich



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Tue Aug 08 18:35:26 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Tue 08 Aug 2006 - 20:17:43 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.