May I specify that the rtest.pcaiv() function is definitely not
confidential and available to anybody from my embryonic Diversity.R
package (http://pelissier.free.fr/Diversity.html).
-- *********************************************** Raphael Pelissier UMR AMAP TA40/PS2 34398 Montpellier cedex 5 France Tel. +33 (0)4 67 61 75 23 Fax. +33 (0)4 67 61 56 68 RP Home Page http://pelissier.free.fr/ ************************************************Stephane Dray a écrit :
> Hi Dominique, > to test the significance of the analysis: > - for between, use randtest.between (see ?randtest.between) > - for pcaiv, (for the moment) there is no permutation test > implemented. However, you can use the R function proposed by R. > Pelissier available at: > > http://pbil.univ-lyon1.fr/ADE-4/adelisthtmlannuel/04/0136.html > (login:ade2001, password:passade) > > The message is in french but I think it is quite undestandable. Note > that you must add an element to the initial element (Dtot is sum of > eigenvalues of the separate analysis) before running the function: > > data(rhone) > pca1 <- dudi.pca(rhone$tab, scan = FALSE, nf = 3) > iv1 <- pcaiv(pca1, rhone$disch, scan = FALSE) > iv1$Dtot<-sum(pca1$eig) > rtest.pcaiv(iv1) > > You can also modify the code of the function to remove this step. Here > is a modified version that Raphael sent to me after some comments on > the first version (hope this was not confidential !!): > > eig.wfit<-function(y,x,sqlw,sqcw) { > n<-nrow(x) > p<-ncol(x) > ny<-NCOL(y) > storage.mode(y)<-"double" > > z <- > .Fortran("dqrls",qr=x*sqlw,n=n,p=p,y=y*sqlw,ny=ny,tol=1e-07,coefficients=mat.or.vec(p,ny),residuals=y,effects=mat.or.vec(n,ny),rank > =integer(1),pivot=1:p,qraux=double(p),work=double(2*p),PACKAGE="base") > tab<-sqlw*y-z$residuals > tab<-sweep(tab,2,sqcw,"*") > eig<-sum(tab^2) > } > > > rtest.pcaiv<-function(dudiv,nrepet=999,...) > { ## from a call to the Fortran routine dqrls{base} in > First.lib{diversity} > cat("Caution: this test is slow.\n") > appel<-as.list(dudiv$call) > dudi1<-eval(appel$dudi,sys.frame(0)) > x<-eval(appel$df,sys.frame(0)) > if(is.factor(x)) x<-as.data.frame(x) > inertot<-sum(dudi1$eig) > y<-as.matrix(dudi1$tab) > x<-as.matrix(x) > sqlw<-sqrt(dudi1$lw) > sqcw<-sqrt(dudi1$cw) > fmla<-as.formula(paste("y ~", paste(dimnames(x)[[2]], collapse = > "+"))) > mf<-model.frame(fmla,data=cbind.data.frame(y,x)) > mt<-attr(mf,"terms") > x<-model.matrix(mt,mf) > obs<-eig.wfit(y,x,sqlw,sqcw)/inertot > isim<-c() > for(i in 1:nrepet) > isim[i]<-eig.wfit(y,x[sample(nrow(x)),],sqlw,sqcw)/inertot > return(as.randtest(isim,obs,call=match.call())) > } > > > Copy the two functions above in a text file, save it and then source > the file in R, and then run: > > data(rhone) > pca1 <- dudi.pca(rhone$tab, scan = FALSE, nf = 3) > iv1 <- pcaiv(pca1, rhone$disch, scan = FALSE) > rtest.pcaiv(iv1) > > Sincerely. > > > > > At 03:10 19/01/2005, Dominique Grüter wrote: > >> Dear list >> I used the "between" - and the "pcaiv"-command. >> How can I find out if the results of these two commands are significant? >> (The distribution of my normal "dudi. pca" - results is not normal ) >> >> Thanks >> >> domi > > > Stéphane DRAY > -------------------------------------------------------------------------------------------------- > > Département des Sciences Biologiques > Université de Montréal, C.P. 6128, succursale centre-ville > Montréal, Québec H3C 3J7, Canada > > Tel : (514) 343-6111 poste 1233 Fax : (514) 343-2293 > E-mail : stephane.dray@umontreal.ca > -------------------------------------------------------------------------------------------------- > > Web > http://www.steph280.freesurf.fr/ > -------------------------------------------------------------------------------------------------- > > >
This archive was generated by hypermail 2b30 : Tue Oct 18 2005 - 10:35:39 MEST