[R] creating 'all' sum contrasts

[Jonathan Christensen]

Michael,

Let c_1 and c_2 be vectors representing contrasts. Then c_1 and c_2
are orthogonal if and only if the inner product is 0. In your example,
you have vectors (1,0,-1) and (0,1,-1). The inner product is 1, so
they are not orthogonal. It's impossible to have more orthogonal
contrasts than you have levels in your factor, a result from basic
linear algebra.

You can get all possible pairwise contrasts, which is different from
orthogonal contrasts (in fact, it's only possible to have floor(n/2)
orthogonal pairwise contrasts). This is probably not the easiest way,
but it works:

n <- 10
M <- matrix(0,nrow=n,ncol=n*(n-1)/2)
comb <- combn(n,2)
M[cbind(comb[1,],1:(n*(n-1)/2))] <- 1
M[cbind(comb[2,],1:(n*(n-1)/2))] <- -1

M is then a matrix containing all pairwise contrasts for n levels of a factor.

Hope that helps,

Jonathan


On Fri, Oct 15, 2010 at 10:30 AM, Michael Hopkins
<hopkins at upstreamsystems.com> wrote:
>
> On 15 Oct 2010, at 13:55, Berwin A Turlach wrote:
>
>> G'day Michael,
>>
>
> Hi Berwin
>
> Thanks for the reply
>
>> On Fri, 15 Oct 2010 12:09:07 +0100
>> Michael Hopkins <hopkins at upstreamsystems.com> wrote:
>>
>>> OK, my last question didn't get any replies so I am going to try and
>>> ask a different way.
>>>
>>> When I generate contrasts with contr.sum() for a 3 level categorical
>>> variable I get the 2 orthogonal contrasts:
>>>
>>>> contr.sum( c(1,2,3) )
>>>  [,1] [,2]
>>> 1    1    0
>>> 2    0    1
>>> 3   -1   -1
>>
>> These two contrasts are *not* orthogonal.
>>
> I'm surprised.  Can you please tell me how you calculated that.
>
>>> This provides the contrasts <1-3> and <2-3> as expected.  But I also
>>> want it to create <1-2> (i.e. <1-3> - <2-3>).  So in general I want
>>> all possible orthogonal contrasts - think of it as the contrasts for
>>> all pairwise comparisons between the levels.
>>
>> You have to decide what you want.  The contrasts for all pairwise
>> comparaisons between the levels or all possible orthogonal contrasts?
>>
>
> Well the pairwise contrasts are the most important as I am looking for evidence of whether they are zero (i.e. no difference between levels) or not.  But see my above comment about orthogonality.
>
>> The latter is actually not well defined.  For a factor with p levels,
>> there would be p orthogonal contrasts, which are only identifiable up to
>> rotation, hence infinitely many such sets. But there are p(p-1)
>> pairwise comparisons. So unless p=2, yo have to decide what you want....
>>
> Well of course the pairwise comparisons are bi-directional so in fact only p(p-1)/2 are of interest to me.
>
>>> Are there are any options for contrast() or other functions/libraries
>>> that will allow me to do this automatically?
>>
>> Look at package multcomp, in particular functions glht and mcp, these
>> might help.
>>
> Thanks I will have a look.
>
> But I want to be able to do this transparently "within" lm() using regsubsets() etc as I am collecting large quantities of summary stats from all possible models to use with a model choice criterion based upon true Bayesian model probabilities.
>
>> Cheers,
>>
>>       Berwin
>>
>> ========================== Full address ============================
>> Berwin A Turlach                      Tel.: +61 (8) 6488 3338 (secr)
>> School of Maths and Stats (M019)            +61 (8) 6488 3383 (self)
>> The University of Western Australia   FAX : +61 (8) 6488 1028
>> 35 Stirling Highway
>> Crawley WA 6009                e-mail: berwin at maths.uwa.edu.au
>> Australia                        http://www.maths.uwa.edu.au/~berwin
>
>
>
>
> Michael Hopkins
> Algorithm and Statistical Modelling Expert
>
> Upstream
> 23 Old Bond Street
> London
> W1S 4PZ
>
> Mob +44 0782 578 7220
> DL   +44 0207 290 1326
> Fax  +44 0207 290 1321
>
> hopkins at upstreamsystems.com
> www.upstreamsystems.com
>
>
>
>
>        [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org 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.
>