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 : Wed Jan 19 2005 - 17:04:22 MET