Le vendredi, 13 aoû 2004, à 07:35 US/Pacific, Stephane DRAY a écrit :
> Bonjour,
>
> At 04:51 13/08/2004, Charline Laurent wrote:
>> Bonjour,
>>
>> Je dois traiter un jeu de données d'acoustique (écho-intégration par
>> bancs). Je dispose de 10 variables quantitatives sur 6236 relevés.
>> Je dispose également de 2 variables quantitatives indiquant pour
>> chaque relevé le pays (2 modalités) et si le relevé a été effectué de
>> jour ou de nuit (variable "diel").
>> Je dois voir s'il existe un effet "pays", un effet "diel" (voire une
>> interaction).
>>
>> J'ai réalisé une analyse inter-classe (pays) : le pourcentage
>> d'inertie expliquée par l'analyse est de 0.8%. L'analyse inter-classe
>> (diel) me donne un pourcentage d'inertie expliquée de 2.3 %.
>> Les test de permutation (fonction randtest.between) me donnent pour
>> ces deux analyses inter-classes me donnent des résultats
>> significatifs.
>>
>> Des pourcentages d'inertie expliquée très faible mais des tests de
>> permutation très significatifs... cela est-il très cohérent ?
>
> Vu le grand nombre d'individus, cela n'est pas tres etonnant. Il y a
> beaucoup de variance intra-classe car on a beaucoup d'individus
>
>> Une autre petite question : Je voudrais réaliser une ACPVI sur ces
>> données, mais je n'arrive pas à trouver dans les sorties de R les
>> pourcentages d'inertie expliquée par chaque facteur, et comment
>> tester la significativité de chaque facteur.
>
> Ca revient au meme !! L'analyse inter-classe est une ACPVI. Pour s'en
> convaincre:
>
> data(meaudret)
> pca1 <- dudi.pca(meaudret$mil, scan = FALSE, nf = 4)
> bet1 <- between(pca1, meaudret$plan$sta, scan = FALSE, nf = 2)
> acpvi1=pcaiv(pca1, meaudret$plan$sta,scan = FALSE, nf = 2)
>
> > acpvi1$eig
> [1] 2.681164426 0.620767949 0.113168924 0.009502576
> > bet1$eig
> [1] 2.681164426 0.620767949 0.113168924 0.009502576
>
> Ces valeurs propres sont des variances expliquees.
> Le pourcentage de variance expliquee est simplement la rapport inertie
> de l'analyse simple/inertie analyse contrainte:
> sum(acpvi1$eig)/sum(pca1$eig)
> [1] 0.3805115
>
> C'est le pourcentage de variance interclasse.
>
> on retrouve ce rapport dans la sortie de between:
> > bet1$ratio
> [1] 0.3805115
>
> Le pourcentage de variance expliquee par chaque axe:
>
> > (acpvi1$eig)/sum(pca1$eig)
> [1] 0.297907158 0.068974217 0.012574325 0.001055842
>
> Pour le tester, c'est une autre affaire... c'est pas possible avec
> ade. Canoco propose un test pour le premier axe ou pour tous les axes
> mais pas de tester chaque axe. On peut envisager de tester par
> permutation axe1, axe1+axe2, axe1+axe2+axe3... mais est-ce bien utile
> ??
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/
************************************************
>> Encore merci
>> Charline L.
>
> De rien
>
> Cordialement.
>
> 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 Sep 07 2004 - 13:30:57 MEST