# Re: between and pcaiv - command

From: Stephane Dray (dray@biomserv.univ-lyon1.fr)
Date: Wed Jan 19 2005 - 17:02:49 MET

• Next message: Raphael Pelissier: "Re: between and pcaiv - command"

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:

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