Search code examples
rvegan

Errors with distance-decay using betapart and ddecay packages

My goal is to create a distance-decay curve for species data vs geographic distance. However, I am running into errors. For the betapart package, this may be due to the lack of columns relative to the number of rows. Is there a way to get past this? If not, is there another method for creating a distance-decay curve (and plotting it)? I also tried the ddecay package but ran into errors there too. Any help is much appreciated. Data is in structure form below.

# BETAPART -------------------------------------------------
library(betapart)
spat.dist<-dist(coords)
dissim.BCI<-beta.pair.abund(spec)$beta.bray.bal

plot(spat.dist, dissim.BCI, ylim=c(0,1), xlim=c(0, max(spat.dist)))

BCI.decay.exp<-decay.model(dissim.BCI, spat.dist, y.type="dissim", model.type="exp", perm=100)
#========================================================================================================

I also tried a few other packages --------------------------

# ddecay package -------------------------------------------
devtools::install_github("chihlinwei/ddecay")
the issue with this method is that it requires the use of a gradient however, I would like to avoid that if possible but I do not see a way around this. Also they do not include their example data in the package.

dd <- beta.decay(gradient=spat.dist, counts=decostand(spec, method="pa"),
             coords=coords, nboots=1000,
             dis.fun = "beta.pair", index.family = "sorensen", dis = 1, like.pairs=T)
x <- vegdist(coords, method = "euclidean")
y <- 1 - dist(decostand(spec, method="pa"), index.family = "sorensen")[[1]]
plot(x, y)
lines(dd$Predictions[, "x"], dd$Predictions[,"mean"], col="red", lwd=2)
#========================================================================================================

# DATA -----------------------------------------------------

    spec <- structure(list(Ccol = c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), Acol = c(0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0), NYcol = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0), Mcol = c(0, 
    0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0), AAcol = c(14, 0, 14, 3, 11, 1, 0, 2, 0, 
    3, 0, 4, 0, 1, 8, 2, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 7), 
        Ncol = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1), ATBcol = c(0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 20, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 3), CVcol = c(0, 0, 0, 0, 0, 0, 1, 20, 
        0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 2, 0, 0, 
        0, 6), AZNcol = c(0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), GBcol = c(0, 
        0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0), KHAcol = c(0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 
        0, 0, 0, 0), AFcol = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 1, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0), AFPcol = c(0, 
        0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 1), TIAcol = c(4, 1, 0, 2, 6, 0, 
        1, 1, 0, 2, 0, 0, 0, 1, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 1, 0), AUcol = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), AScol = c(0, 
        4, 0, 2, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 5, 0, 0), NSAcol = c(0, 0, 0, 0, 0, 0, 
        0, 0, 0, 7, 0, 0, 3, 0, 0, 0, 4, 0, 2, 0, 1, 0, 9, 5, 1, 
        0, 0, 2, 0), WZcol = c(0, 0, 0, 0, 0, 0, 1, 0, 0, 10, 4, 
        0, 0, 0, 0, 0, 0, 1, 5, 0, 0, 0, 17, 4, 0, 0, 0, 0, 0), AJcol = c(0, 
        3, 6, 0, 0, 1, 0, 4, 0, 0, 0, 0, 39, 12, 0, 0, 0, 0, 0, 0, 
        0, 4, 5, 1, 12, 13, 16, 0, 5), EADcol = c(4, 1, 2, 1, 2, 
        0, 0, 0, 0, 4, 0, 2, 1, 1, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 
        0, 0, 0, 0, 1), CAcol = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0), Pcol = c(0, 
        0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 60, 0, 0, 
        13, 0, 8, 1, 0, 0, 0, 0, 0), ASDcol = c(3, 5, 6, 17, 3, 5, 
        26, 2, 0, 17, 3, 10, 6, 3, 2, 4, 0, 0, 5, 25, 0, 0, 0, 2, 
        2, 9, 0, 2, 8), RMAcol = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 
        OUcol = c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), KAcol = c(0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 12, 
        0, 0, 0, 0, 0, 8, 1), PACcol = c(0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 11, 2, 0, 37, 0, 24, 
        1, 0, 0), LAAcol = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 
        0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0), GAcol = c(1, 
        0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 
        0, 0, 3, 0, 0, 0, 2, 0, 0), AAcol = c(1, 0, 1, 0, 0, 0, 0, 
        0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 1, 0), EVAcol = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0), EAcol = c(0, 
        0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0), AKcol = c(0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 
        0, 0, 0), Acol = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 1, 0), QAcol = c(0, 
        0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0), YAcol = c(11, 24, 21, 63, 44, 
        95, 12, 43, 0, 5, 26, 22, 25, 48, 86, 2, 0, 0, 13, 0, 0, 
        2, 0, 0, 60, 6, 7, 0, 45), BANcol = c(0, 0, 0, 3, 0, 0, 0, 
        0, 0, 0, 0, 0, 24, 0, 6, 0, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 
        9, 17, 17), VCcol = c(0, 38, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), Vcol = c(0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 1, 0, 0, 0, 0, 0, 0), Ocol = c(0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0), AVcol = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), JXcol = c(0, 
        3, 3, 0, 0, 0, 0, 0, 8, 0, 0, 10, 3, 0, 0, 5, 0, 0, 0, 1, 
        0, 0, 0, 2, 4, 1, 0, 0, 0)), class = "data.frame", row.names = c(NA, 
    -29L))

   coords <- structure(list(Lat.x = c(34.43363, 34.36784, 34.32587, 34.19891, 
    34.24217, 34.24863, 34.18137, 34.16838, 34.10961, 34.08329, 34.40571, 
    34.39591, 34.39292, 34.37466, 34.28948, 34.26146, 34.04687, 34.0409, 
    34.068339, 34.34679, 34.17161, 34.23308, 34.21544, 34.14922, 
    34.27539, 34.2323, 34.19057, 34.07042, 34.06289), Lon.x = c(-94.94494, 
    -94.92512, -94.94429, -94.84497, -94.8573, -94.85641, -94.887, 
    -94.91322, -94.92913, -94.93276, -95.02622, -95.04382, -94.96295, 
    -94.83733, -94.81071, -94.79161, -95.03968, -95.0608, -95.086986, 
    -95.03345, -95.23862, -95.25619, -95.1041, -95.02286, -95.02672, 
    -95.02626, -95.02941, -95.01746, -94.98786)), class = "data.frame", row.names = c(NA, 
    -29L))
Solution

  • You can get more answers, if you tell what was the problem. For instance, which functions failed and what was the error message. I had a look at betapart::decay.model(), where I could get this error message:

    Error in eval(family$initialize) : 
  cannot find valid starting values: please specify some

    I cut the long story short: you cannot use this function with your data because you have dissimilarities of 1 in your data, dissimilarities are turned into similarities with 1-dissimilarity and this makes these values zero similarities (that is, these pairs of sampling unit have nothing in common, they share no species). Function decay.model uses glm with gaussian family with log-link, and log-link requires that you give the starting values, if you have zeros in the y-variate.

    I think that you have four alternatives:

    1. You do not use the method as it does not suit your data.
    2. You modify the decay.model function so that you can specify the starting values, like the error message suggested. This means that you add mustart to the function call so that it reads, e.g., glm(y ~ x, family=gaussian(link="log"), mustart=pmax(y, 0.01)). This replaces zeros with 0.01 as starting values.
    3. You change maximum distances from 1 to something smaller, for instance, 0.99: dissim.BCI[dissim.BCI==1] <- 0.99. However, this changes the data, and also changes the results from those you get with alternative 2 (which only changes starting values, but data are unmodified). However, the effect is not very large and any Bayesian would claim that dissimilarity 1 is just a frequentist folly (you just haven't seen the case that is in common with these sampling units).
    4. You change the maximum distance to missing values. This will change data more than alternative 3. It removes maximum dissimilarities and these no longer influence the decay curve. The effect is the same as censoring greatest dissimilarities. The results change more than in alternative 3.