[R] mixtures as outcome variables
James Reilly
reilly at stat.auckland.ac.nz
Fri Mar 25 17:43:10 CET 2005
A collection of functions for compositional data analysis were posted on
the S-news mailing list about a year ago.
Basic Compositional Data Analysis functions for S+/R
http://www.biostat.wustl.edu/archives/html/s-news/2003-12/msg00139.html
James
> Date: Thu, 24 Mar 2005 20:28:51 -0400
> From: Kjetil Brinchmann Halvorsen <kjetil at acelerate.com>
> Subject: Re: [R] mixtures as outcome variables
> Cc: r-help at stat.math.ethz.ch, "Jason W. Martinez" <jmartinez5 at verizon.net>
>
> Kjetil Brinchmann Halvorsen wrote:
>
>>> Jason W. Martinez wrote:
>>>
>>
>>>>> Dear R-users,
>>>>>
>>>>> I have an outcome variable and I'm unsure about how to treat it. Any
>>>>> advice?
>>>>>
>>
>>> If you only concentrate on the relative proportions, this are called
>>> compositional data. I f your data are in
>>> mydata (n x 4), you obtain compositions by
>>> sweep(mydata, 1, apply(mydata, 1, sum), "/")
>>>
>>> There are not (AFAIK) specific functions/packages for R for
>>> compositional data AFAIK, but you
>>> can try googling. Aitchison has a monography (Chapman & Hall) and a
>>> paper in JRSS B.
>>>
>>> One way to start might be lm's or anova on the symmetric logratio
>>> transform of the
>>> compositons. The R function lm can take a multivariate response, but
>>> some extra programming will be needed
>>> for interpretation. With simulated data:
>>>
>>
>>>> > slr
>>
>>> function(y) { # y should sum to 1
>>> v <- log(y)
>>> return( v - mean(v) ) }
>>
>>>> > testdata <- matrix( rgamma(120, 2,3), 30, 4)
>>>> > str(testdata)
>>
>>> num [1:30, 1:4] 0.200 0.414 0.311 2.145 0.233 ...
>>
>>>> > comp <- sweep(testdata, 1, apply(testdata,1,sum), "/")
>>
>>> # To get the symmetric logratio transform:
>>> comp <- t(apply(comp, 1, slr))
>>> # Observe:
>>> apply(cov(comp), 1, sum)
>>> [1] -5.551115e-17 2.775558e-17 5.551115e-17 -2.775558e-17
>>
>>>> > lm( comp ~ 1)
>>
>>>
>>> Call:
>>> lm(formula = comp ~ 1)
>>>
>>> Coefficients:
>>> [,1] [,2] [,3] [,4] (Intercept)
>>> 0.17606 0.06165 -0.03783 -0.19988
>
>
> Followup:
>
> > mmod <- manova(comp ~ x)
> > summary(mmod)
> Error in summary.manova(mmod) : residuals have rank 3 < 4
> >
>
> So the manova() function cannot be used. I guess MANOVA for
> compositional data should be
> a straight extension, but it must be programmed , standard manova cannot
> be used.
>
> Kjetil
>
> -- Kjetil Halvorsen. Peace is the most effective weapon of mass construction. -- Mahdi Elmandjra
More information about the R-help
mailing list