[R] help with aggregate(): tables of means for terms in an mlm

Michael Friendly friendly at yorku.ca
Tue Aug 28 22:04:49 CEST 2007


I'm trying to extend some work in the car and heplots packages
that requires getting a table of multivariate means for one
(or later, more) terms in an mlm object.  I can do this for
concrete examples, using aggregate(), but can't figure out how to 
generalize it.  I want to return a result that has the factor-level
combinations as rownames, and the means as the body of the table
(aggregate returns the factors as initial columns).

# Examples: m1 & m2 are desired results

 > library(car)
 > soils.mod <- lm(cbind(pH,N,Dens,P,Ca,Mg,K,Na,Conduc) ~ Block + 
Contour*Depth, data=Soils)
 > term.names(soils.mod)
[1] "(Intercept)"   "Block"         "Contour"       "Depth"
[5] "Contour:Depth"
 >

 > # response variables
 > resp<- model.response(model.frame(soils.mod))
 > # 1-factor means: term="Contour"
 > m1<-aggregate(resp, list(Soils$Contour), mean)
 > rownames(m1) <- m1[,1]
 > ( m1 <- m1[,-1] )
               pH       N  Dens     P    Ca    Mg      K    Na Conduc
Depression 4.692 0.08731 1.343 188.2 7.101 8.986 0.3781 5.823  6.946
Slope      4.746 0.10594 1.333 159.4 8.109 8.320 0.4156 6.103  6.964
Top        4.570 0.11256 1.272 150.9 8.877 8.088 0.6050 4.872  5.856
 >

 > # 2-factor means:  term="Contour:Depth"
 > m2<-aggregate(resp, list(Soils$Contour, Soils$Depth), mean)
 > rownames(m2) <- paste(m2[,1], m2[,2],sep=":")
 > ( m2 <- m2[,-(1:2)] )
                     pH       N   Dens     P     Ca    Mg      K     Na 
Conduc
Depression:0-10  5.353 0.17825 0.9775 333.0 10.685 7.235 0.6250 1.5125 
1.473
Slope:0-10       5.508 0.21900 1.0500 258.0 12.248 7.232 0.6350 1.9900 
2.050
Top:0-10         5.332 0.19550 1.0025 242.8 13.385 6.590 0.8000 0.9225 
1.373
Depression:10-30 4.880 0.08025 1.3575 187.5  7.548 9.635 0.4500 4.6400 
5.480
Slope:10-30      5.283 0.10100 1.3475 160.2  9.515 8.980 0.4800 4.9350 
4.910
Top:10-30        4.850 0.11750 1.3325 147.5 10.238 8.090 0.6500 2.9800 
3.583
Depression:30-60 4.362 0.05050 1.5350 124.2  5.402 9.918 0.2400 7.5875 
9.393
Slope:30-60      4.268 0.06075 1.5100 114.5  5.877 8.968 0.3000 7.6300 
8.925
Top:30-60        4.205 0.07950 1.3225 116.2  6.620 8.742 0.5450 6.2975 
7.440
Depression:60-90 4.173 0.04025 1.5025 108.0  4.770 9.157 0.1975 9.5525 
11.438
Slope:60-90      3.927 0.04300 1.4225 105.0  4.798 8.100 0.2475 9.8575 
11.970
Top:60-90        3.893 0.05775 1.4300  97.0  5.268 8.928 0.4250 9.2900 
11.030
 >

Here is the current version of a function that doesn't work, because I
can't supply the factor names to aggregate in the proper way.
Can someone help me make it work?

termMeans.mlm <- function( object, term ) {
    resp<- model.response(model.frame(object))
    terms <- term.names(soils.mod)
    terms <- terms[terms != "(Intercept)"]
    factors <- strsplit(term, ":")
    # browser()
    means <- aggregate(resp, factors, mean)
    # rownames(means) <- ...
    # means <- means[, -(1:length(factors)]
}

 > termMeans.mlm(soils.mod, "Contour")
Error in FUN(X[[1L]], ...) : arguments must have same length

thanks,
-Michael

-- 
Michael Friendly     Email: friendly AT yorku DOT ca
Professor, Psychology Dept.
York University      Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Street    http://www.math.yorku.ca/SCS/friendly.html
Toronto, ONT  M3J 1P3 CANADA



More information about the R-help mailing list