Re: summary(object, test=c("Roy", "Wilks", "Pillai", ....)

[Ray Haraf <rayharaf@rogers.com>]

Thanks, Professor Venables, for your prompt and valuable reply. I tried manova() function and got 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

I am still struggling 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

 

Lastly but not the least, would be too ambitious to draw the axes (i.e, the eigenvalues) to  the ellipse? fear God, honour the king and work very hard!

 

Thanks and very kind regards,

Ray

 

Bill.Venables@csiro.au wrote:

With a matris response lm(...) doesn't really fit a multivariate
regression. It fits multiple univariate linear regressions. You can
build a multivariate linear regression from it, of course.

I suspect you need to looking at manova(...) to fit the model. Then the
summary method, and the other model comparison tools of course, will
give you multivariate tests.

As to simultaneous prediction intervals and ellipsoids, the only advice
I can give you is the standard: fear God, honour the king and work very
hard. Best of luck! :-)

Bill Venables.


-----Original Message-----
From: s-news-owner@lists.biostat.wustl.edu
[mailto:s-news-owner@lists.biostat.wustl.edu] On Behalf Of Ray Haraf
Sent: Sunday, 30 March 2008 3:32 PM
To: s-news@lists.biostat.wustl.edu
Subject: [S] summary(object, test=c("Roy", "Wilks", "Pillai", ....)

Dear All,

I am running multivariate multiple regression using lm() and would be
very appreciative of your kind and prompt help (as usual) with how to
obtain the usual statistics of multivariate analysis, i.e., Wilks'
lambda, Pillai's trace, Hotelling-Lawley trace, and Roy's greatest root.

I unsuccessfully tried this
> ex7.8 <- data.frame(z1 = c(0, 1, 2, 3, 4),
+ y1 = c(1, 4, 3, 8, 9),
+ y2 = c(-1, -1, 2, 3, 2))
> summary(lm(cbind(y1, y2) ~ z1, data ="">ex7.8),test=c("Wilks","Roy","Pillai"))

Also, any suggestion on how obtain simultaneous prediction intervals and
prediction ellipsoid would be very valued. Thanks for your kind care

Here is the result of the above summary() function

Response y1 :
Call:
lm(formula = y1 ~ z1, data = "">Residuals:
1 2 3 4 5
6.880e-17 1.000e+00 -2.000e+00 1.000e+00 -1.326e-16
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.0000 1.0954 0.913 0.4286
z1 2.0000 0.4472 4.472 0.0208 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.414 on 3 degrees of freedom
Multiple R-squared: 0.8696, Adjusted R-squared: 0.8261
F-statistic: 20 on 1 and 3 DF, p-value: 0.02084

Response y2 :
Call:
lm(formula = y2 ~ z1, data = "">Residuals:
1 2 3 4 5
9.931e-17 -1.000e+00 1.000e+00 1.000e+00 -1.000e+00
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.0000 0.8944 -1.118 0.3450
z1 1.0000 0.3651 2.739 0.0714 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.155 on 3 degrees of freedom
Multiple R-squared: 0.7143, Adjusted R-squared: 0.619
F-statistic: 7.5 on 1 and 3 DF, p-value: 0.07142

>
Kind regards,
Ray