paired.x <- rbind(splitx[[group.combinations[1,i]]],
splitx[[group.combinations[2,i]]])
# remove empty species
ifzero <- which(apply(paired.x, 2, sum) == 0)
if(length(ifzero > 0)){
paired.x <- paired.x[,-which(colSums(paired.x)==0)]}
if(length(which(rowSums(paired.x) == 0)) != 0){stop("ERROR : A row/sample is empty")}
group.pair <- c(rep(group.combinations[1,i], NROW(splitx[[group.combinations[1,i]]])),
rep(group.combinations[2,i], NROW(splitx[[group.combinations[2,i]]])))
ddelta <- rbind(ddelta, DNCImper:::DNCI.ses(x=paired.x,grouping=group.pair,id=id, Nperm = Nperm,
count = count, plotSIMPER = plotSIMPER, dataTYPE = dataTYPE, parallelComputing = parallelComputing)) #here is the part that calculates the index based on PERSIMPER
if(symmetrize == TRUE)
{
Add <- which(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]])) == max(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]]))))
if(Add == 1)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[1,i]]][,1]),
length(splitx[[group.combinations[2,i]]][,1]))
splitx[[group.combinations[1,i]]] <- splitx[[group.combinations[1,i]]][sampled_lines,]
}
if(Add == 2)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[2,i]]][,1]),
length(splitx[[group.combinations[1,i]]][,1]))
splitx[[group.combinations[2,i]]] <- splitx[[group.combinations[2,i]]][sampled_lines,]
}
}
i=3
if(symmetrize == TRUE)
{
Add <- which(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]])) == max(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]]))))
if(Add == 1)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[1,i]]][,1]),
length(splitx[[group.combinations[2,i]]][,1]))
splitx[[group.combinations[1,i]]] <- splitx[[group.combinations[1,i]]][sampled_lines,]
}
if(Add == 2)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[2,i]]][,1]),
length(splitx[[group.combinations[1,i]]][,1]))
splitx[[group.combinations[2,i]]] <- splitx[[group.combinations[2,i]]][sampled_lines,]
}
}
symmetrize
if(symmetrize == TRUE)
{
Add <- which(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]])) == max(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]]))))
if(Add == 1)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[1,i]]][,1]),
length(splitx[[group.combinations[2,i]]][,1]))
splitx[[group.combinations[1,i]]] <- splitx[[group.combinations[1,i]]][sampled_lines,]
}
if(Add == 2)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[2,i]]][,1]),
length(splitx[[group.combinations[1,i]]][,1]))
splitx[[group.combinations[2,i]]] <- splitx[[group.combinations[2,i]]][sampled_lines,]
}
}
i=2
if(symmetrize == TRUE)
{
Add <- which(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]])) == max(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]]))))
if(Add == 1)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[1,i]]][,1]),
length(splitx[[group.combinations[2,i]]][,1]))
splitx[[group.combinations[1,i]]] <- splitx[[group.combinations[1,i]]][sampled_lines,]
}
if(Add == 2)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[2,i]]][,1]),
length(splitx[[group.combinations[1,i]]][,1]))
splitx[[group.combinations[2,i]]] <- splitx[[group.combinations[2,i]]][sampled_lines,]
}
}
i=1
if(symmetrize == TRUE)
{
Add <- which(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]])) == max(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]]))))
if(Add == 1)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[1,i]]][,1]),
length(splitx[[group.combinations[2,i]]][,1]))
splitx[[group.combinations[1,i]]] <- splitx[[group.combinations[1,i]]][sampled_lines,]
}
if(Add == 2)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[2,i]]][,1]),
length(splitx[[group.combinations[1,i]]][,1]))
splitx[[group.combinations[2,i]]] <- splitx[[group.combinations[2,i]]][sampled_lines,]
}
}
x=comm_mat
grouping=group$Group
group.combinations <- combn(unique(sort(grouping)),2)
group.combinations
i=1
symmetrize
symmetrize=TRUE
symmetrize
if(symmetrize == TRUE)
symmetrize == TRUE
if(symmetrize == TRUE)
{
Add <- which(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]])) == max(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]]))))
if(Add == 1)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[1,i]]][,1]),
length(splitx[[group.combinations[2,i]]][,1]))
splitx[[group.combinations[1,i]]] <- splitx[[group.combinations[1,i]]][sampled_lines,]
}
if(Add == 2)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[2,i]]][,1]),
length(splitx[[group.combinations[1,i]]][,1]))
splitx[[group.combinations[2,i]]] <- splitx[[group.combinations[2,i]]][sampled_lines,]
}
}
i=2
if(symmetrize == TRUE)
{
Add <- which(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]])) == max(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]]))))
if(Add == 1)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[1,i]]][,1]),
length(splitx[[group.combinations[2,i]]][,1]))
splitx[[group.combinations[1,i]]] <- splitx[[group.combinations[1,i]]][sampled_lines,]
}
if(Add == 2)
{
sampled_lines <- sample(1:length(splitx[[group.combinations[2,i]]][,1]),
length(splitx[[group.combinations[1,i]]][,1]))
splitx[[group.combinations[2,i]]] <- splitx[[group.combinations[2,i]]][sampled_lines,]
}
}
Add <- which(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]])) == max(c(NROW(splitx[[group.combinations[1,i]]]),
NROW(splitx[[group.combinations[2,i]]]))))
Add
NROW(splitx[[group.combinations[1,i]]])
NROW(splitx[[group.combinations[2,i]]]))
NROW(splitx[[group.combinations[2,i]]])
n_reps <- 100
dnci_values <- numeric(n_reps)
dnci_values
library(MVNH)
??MVNH_det()
??MVNH_dissimilarity()
??MVNH_det()
var.partition <- function(metacomm_tsdata){
## The function "var.partition" performs the partitioning of variability
## across hierarchical levesl within a metacommunity.
## The input array "metacomm_tsdata" is an N*T*M array. The first dimension represents N species,
## the second represents time-series observations of length T, and the third represents M local communities.
## The output includes four variability and four synchrony metrics as defined in the main text.
## Note that, to be able to handle large metacommunities, this code has avoided calculating all covariance.
ts_metacom <- apply(metacomm_tsdata,2,sum) #summing the total biovolume across all species and patch at every time point
ts_patch <- apply(metacomm_tsdata,c(2,3),sum) #summing the total biovolume across all species in every patch at every time point
ts_species <- apply(metacomm_tsdata,c(1,2),sum)#summing the total biovolume across all patches for every species at every time point
sd_metacom <- sd(ts_metacom)
sd_patch_k <- apply(ts_patch,2,sd)
sd_species_i <- apply(ts_species,1,sd)
sd_species_patch_ik <- apply(metacomm_tsdata,c(1,3),sd)
mean_metacom <- mean(ts_metacom)
CV_S_L <- sum(sd_species_patch_ik)/mean_metacom
CV_C_L <- sum(sd_patch_k)/mean_metacom
CV_S_R <- sum(sd_species_i)/mean_metacom
CV_C_R <- sd_metacom/mean_metacom
phi_S_L2R <- CV_S_R/CV_S_L
phi_C_L2R <- CV_C_R/CV_C_L
phi_S2C_L <- CV_C_L/CV_S_L
phi_S2C_R <- CV_C_R/CV_S_R
partition_3level <- c(CV_S_L=CV_S_L, CV_C_L=CV_C_L, CV_S_R=CV_S_R, CV_C_R=CV_C_R,
phi_S_L2R=phi_S_L2R, phi_C_L2R=phi_C_L2R, phi_S2C_L=phi_S2C_L,
phi_S2C_R=phi_S2C_R)
return(partition_3level)
}
#The master funciton to be benchmarked
CV_simple_function <- function(species, time, plot, abundance){
# Extract unique values for Species, Sampling_date_order, and plot
species_vals <- unique(species)
date_vals <- unique(time)
plot_vals <- unique(plot)
# Extract unique values for Species, Sampling_date_order, and plot
species_vals <- unique(species)
date_vals <- unique(time)
plot_vals <- unique(plot)
# Create the array with dimensions based on unique Species, Sampling_date_order, and plot
metacomm_tsdata <- array(0, dim = c(length(species_vals), length(date_vals), length(plot_vals)))
# Populate the array with values from the data frame
for(i in 1:nrow(df)){
# Map Species, Sampling_date_order, and plot to their corresponding indices
species_index <- which(species_vals == species[i])
date_index <- which(date_vals == time[i])
plot_index <- which(plot_vals == plot[i])
# Assign the value from the data frame (Abundance) to the array
metacomm_tsdata[species_index, date_index, plot_index] <- abundance[i]
}
result <- var.partition(metacomm_tsdata)
return(result)
}
var.partition <- function(metacomm_tsdata){
## The function "var.partition" performs the partitioning of variability
## across hierarchical levesl within a metacommunity.
## The input array "metacomm_tsdata" is an N*T*M array. The first dimension represents N species,
## the second represents time-series observations of length T, and the third represents M local communities.
## The output includes four variability and four synchrony metrics as defined in the main text.
## Note that, to be able to handle large metacommunities, this code has avoided calculating all covariance.
ts_metacom <- apply(metacomm_tsdata,2,sum) #summing the total biovolume across all species and patch at every time point
ts_patch <- apply(metacomm_tsdata,c(2,3),sum) #summing the total biovolume across all species in every patch at every time point
ts_species <- apply(metacomm_tsdata,c(1,2),sum)#summing the total biovolume across all patches for every species at every time point
sd_metacom <- sd(ts_metacom)
sd_patch_k <- apply(ts_patch,2,sd)
sd_species_i <- apply(ts_species,1,sd)
sd_species_patch_ik <- apply(metacomm_tsdata,c(1,3),sd)
mean_metacom <- mean(ts_metacom)
CV_S_L <- sum(sd_species_patch_ik)/mean_metacom
CV_C_L <- sum(sd_patch_k)/mean_metacom
CV_S_R <- sum(sd_species_i)/mean_metacom
CV_C_R <- sd_metacom/mean_metacom
phi_S_L2R <- CV_S_R/CV_S_L
phi_C_L2R <- CV_C_R/CV_C_L
phi_S2C_L <- CV_C_L/CV_S_L
phi_S2C_R <- CV_C_R/CV_S_R
partition_3level <- c(CV_S_L=CV_S_L, CV_C_L=CV_C_L, CV_S_R=CV_S_R, CV_C_R=CV_C_R,
phi_S_L2R=phi_S_L2R, phi_C_L2R=phi_C_L2R, phi_S2C_L=phi_S2C_L,
phi_S2C_R=phi_S2C_R)
return(partition_3level)
}
#The master funciton to be benchmarked
CV_simple_function <- function(species, time, plot, abundance){
# Extract unique values for Species, Sampling_date_order, and plot
species_vals <- unique(species)
date_vals <- unique(time)
plot_vals <- unique(plot)
# Create the array with dimensions based on unique Species, Sampling_date_order, and plot
metacomm_tsdata <- array(0, dim = c(length(species_vals), length(date_vals), length(plot_vals)))
# Populate the array with values from the data frame
for(i in 1:nrow(df)){
# Map Species, Sampling_date_order, and plot to their corresponding indices
species_index <- which(species_vals == species[i])
date_index <- which(date_vals == time[i])
plot_index <- which(plot_vals == plot[i])
# Assign the value from the data frame (Abundance) to the array
metacomm_tsdata[species_index, date_index, plot_index] <- abundance[i]
}
result <- var.partition(metacomm_tsdata)
return(result)
}
df<-read.csv("/Users/yc2864/Downloads//Users/yc2864/Downloads/model_01_03_df_dispersal_0.00010051545923590379_niche_0.034879574602545566.csv")
df<-read.csv("/Users/yc2864/Downloads/model_01_03_df_dispersal_0.00010051545923590379_niche_0.034879574602545566.csv")
data<-df%>%filter(Rep==5)
data
CV_simple_function(data$Species, data$Sampling_date_order, data$plot, data$Abundance)
species_vals <- unique(data$Species)
date_vals <- unique(data$Time)
plot_vals<- unique(data$Patch)
metacomm_tsdata <- array(0, dim = c(length(species_vals), length(date_vals), length(plot_vals)))
for(i in 1:nrow(df)){
# Map Species, Sampling_date_order, and plot to their corresponding indices
species_index <- which(species_vals == species[i])
date_index <- which(date_vals == time[i])
plot_index <- which(plot_vals == plot[i])
# Assign the value from the data frame (Abundance) to the array
metacomm_tsdata[species_index, date_index, plot_index] <- abundance[i]
}
data<-df%>%filter(Rep==5)%>%
filter(Species == 212)
data<-df%>%filter(Rep==5)%>%
filter(Species == 212, Species == 236, Species == 289, Species == 426, Patch ==8)
data<-df%>%filter(Rep==5)%>%
filter(Species == 212, Species == 236, Species == 289, Species == 426,
Patch == 8, Patch == 11, Patch == 12, Patch == 19, Patch == 26,
Patch == 28, Patch == 31, Patch == 32, Patch == 20)
result<-CV_simple_function(data$Species, data$Sampling_date_order, data$plot, data$Abundance)
result<-CV_simple_function(data$Species, data$Sampling_date_order, data$plot, data$Abundance)
result
data<-df%>%filter(Rep==5)%>%
filter(Species == 212, Species == 236, Species == 289, Species == 426,
Patch == 8, Patch == 11, Patch == 12, Patch == 19, Patch == 26,
Patch == 28, Patch == 31, Patch == 32, Patch == 20)
data
data<-df%>%filter(Rep==5)
data
df%>%filter(Rep==5)%>%
filter(Species == 212, Species == 236, Species == 289, Species == 426)
df%>%filter(Rep==5)%>%
filter(Species == 212)
df%>%filter(Rep==5)%>%
filter(Species == 212, Species == 236)
df%>%filter(Rep==5)%>%
filter(Species == 212 | Species == 236)
df%>%filter(Rep==5)%>%
filter(Species == 212 | Species == 236 | Species == 289)
data<-df%>%filter(Rep==5)%>%
filter(Species == 212 | Species == 236 | Species == 289 | Species == 426 |
Patch == 8 | Patch == 11 | Patch == 12 | Patch == 19 | Patch == 26 |
Patch == 28 | Patch == 31 | Patch == 32 | Patch == 20)
data
data<-df%>%filter(Rep==5)%>%
filter(Species == 212 | Species == 236 | Species == 289 | Species == 426)%>%
filter(Patch == 8 | Patch == 11 | Patch == 12 | Patch == 19 | Patch == 26 |
Patch == 28 | Patch == 31 | Patch == 32 | Patch == 20)
data
853+83
result<-CV_simple_function(data$Species, data$Sampling_date_order, data$plot, data$Abundance)
result
species_vals <- unique(data$Species)
date_vals <- unique(data$Time)
plot_vals <- unique(data$Patch)
plot_vals <- unique(plot)
plot_vals <- unique(plot)
species_vals
date_vals
plot_vals
plot_vals <- unique(data$Patch)
plot_vals
metacomm_tsdata <- array(0, dim = c(length(species_vals), length(date_vals), length(plot_vals)))
abundance<-data$N
abundance
metacomm_tsdata
metacomm_tsdata <- array(0, dim = c(length(species_vals), length(date_vals), length(plot_vals)))
metacomm_tsdata
for(i in 1:nrow(df)){
# Map Species, Sampling_date_order, and plot to their corresponding indices
species_index <- which(species_vals == species[i])
date_index <- which(date_vals == time[i])
plot_index <- which(plot_vals == plot[i])
# Assign the value from the data frame (Abundance) to the array
metacomm_tsdata[species_index, date_index, plot_index] <- abundance[i]
}
result <- var.partition(metacomm_tsdata)
result
metacomm_tsdata
abundance
plot_index
time
i
i=1
1:nrow(data)
i
species_vals == species[i]
species[i]
species
species=data$Species
time=data$Time
plot=data$Patch
species_vals
date_vals
plot_vals
specie
species
metacomm_tsdata <- array(0, dim = c(length(species_vals), length(date_vals), length(plot_vals)))
for(i in 1:nrow(data)){
# Map Species, Sampling_date_order, and plot to their corresponding indices
species_index <- which(species_vals == species[i])
date_index <- which(date_vals == time[i])
plot_index <- which(plot_vals == plot[i])
# Assign the value from the data frame (Abundance) to the array
metacomm_tsdata[species_index, date_index, plot_index] <- abundance[i]
}
result <- var.partition(metacomm_tsdata)
result
result<-DNCImper:::DNCI_multigroup(comm_mat,
group$Group, Nperm = 100,
symmetrize = TRUE,
plotSIMPER = FALSE)
result<-DNCImper:::DNCI_multigroup(comm_mat,
group$Group, Nperm = 10,
symmetrize = FALSE,
plotSIMPER = FALSE)
result
result[,2:3]
result[,2:4]
rbind(result[,2:4], result[,2:4])
dnci_values_df<-[]
dnci_values_df<-()
dnci_values_df<-()
dnci_values_df<-data.frame()
rbind(dnci_values_df, result[,2:4])
dnci_values_df<-data.frame()
rbind(dnci_values_df, result[,2:4])
rbind(dnci_values_df, result[,2:4])
rbind(dnci_values_df, result[,2:4])
rbind(dnci_values_df, result[,2:4])
rbind(dnci_values_df, result[,2:4])
dnci_values_df<-rbind(dnci_values_df, result[,2:4])
dnci_values_df<-rbind(dnci_values_df, result[,2:4])
dnci_values_df
n_reps <- 2
dnci_values_df<-data.frame()
set.seed(123)  # For reproducibility
for (i in 1:n_reps) {
result <- DNCImper:::DNCI_multigroup(comm_mat,
group$Group, Nperm = 10,
symmetrize = FALSE,
plotSIMPER = FALSE)
# Save the DNCI values for all first group comparison
dnci_values[i] <- result[,2:4]
dnci_values_df<-rbind(dnci_values_df, dnci_values[i])
}
dnci_values_df
n_reps <- 2
dnci_values_df<-data.frame()
set.seed(123)  # For reproducibility
for (i in 1:n_reps) {
result <- DNCImper:::DNCI_multigroup(comm_mat,
group$Group, Nperm = 5,
symmetrize = FALSE,
plotSIMPER = FALSE)
# Save the DNCI values for all first group comparison
dnci_values_df<-rbind(dnci_values_df, result[,2:4] )
}
dnci_values_df
n_reps <- 2
dnci_values_df<-data.frame()
set.seed(123)  # For reproducibility
for (i in 1:n_reps) {
result <- DNCImper:::DNCI_multigroup(comm_mat,
group$Group, Nperm = 5,
symmetrize = FALSE,
plotSIMPER = FALSE)
# Save the DNCI values for all first group comparison
dnci_values_df<-rbind(dnci_values_df, result[,2:4] )
}
dnci_values_df
library(data.table)
library(tidyverse)
library(bench)
df <- read.csv("~/.julia/dev/MetaCommunityMetrics/data/metacomm_rodent_df.csv",na.strings = c(""), stringsAsFactors = FALSE) #There is a species called "NA", need to handle NAs with caution
df <- read.csv("~/Documents/research/MetaCommunityMetrics.jl/data/metacomm_rodent_df.csv",na.strings = c(""), stringsAsFactors = FALSE) #There is a species called "NA", need to handle NAs with caution
df
unique(df$Species)
comm_df_for_DNCI<-read.csv("~/Documents/research/MetaCommunityMetrics.jl/benchmarks/benchmark_r/data/DNCI_comm_t50.csv")
grouping_for_DNCI<-read.csv("~/Documents/research/MetaCommunityMetrics.jl/benchmarks/benchmark_r/data/cluster_list_t50.csv")
community_matrix<-df %>%
select(-Presence) %>%
pivot_wider(names_from = Species, values_from = Abundance, values_fill=0) %>%
select(-c(Year, Month, Day, Sampling_date_order, plot, Longitude, Latitude)) %>%
select(which(colSums(.) !=0)) %>%
filter(rowSums(.) != 0)
matrix_with_abundance <- df %>%
filter(Sampling_date_order == 50) %>%
select(-Presence) %>%
pivot_wider(names_from = Species, values_from = Abundance, values_fill=0) %>%
select(-c(Year, Month, Day, Sampling_date_order, plot, Longitude, Latitude, normalized_temperature, normalized_precipitation)) %>%
select(which(colSums(.) !=0)) %>%
filter(rowSums(.) != 0)
### The binary matrix
matrix_with_presence <- df %>%
filter(Sampling_date_order == 50) %>%
select(-Abundance) %>%
pivot_wider(names_from = Species, values_from = Presence, values_fill=0) %>%
select(-c(Year, Month, Day, Sampling_date_order, plot, Longitude, Latitude)) %>%
select(which(colSums(.) !=0)) %>%
filter(rowSums(.) != 0)
beta_diversity_1 <- mark(beta.div.comp(matrix_with_abundance, coef = "J", quant = TRUE),
iterations = 10000,
check = TRUE,
time_unit = "us")
library(bench)
install.packages("bench")
beta_diversity_1 <- mark(beta.div.comp(matrix_with_abundance, coef = "J", quant = TRUE),
iterations = 10000,
check = TRUE,
time_unit = "us")
library(bench)
beta_diversity_1 <- mark(beta.div.comp(matrix_with_abundance, coef = "J", quant = TRUE),
iterations = 10000,
check = TRUE,
time_unit = "us")
library(adespatial)
library(DNCImper)
beta_diversity_1 <- mark(beta.div.comp(matrix_with_abundance, coef = "J", quant = TRUE),
iterations = 10000,
check = TRUE,
time_unit = "us")
beta_diversity_1
beta_diversity_1 $time
library(profvis)
beta_diversity_1_prof <- profvis({
beta.div.comp(matrix_with_abundance, coef = "J", quant = TRUE)
})
beta_diversity_1_prof
beta_diversity_1_prof$x
full_df <- read.csv("~/Documents/research/MetaCommunityMetrics.jl/data/metacomm_rodent_df.csv",na.strings = c(""), stringsAsFactors = FALSE) #There is a species called "NA", need to handle NAs with caution
comm_full_df<-read.csv("~/Documents/research/MetaCommunityMetrics.jl/data/data_for_testing/DNCI_comm_t1_full_df.csv")
groups_full_df<-read.csv("~/Documents/research/MetaCommunityMetrics.jl/data/data_for_testing/cluster_list_t1_full_df.csv")
community_matrix<-df %>%
select(-Presence) %>%
pivot_wider(names_from = Species, values_from = Abundance, values_fill=0) %>%
select(-c(Year, Month, Day, Sampling_date_order, plot, Longitude, Latitude)) %>%
select(which(colSums(.) !=0)) %>%
filter(rowSums(.) != 0)
community_matrix
