Stéphane Dray me fait remarquer que mon clavier a fourché dans mon
message d'hier. Le test le plus naturel pour des variables explicatives
qualitatives est bien entendu 'randtest.between'. Mais vous aurez tous
corrigé ... hum, hum !
> Si les variables explicatives sont qualitatives 'randtest.discrimin'
> convient très bien. Pour les variables quantitatives, j'ai écrit une
> sorte de 'rtest' pour l'acpvi, qui marche bien mais reste un peu lent.
> Il fait appel au code Fortran de la routine 'dqrls{base}' via la
> fonction 'eig.wfit'. Il est clair qu'un vrai 'randtest.pcaiv' dans
> ade4 serait très utile ... ;-) à Stéphane.
>
> Voici le code en question :
>
> ## Fast function for computing the sum of eigenvalues of the fitted
> table from routines 'dqrls{base}'.
> 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(svd(tab,nu=0,nv=0)$d^2)
> }
>
> rtest.pcaiv<-function(dudiv,nrepet=999)
> { ## from a call to the Fortran routine dqrls{base} via function
> 'eig.wfit'.
> cat("Caution: this test is slow.\n")
> y<-as.matrix(dudiv$Y)
> x<-as.matrix(dudiv$X)
> sqlw<-sqrt(dudiv$lw)
> sqcw<-sqrt(dudiv$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)/dudiv$Dtot
> isim<-c()
> for(i in 1:nrepet)
> isim[i]<-eig.wfit(y,x[sample(nrow(x)),],sqlw,sqcw)/dudiv$Dtot
> return(as.randtest(isim,obs,call=match.call()))
> }
>
> Une petite précaution : il faut rajouter à l'objet 'dudiv' la variable
> 'Dtot' qui correspond à l'inertie totale de l'analyse initiale.
>
> Si l'on reprend l'exemple de '?pcaiv' d'ade4 :
>
> 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)
>
> #Monte-Carlo test
> #Observation: 0.5028486
> #Call: rtest.pcaiv(dudiv = iv1)
> #Based on 999 replicates
> #Simulated p-value: 0.001
>
> Tourne en moins de 10 secondes sur un Mac G4 avec R 1.9 pour OSX.
>
> ***********************************************
> Raphaël Pélissier
> 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/
> ************************************************
***********************************************
Raphaël Pélissier
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/
************************************************
This archive was generated by hypermail 2b30 : Tue Sep 07 2004 - 13:30:57 MEST