Francisco Mora Ardila
2011-Nov-05 21:01 UTC
[R] testing significance of axis loadings from multivariate dudi.mix
Hi all I?m trying to tests the significance of loadings from a ordination of 46 variables (caategorical, ordinal and nominal). I used dudi.mix from ade4 for the ordination. A years ago Jari Oksanen wrote this script implementing Peres-Neto et al. 2003 (Ecology) bootstraping method: netoboot <- function (x, permutations=1000, ...) { pcnull <- princomp(x, cor = TRUE, ...) res <- pcnull$loadings out <- matrix(0, nrow=nrow(res), ncol=ncol(res)) N <- nrow(x) for (i in 1:permutations) { pc <- princomp(x[sample(N, replace=TRUE), ], cor = TRUE ...) pred <- predict(pc, newdata = x) r <- cor(pcnull$scores, pred) k <- apply(abs(r), 2, which.max) reve <- sign(diag(r[k,])) sol <- pc$loadings[ ,k] sol <- sweep(sol, 2, reve, "*") out <- out + ifelse(res > 0, sol <= 0, sol >= 0) } out/permutations } I tried to aply it to the case of dudi.mix instead of princomp this way: netoboot1<-function (x, permutations=1000,...) { dudinull <- dudi.mix(x, scannf = FALSE, nf = 3) res <- dudinull$c1 out <- matrix(0, nrow=nrow(res), ncol=ncol(res)) N <- nrow(x) for (i in 1:permutations) { dudi <- dudi.mix(x[sample(N, replace=TRUE), ], scannf = FALSE, nf = 3) pred <- predict(dudi, newdata = x) r <- cor(dudinull$li, pred) k <- apply(abs(r), 2, which.max) reve <- sign(diag(r[k,])) sol <- dudi$c1[ ,k] sol <- sweep(sol, 2, reve, "*") out <- out + ifelse(res > 0, sol <= 0, sol >= 0) } out/permutations } But a problem arised with the predict function: it doesn?t seem to work with an object from dudi.mix and I dont understand why. Can somebody tell me why? Any suggestions to modify the script or to use other method? Thanks in advance. Francisco Francisco Mora Ardila Laboratorio de Biodiversidad y Funcionamiento del Ecosistema Centro de Investigaciones en Ecosistemas UNAM-Campus Morelia Tel 3222777 ext. 42621 Morelia , MIchoac?n, M?xico. -- Open WebMail Project (http://openwebmail.org)
Mark Difford
2011-Nov-06 09:33 UTC
[R] testing significance of axis loadings from multivariate dudi.mix
On Nov 05, 2011 at 11:01pm Francisco Mora Ardila wrote:> But a problem arised with the predict function: it doesn?t seem to work > with an object > from dudi.mix and I dont understand why.Francisco, There is no predict() method for dudi.mix() or for any of the dudi objects in ade4. I don't see why you can't get around this by doing something like the following, but you need to take account of any scaling/centring that you might do to your data before calling dudi.mix(). ## Does a dudi.mix on continuous data, so really equals a dudi.pca/princomp/PCA library(ade4) data(deug) deug.dudi <- dudi.mix(deug$tab, scann=F, nf=2) tt <- as.matrix(deug.dudi$tab) %*% as.matrix(deug.dudi$c1) ## see note below qqplot(deug.dudi$li[,1], tt[,1]) qqplot(deug.dudi$li[,2], tt[,2]) deug.princ <- princomp(deug$tab, cor=T) qqplot(predict(deug.princ)[,1], tt[,1]) ## scaling not accounted for: deug.princ <- princomp(deug$tab, cor=F) qqplot(predict(deug.princ)[,1], tt[,1]) rm(tt, deug.dudi, deug.princ) Note that in the code given above, "as.matrix(deug.dudi$tab) %*% as.matrix(deug.dudi$c1)" is based on how stats:::predict.princomp does it. Regards, Mark. ----- Mark Difford (Ph.D.) Research Associate Botany Department Nelson Mandela Metropolitan University Port Elizabeth, South Africa -- View this message in context: http://r.789695.n4.nabble.com/testing-significance-of-axis-loadings-from-multivariate-dudi-mix-tp3994281p3995350.html Sent from the R help mailing list archive at Nabble.com.