Need Help Making an Ordihull

Question

I have been collaborating on this code that creates an NMDS plot and I want to add shaded polygons of the points. However, the ordihull code keeps returning the following error. Why would the argument be of length zero?

Error in if (n < 4) return(colMeans(x[-n, , drop = FALSE])) : argument is of length zero

> m1 <- metaMDS(d1)
> m2 <- metaMDS(d2)
> m3 <- metaMDS(d3)

> mdat <- data.frame(m3$points)
> mdat$site <- substr(rownames(mdat), 1, 1) mdat$col <- ifelse(mdat$site == "D", "red", 
               ifelse(mdat$site == "H", "blue", "green"))

> plot(mdat[,1], mdat[,2], pch=16, col=mdat$col, display = "sites",
 xlab="NMDS1", ylab="NMDS2", xlim=c(-0.2, 0.2), 
 ylim=c(-0.2, 0.2), main= "Phylum")

> ordihull(mdat[,1], mdat[,2], display="sites", label=T,
     lwd=2, draw="polygon",col= c("blue", "red", "green"))

---

Here is the Dput:

> structure(list(p__Proteobacteria = c(44.807, 40.907, 36.558,36.811, 
 39.401, 40.114, 45.911, 43.133, 30.137, 27.734, 26.722,
 31.261), p__Actinobacteria = c(26.819, 34.651, 40.904, 38.847,
 39.446, 37.523, 29.881, 29.251, 31.783, 23.641, 34.918, 31.308
 ), p__Acidobacteria = c(8.48, 6.6, 5.934, 6.609, 5.89, 7.567,
 5.795, 6.666, 10.616, 10.709, 8.988, 11.794), p__Bacteroidetes = 
 c(7.56, 8.189, 5.363, 6.223, 4.716, 3.613, 4.65, 5.2, 4.281, 2.785, 
 2.808, 3.271), p__Gemmatimonadetes = c(3.529, 2.108, 1.213, 1.193, 
 1.541, 1.439, 1.006, 1.171, 5.794, 4.107, 4.001, 2.747), 
 p__Chloroflexi = c(2.686, 2.987, 2.979, 3.049, 4.128, 4.564, 5.304, 
 4.624, 3.669, 2.775, 4.534, 4.94), p__Bacteria_unclassified = 
 c(2.38, 1.869, 1.579, 1.247, 2.3, 2.108, 1.36, 1.193, 3.126, 1.885, 
 2.987, 2.37), p__Firmicutes = c(0.998, 0.807, 2.76, 2.962, 0.866, 
 1.32, 1.651, 2.073, 1.099, 1.046, 1.3, 1.302), p__Verrucomicrobia = 
 c(0.676, 0.404, 0.32, 0.35, 0.293, 0.239, 0.188, 0.261, 0.521, 
 0.726, 0.52, 0.397), p__Nitrospirae = c(0.464, 0.244, 0.198, 0.208, 
 0.016, 0.032, 0.024, 0.042, 0.296, 0.103, 0.229, 0.211), 
 p__Candidatus_Saccharibacteria = c(0.421, 0.511, 0.456, 0.552, 
 0.523, 0.6, 0.842, 1.016, 0.672, 0.636, 0.465, 0.736), 
 p__Planctomycetes = c(0.392, 0.267, 0.354, 0.285, 0.275, 0.356, 
 0.285, 0.276, 0.33, 0.438, 0.552, 0.365), p__Fibrobacteres = c(0.14,
 0.074, 0.007, 0.009, 0.072, 0.044, 0.136, 0.079, 0.117, 0.018,
 0.167, 0.065), p__Candidatus_Latescibacteria = c(0.113, 0.059,
 0.017, 0.005, 0.004, 0.017, 0.015, 0.009, 0, 0.011, 0.007, 0.018
 ), p__Latescibacteria = c(0.085, 0.04, 0.01, 0.004, 0.012, 0.015,
 0.033, 0.015, 0.012, 0.016, 0.011, 0.018), p__Cyanobacteria = 
 c(0.079, 0.048, 1.071, 1.372, 0.32, 0.19, 2.629, 4.689, 7.133, 
 22.963, 11.417, 8.767), p__Thermodesulfobacteria = c(0.068, 0.057, 
 0.115, 0.103, 0.008, 0.01, 0.015, 0.007, 0.01, 0.003, 0.002, 0.013),
 p__Elusimicrobia = c(0.059, 0.021, 0.012, 0.001, 0.004, 0.002, 
 0.015, 0.017, 0, 0.002, 0.005, 0.006), p__Chlorobi = c(0.052,
 0.025, 0.002, 0.012, 0.029, 0.046, 0.033, 0.04, 0.05, 0.02,
 0.046, 0.025), p__Armatimonadetes = c(0.046, 0.053, 0.051,
 0.072, 0.076, 0.095, 0.048, 0.053, 0.197, 0.159, 0.128, 0.125
 ), p__Spirochaetes = c(0.035, 0.021, 0.002, 0.001, 0, 0.002,
 0.024, 0.039, 0, 0, 0, 0), p__Parcubacteria = c(0.03, 0.013,
 0, 0, 0.01, 0.015, 0.042, 0.037, 0.032, 0.059, 0.053, 0.011
 ), p__Chlamydiae = c(0.028, 0.017, 0.046, 0.05, 0.014, 0.007,
 0.021, 0.022, 0.07, 0.074, 0.08, 0.152)), class = "data.frame", 
 row.names = c("D15B", "D610B", "D15F", "D610F", "HR15B", "HR610B", 
 "HR15F", "HR610F", "C15B", "C610B", "C15F", "C610F"))

Here are the codes:

> phylum.dat <- dput

> x <- data.frame(tax=names(phylum.dat), nsites=apply(phylum.dat, 2, function(x){length(which(x>0))}))

> d1 <- vegdist(phylum.dat, method = "jaccard", binary = TRUE)
> d2 <- vegdist(log1p(phylum.dat, method = "jaccard"))

> logit_phylum <- as.matrix(phylum.dat+1)/100
> d3 <- qlogis(logit_phylum)
> d3 <- d3+abs(min(d3))
> d3 <- vegdist(d3, method = "jaccard")

> m1 <- metaMDS(d1)
> m2 <- metaMDS(d2)
> m3 <- metaMDS(d3)

> e1 <- envfit(m3, phylum.dat)
> exy <- data.frame(tax=names(phylum.dat),
> x=e1$vectors$arrows[,1],
> y=e1$vectors$arrows[,2],
> pval=e1$vectors$pvals,
> r=e1$vectors$r)
> rownames(exy) <- NULL

> exy <- exy[order(-exy$r),]

> mdat <- data.frame(m3$points)
> mdat$site <- substr(rownames(mdat), 1, 1)
> mdat$col <- ifelse(mdat$site == "D", "red",
> ifelse(mdat$site == "H", "blue", "green"))

> mdat$rad <- sqrt((mdat$MDS1^2) + (mdat$MDS2^2))
> max(mdat$rad)

> exy$x2 <- 0.17 * exy$r * exy$x
> exy$y2 <- 0.17 * exy$r * exy$y
> exy$adj <- ifelse(exy$x < 0, 1, 0)

> plot(mdat[,1], mdat[,2], pch=16, col=mdat$col,
> xlab="NMDS1", ylab="NMDS2", xlim=c(-0.2, 0.2),
> ylim=c(-0.2, 0.2), main= "Phylum")

> ordihull(mdat[,1], mdat[,2], display="sites", label=T,
> lwd=2, draw="polygon",col= c("blue", "red", "green"))