Hi All! I wish to use R for MANOVA analysis. Where could I get info on that? __________________________________________________ Do You Yahoo!? Send instant messages & get email alerts with Yahoo! Messenger. http://im.yahoo.com/ -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On Mon, 8 May 2000, mandar oak wrote:> I wish to use R for MANOVA analysis. Where could I get > info on that?Depends exactly what you mean by MANOVA. Try ?aov which handles multiple responses. However, the MANOVA multivariate tests (Pillai, Wilks lambda ...) are not implemented in base R. We would welcome submission of code for those, of course. -- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Prof Brian D Ripley <ripley at stats.ox.ac.uk> writes:> On Mon, 8 May 2000, mandar oak wrote: > > > I wish to use R for MANOVA analysis. Where could I get > > info on that? > > Depends exactly what you mean by MANOVA. Try ?aov which handles multiple > responses. However, the MANOVA multivariate tests (Pillai, Wilks lambda > ...) are not implemented in base R. We would welcome submission of code > for those, of course.At lest some of that sits in contributed packages "multilm" and "norm". (Both of those are on my list of "things to try when I get the time"... I believe multilm more or less *is* MANOVA, whereas norm handles estimation in the presence of missing values. One interesting question is whether you can do MANOVA with missings by combining the two.) -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /''_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On 9 May 2000, Peter Dalgaard BSA wrote:> Prof Brian D Ripley <ripley at stats.ox.ac.uk> writes: > > > On Mon, 8 May 2000, mandar oak wrote: > > > > > I wish to use R for MANOVA analysis. Where could I get > > > info on that? > > > > Depends exactly what you mean by MANOVA. Try ?aov which handles multiple > > responses. However, the MANOVA multivariate tests (Pillai, Wilks lambda > > ...) are not implemented in base R. We would welcome submission of code > > for those, of course. > > At lest some of that sits in contributed packages "multilm" and > "norm". (Both of those are on my list of "things to try when I get the > time"... I believe multilm more or less *is* MANOVA, whereas norm > handles estimation in the presence of missing values. One interesting > question is whether you can do MANOVA with missings by combining the > two.)Not really any of it! multilm fits multivariate linear models and does Hotelling T^2 tests. It does not group terms for AOV, and I think is only applicable for continuous variates (not factor explantory variables). Certainly the help and examples only discuss that case. norm is about multivariate normal data, no regressors, AFAIK. Various people have pointed out the dangers of that: how often is iid on a multivriate normal population appropriate? -- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Prof Brian D Ripley <ripley at stats.ox.ac.uk> writes:> Not really any of it! multilm fits multivariate linear models and does > Hotelling T^2 tests. It does not group terms for AOV, and I think is only > applicable for continuous variates (not factor explantory variables). > Certainly the help and examples only discuss that case.Had a closer look. The factor expl.var. would seem to be there case-wise, so you can build the models, but Wilk''s Lambda and friends are absent, so you cannot (easily) do multi-df model reduction tests. However, it''s no too far off. Is the author planning to develop it further?> norm is about multivariate normal data, no regressors, AFAIK. Various > people have pointed out the dangers of that: how often is iid on a > multivriate normal population appropriate?Not often, but sometimes things fall apart into several homogeneous groups and then it can be used as a building block. I''ve used something like this for an RCT a while back back. -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /''_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On 9 May 2000, Peter Dalgaard BSA wrote:> Prof Brian D Ripley <ripley at stats.ox.ac.uk> writes: > > > Not really any of it! multilm fits multivariate linear models and does > > Hotelling T^2 tests. It does not group terms for AOV, and I think is only > > applicable for continuous variates (not factor explantory variables). > > Certainly the help and examples only discuss that case. > > Had a closer look. The factor expl.var. would seem to be there > case-wise, so you can build the models, but Wilk''s Lambda and friends > are absent, so you cannot (easily) do multi-df model reduction tests. > However, it''s no too far off. Is the author planning to develop it > further?Now it is my turn, I think :-) 1) "multilm" is on my list, but there are some theoretical problems (e.g. is the use of Moore-Penrose numerically appropriate ?) Any help by other distributors would be very welcome! 2) "multilm" was written for the use of the stabilized multivariate tests by Laeuter and Kropf, T^2 is only a by-product. Therefore the documentation is heavily biased on this. 3) I do think that one can take factors as design variables (there is an Iris example, I think). 4) I just had my third move this year (this time to Erlangen, Bavaria), so time is a problem too :-) Torsten -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On Tue, 9 May 2000, Torsten Hothorn wrote:> On 9 May 2000, Peter Dalgaard BSA wrote: > > > Prof Brian D Ripley <ripley at stats.ox.ac.uk> writes: > > > > > Not really any of it! multilm fits multivariate linear models and does > > > Hotelling T^2 tests. It does not group terms for AOV, and I think is only > > > applicable for continuous variates (not factor explantory variables). > > > Certainly the help and examples only discuss that case. > > > > Had a closer look. The factor expl.var. would seem to be there > > case-wise, so you can build the models, but Wilk''s Lambda and friends > > are absent, so you cannot (easily) do multi-df model reduction tests. > > However, it''s no too far off. Is the author planning to develop it > > further? > > Now it is my turn, I think :-) > > 1) "multilm" is on my list, but there are some theoretical problems > (e.g. is the use of Moore-Penrose numerically appropriate ?) Any help by > other distributors would be very welcome! > > 2) "multilm" was written for the use of the stabilized multivariate tests > by Laeuter and Kropf, T^2 is only a by-product. Therefore the > documentation is heavily biased on this. > > 3) I do think that one can take factors as design variables (there is an > Iris example, I think).One can fit the model, but as the Iris example shows, it assumes three continuous vars, NOT one factor in a model formula, and so the information given is not appropriate for a manova as I understand that term. Manova includes a series of tests of submodels. Here''s an S example of manova: wafer.manova <- manova(cbind(pre.mean, post.mean) ~ maskdim + visc.tem + spinsp, wafer) summary(wafer.manova) # manova table with Pillai''s trace Df Pillai Trace approx. F num df den df P-value maskdim 1 0.42662 4.09229 2 11 0.04693 visc.tem 2 0.5809 2.45605 4 24 0.07305 spinsp 2 0.47904 1.88978 4 24 0.14495 Residuals 12 so factors in the model are grouped and several submodels have been fitted. -- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
I hope this could help. It was originally written in Spanish, so several errors migth have appeared in the translation process. If so, let me know. # Let?s define some kind of determinant function> det <- function (x) { prod(eigen(x)$values) }# It coul be useful a function to make void lists with as many components as the parameter ncomp states.> DoClassList <- function(ncomp)+{ + aux list(''name''=NULL) + if (ncomp>1) + for (i in 2:ncomp) + aux <- c(aux, list(''name''=NULL)) + for (i in 1:ncomp) + attributes(aux)$names[i] <- paste(''Class'',i) + return(aux) +} # Let?s define a function to evaluate Wilks statistic.> TestWilk function (data.in)+ { + nclass <-length(data.in) + nvars <- ncol(data.in[[1]]) + aveQ <- DoClassList(nclass) + for (i in 1:nclass) + aveQ[[i]] <- apply(data.in[[i]],2,mean) + SQ <- lapply(data.in,var) + W <- SQ[[1]]*(nrow(data.in[[1]])-1) + for (i in 2:nclass) + W <- W + SQ[[i]]*(nrow(data.in[[i]])-1) + aveTotal <- aveQ[[1]]*nrow(data.in[[1]]) + for (i in 2:nclass) + aveTotal <- aveTotal + aveQ[[i]] * + nrow(data.in[[i]]) + aveTotal <- aveTotal / sum(unlist( lapply(data.in,nrow))) + B <- nrow(data.in[[1]]) * as.matrix(aveQ[[1]] - aveTotal) + %*% t(as.matrix(aveQ[[1]] - aveTotal)) + for (i in 2:nclass) + B <- B + nrow(data.in[[i]]) * as.matrix(aveQ[[i]] - + aveTotal) %*% t(as.matrix(aveQ[[i]] - aveTotal)) + n <- nvars + s <- sum(unlist(lapply(data.in,nrow))) - nclass + t <- nclass -1 + Lambda <- det(W)/det(B+W) + w2 <- 1- (nrow(data.in[[i]])*Lambda)/(s+Lambda) + w <- w2 - (n^2+t^2)/(nrow(data.in[[i]])*3)*(1-w2) + if ((n^2+t^2-5) != 0) + { + a <- s + t - (n + t + 1)/2 + b <- sqrt((n^2*t^2-4)/(n^2+t^2-5)) + c <- (n*t-2)/4 + doubt <- (a*b - 2*c) /(n*t) * (Lambda^(-1/b)-1) + return (list( valor=doubt, F=qf(0.99, n*t, (a*b - 2*c)), + L=(doubt < qf(0.99, n*t, (a*b - 2*c))), w=w)) + } + else + { + doubt <- -( s - 0.5*(n-t+1))*log(Lambda) + return (list( value=doubt, chi2=qchisq(0.99, n*t), + L=(doubt < qchisq(0.99, n*t)),w=w)) + } +} As an example of use let?s consider two binormal distributions class1 and class2, with different mean vector but identical covariance matrix.> class1_v1 <- rnorm(5000, mean=0, sd=1) > class1_v2 <- rnorm(5000, mean=0, sd=1) > class1 <- cbind(clase1_v1,clase1_v2) > class2_v1 <- rnorm(5000, mean= 1.5, sd=1) > class2_v2 <- rnorm(5000, mean=2, sd=1) > class2 <- cbind(clase2_v1, clase2_v2)# Lets represent the two classes> plot(rbind(class1,class2),+ col=c(rep("Blue",length(clase1[,1])), + rep("Red",length(clase2[,1]))) , pch=".")> vcellipse(apply(clase1,2,mean), cov(clase1), col="Blue") > vcellipse(apply(clase2,2,mean), cov(clase2), col="Red")# I hope you too have the vcellipse code that was sent to R-help. Nice function ideed. # I would have sent the picture but I?m sure that to be allowed in the list. ?? # And apply the test> TwoClasses <- list(class1, class2) > TestTwoClasses <- TestWilk(TwoClasses) > unlist(TestTwoClasses)value chi2 L w 9403.4041430 9.2103404 FALSE 0.8047111 As for the example distributions normality condition can be applied, this result says NO to the null hipothesis and two different classes of behaviour have been identified. Of course, in a more general case, you must check that normality and homogeneity of covariance matrix can be assesed. But that is another war.... :-) ---------------------------------------------------------------------------- Manuel Castej?n Limas. Project Management Area. e-mail:manuel.castejon at dim.unirioja.es Mechanical Engineering Dept. http://www.unirioja.es/dptos/dim University of La Rioja c/ Luis de Ulloa 20, 26004 Logro?o La Rioja Spain. ---------------------------------------------------------------------------- -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>>>>> Manuel Castejon Limas <manuel.castejon at dim.unirioja.es> writes:... Manuel> # Let´s define some kind of determinant function >> det <- function (x) { prod(eigen(x)$values) } .... We have been there before (a year ago on R-devel, see "eigenvalue"), and since then had the following in ?qr :>> Examples: >> >> ## The determinant of a matrix -- if you really must have it >> det2 <- function(x) prod(diag(qr(x)$qr))*(-1)^(ncol(x)-1)But since people seem to re-invent that teeny part of the wheel, something like the following will be in R 1.1.x : ### From Doug Bates'' 20 Apr 1999 post to R-devel; ### "method" idea from J.K.Lindsey''s rmutil det <- function(x, method = c("qr","eigenvalues")) { if(!is.matrix(x)) stop("x must be a matrix") method <- match.arg(method) # ensures a method from above if(method == "qr") prod(diag(qr(x)$qr)) *(-1)^(ncol(x)-1) else ## if(method == "eigenvalues") Re(prod(eigen(x, only.values=TRUE)$values)) } det package:base R Documentation Calculate the Determinant of a Matrix Description: `det'' calculates the determinant of a matrix either by `qr'' decomposition or from the `eigen''values. Usage: det(x, method = c("qr","eigenvalues")) Arguments: z: numeric matrix. method: `"qr"'' (default) or `"eigenvalues"'' Note: Often, computing the determinant is not what you should be doing to solve a given problem. The `"qr"'' method is much faster for large matrices. See Also: `eigen'', `qr'', `svd'' Examples: (x <- matrix(1:4, ncol=2)) det(x) det(x, method="eigenvalues") det(print(cbind(1,1:3,c(2,0,1)))) -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Dear people, Besides of being looking for outlier detection methods in non normal multivariate distributions I am also interested in knowing about the different ways people manage the bad behaviour of MANOVA as soon as the algorithm is used with irregular number of samples in each class. Any help would be appreciated. And excuse me if I tried to write a det function when I needed to calcute one and found none. I promise not to do it again :-) I hope you understand that a year ago I wasn?t in the list, so I missed that part. Anyway it was not mandatory to use my det function. If it was able to make you all know that the det function in TestWilk meant determinant, I?m happy enough. -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html Send "info", "help", or "[un]subscribe" (in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._