Thomas, 1) Presence absence data means that you have cell probabilities 1/S_i for detections and 0 for missing species in a given community i. As Zoltán also pointed out, it is meaningful to use this, as it has the interpretation of choosing different species from a species list (and not from a pool of individuals).
2) Calculations never end. You might have a bad luck with your data. The underlying algorithm (see refs on ade4:::divcmax help page, I am referring to the underlying functions and not the wrapper from FD) is not very robust and it often fails to find a solution for the max Rao's quadratic entropy value. It is legitimate to use this with 0/1 data, and it hangs up even with count based probabilities quite often according to my experience. I haven't spend much time figuring out a fix, but there are 2 loops in the code without a break. Below is a function where you can specify maxit argument so that it quits before the end of times. I am cc-in to Simon Penel, maintainer of ade4, as this seems to me as a possible feature they might want to add. Cheers, Peter -- Péter Sólymos, Dept Biol Sci, Univ Alberta, T6G 2E9, Canada AB soly...@ualberta.ca, Ph 780.492.8534, http://psolymos.github.com Alberta Biodiversity Monitoring Institute, http://www.abmi.ca Boreal Avian Modelling Project, http://www.borealbirds.ca ## -- code starts here divcmax_2 <- function (dis, epsilon = 1e-08, comment = FALSE, maxit = 1000) { if (!inherits(dis, "dist")) stop("Distance matrix expected") if (epsilon <= 0) stop("epsilon must be positive") if (!is.euclid(dis)) stop("Euclidean property is expected for dis") D2 <- as.matrix(dis)^2/2 n <- dim(D2)[1] result <- data.frame(matrix(0, n, 4)) names(result) <- c("sim", "pro", "met", "num") relax <- 0 x0 <- apply(D2, 1, sum)/sum(D2) result$sim <- x0 objective0 <- t(x0) %*% D2 %*% x0 if (comment == TRUE) print("evolution of the objective function:") xk <- x0 COUNTER <- 0 repeat { if (COUNTER > maxit) stop("maxit reached") repeat { if (COUNTER > maxit) stop("maxit reached") maxi.temp <- t(xk) %*% D2 %*% xk if (comment == TRUE) print(as.character(maxi.temp)) deltaf <- (-2 * D2 %*% xk) sature <- (abs(xk) < epsilon) if (relax != 0) { sature[relax] <- FALSE relax <- 0 } yk <- (-deltaf) yk[sature] <- 0 yk[!(sature)] <- yk[!(sature)] - mean(yk[!(sature)]) if (max(abs(yk)) < epsilon) { break } alpha.max <- as.vector(min(-xk[yk < 0]/yk[yk < 0])) alpha.opt <- as.vector(-(t(xk) %*% D2 %*% yk)/(t(yk) %*% D2 %*% yk)) if ((alpha.opt > alpha.max) | (alpha.opt < 0)) { alpha <- alpha.max } else { alpha <- alpha.opt } if (abs(maxi.temp - t(xk + alpha * yk) %*% D2 %*% (xk + alpha * yk)) < epsilon) { break } xk <- xk + alpha * yk COUNTER <- COUNTER + 1 } if (prod(!sature) == 1) { if (comment == TRUE) print("KT") break } vectD2 <- D2 %*% xk u <- 2 * (mean(vectD2[!sature]) - vectD2[sature]) if (min(u) >= 0) { if (comment == TRUE) print("KT") break } else { if (comment == TRUE) print("relaxation") satu <- (1:n)[sature] relax <- satu[u == min(u)] relax <- relax[1] } COUNTER <- COUNTER + 1 } if (comment == TRUE) print(list(objective.init = objective0, objective.final = maxi.temp)) result$num <- as.vector(xk, mode = "numeric") result$num[result$num < epsilon] <- 0 xk <- x0/sqrt(sum(x0 * x0)) COUNTER <- 0 repeat { if (COUNTER > maxit) stop("maxit reached") yk <- D2 %*% xk yk <- yk/sqrt(sum(yk * yk)) if (max(xk - yk) > epsilon) { xk <- yk } else break COUNTER <- COUNTER + 1 } x0 <- as.vector(yk, mode = "numeric") result$pro <- x0/sum(x0) result$met <- x0 * x0 restot <- list() restot$value <- divc(cbind.data.frame(result$num), dis)[, 1] restot$vectors <- result return(restot) } ## -- code ends here On Wed, Jun 12, 2013 at 8:13 AM, Zoltan Botta-Dukat < botta-dukat.zol...@okologia.mta.hu> wrote: > Dear Thomas, > > > > >> Is it possible, and meaningful, to calculate Rao`s quadratic entropy >> index (function dbFD from package FD) with presence-absence data? >> > Yes, it is possible and meaningful. In this case it is simply the mean > distance between the occurring species. > > > > I don't get a error message but when I use the function with >> (...,scale.Rao=TRUE) the calculations never end. >> > I think this problem is not related to p/a data. Calculation may demand > long time for large datasets. > Anyway, I think scaling by the possible maximum is meaningless if you > don't have abundance data. > > Best wishes > > Zoltan > > > >> Thank you for your answer. >> >> /Thomas >> >> > > -- > Botta-Dukát Zoltán > ------------------------------**-- > Ãkológiai és Botanikai Intézet > Magyar Tudományos Akadémia > Ãkológiai Kutatóközpont > ------------------------------**-- > 2163. Vácrátót, Alkotmány u. 2-4. > tel: +36 28 360122/157 > fax: +36 28 360110 > botta-dukat.zoltan@okologia.**mta.hu <botta-dukat.zol...@okologia.mta.hu> > www.okologia.mta.hu > > > Zoltán BOTTA-Dukát > ------------------------------**-- > Institute of Ecology and Botany > Hungarian Academy of Sciences > Centre for Ecological Research > ------------------------------**-- > H-2163 Vácrátót, Alkomány u. 2-4. > HUNGARY > Phone: +36 28 360122/157 > Fax..: +36 28 360110 > botta-dukat.zoltan@okologia.**mta.hu <botta-dukat.zol...@okologia.mta.hu> > www.okologia.mta.hu > > > ______________________________**_________________ > R-sig-ecology mailing list > R-sig-ecology@r-project.org > https://stat.ethz.ch/mailman/**listinfo/r-sig-ecology<https://stat.ethz.ch/mailman/listinfo/r-sig-ecology> > > [[alternative HTML version deleted]]
_______________________________________________ R-sig-ecology mailing list R-sig-ecology@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
