[R] summary(object, test=c("Roy", "Wilks", "Pillai", ....) AND ellipse(object, center=....)

Sun Apr 6 10:37:41 CEST 2008

Yes, it does mean you can only get one kind of test at a time.

If you want more than one kind of test, you have to do it separately.
The simplest way to do this is probably something like this:

M_obj <- manova(cbind(y1, y2) ~ z1, data = ex7.8)

for(test in c("Wilks", "Roy", "Pillai"))
	print(summary(M_obj, test = test)) 

Further, I do not anticipate R will change to make this sort of thing
any easier.  
I know some of the big guys offer multiple kinds of test all at once, no
questions asked, presumably to allow the user to pick and choose post
hoc.  Perhaps we can persuade them to adopt better practices in future,
but somehow I doubt it.

Bill Venables
CSIRO Laboratories
PO Box 120, Cleveland, 4163
AUSTRALIA
Office Phone (email preferred): +61 7 3826 7251
Fax (if absolutely necessary):  +61 7 3826 7304
Mobile:                         +61 4 8819 4402
Home Phone:                     +61 7 3286 7700
mailto:Bill.Venables at csiro.au
http://www.cmis.csiro.au/bill.venables/ 

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On Behalf Of Ray Haraf
Sent: Sunday, 6 April 2008 3:16 AM
To: r-help at r-project.org
Subject: [R] summary(object, test=c("Roy", "Wilks", "Pillai",....) AND
ellipse(object, center=....)

  Dear All,

  I would be very appreciative of your help with the following

  1). I am running multivariate multiple regression through the manova()
function (kindly suggested by Professor Venables) and getting two
different answers for test=c("Wilks","Roy","Pillai") and
tests=c("Wilks","Roy",'"Pillai") as shown below. In the first case
(test=c(list)) I got error message which probably means I can only call
one test at a time. I thought I could get ride of this by adding "s" to
test; in this case (tests=c(list)), I got Pillai test. Does this mean
that Pillai would be the default test and summary(manova()) can only
post one test at a time?

> summary(manova(cbind(y1, y2) ~ z1, data =
+ ex7.8),test=c("Wilks","Roy","Pillai"))
Error in match.arg(test) : 'arg' must be of length 1
> summary(manova(cbind(y1, y2) ~ z1, data =
+ ex7.8),tests=c("Wilks","Roy","Pillai"))
          Df  Pillai approx F num Df den Df Pr(>F)  
z1         1  0.9375  15.0000      2      2 0.0625 .
Residuals  3                                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

2). My next struggle is to construct prediction ellipse. Both ellipse()
and ellipse.lm() are not giving me the solution to "Sampling from
multivariate multiple regression prediction regions" posted by Iain
Pardoe, Mon May 9 18:43:46 2005. I am working on the same problem and
performed all the steps he suggested

> ex7.10 <- 
+ data.frame(y1 = c(141.5, 168.9, 154.8, 146.5, 172.8, 160.1, 108.5), 
+ y2 = c(301.8, 396.1, 328.2, 307.4, 362.4, 369.5, 229.1), 
+ z1 = c(123.5, 146.1, 133.9, 128.5, 151.5, 136.2, 92), 
+ z2 = c(2.108, 9.213, 1.905, .815, 1.061, 8.603, 1.125)) 
> attach(ex7.10) 
> f.mlm <- lm(cbind(y1,y2)~z1+z2) 
> y.hat <- c(1, 130, 7.5) %*% coef(f.mlm) 
> round(y.hat, 2) 
         y1 y2 
[1,] 151.84 349.63 
> qf.z <- t(c(1, 130, 7.5)) %*% 
+ solve(t(cbind(1,z1,z2)) %*% cbind(1,z1,z2)) %*% 
+ c(1, 130, 7.5) 
> round(qf.z, 5) 
        [,1] 
[1,] 0.36995 
> n.sigma.hat <- SSD(f.mlm)$SSD # same as t(resid(f.mlm)) %*% 
resid(f.mlm) 
> round(n.sigma.hat, 2) 
     y1 y2 
y1 5.80 5.22 
y2 5.22 12.57 
> F.quant <- qf(.95,2,3) 
> round(F.quant, 2) 
[1] 9.55 

>From here how could I calculate a 95% prediction ellipse for y=(y1,y2)
at (z1,z2)=(130,7.5) using either ellipse or ellipse.lm? y1 would be the
x-axis and y2, the y-axis. The center is different from (0,0) and I
don't know what would be the appropriate x (the lm object). Should I
used predicted values or residuals? In both cases I have vectors which
is different from the example given with ellipse.lm

3). Lastly but not the least, would be too ambitious to draw the axes
(i.e, the eigenvalues) to  the ellipse?

Thanks and very kind regards,
Ray

	[[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.