03.2_PhysiologySIBER.Rmd

---
title: "SIBER Analysis of Coral Physiology"
author: "Peter Flood and Serena Hackerott"
date: "2/8/2024"
output:
  html_notebook:
    toc: yes
    toc_float: yes
  html_document:
    toc: yes
    df_print: paged
---

# Setup

```{r Setup, include = FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
getwd()
```


### Load Packages
```{r}
##Install Packages if Needed
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("cowplot")) install.packages("cowplot")
if (!require("SIBER")) install.packages("SIBER")
if (!require("vegan")) install.packages("vegan")
if (!require("Hmisc")) install.packages("Hmisc")
if (!require("corrplot")) install.packages("corrplot")
if (!require("tidyverse")) install.packages("tidyverse")
if (!require("Rmisc")) install.packages("Rmisc")
if (!require("ggpubr")) install.packages("ggpubr")
if (!require("lme4")) install.packages("lme4")
if (!require("lmerTest")) install.packages("lmerTest")
if (!require("DHARMa")) install.packages("DHARMa")
if (!require("emmeans")) install.packages("emmeans")

##Load Packages
library(ggplot2) #Required for ggplots
library(cowplot) #Required for plotting panel figures
library(SIBER) #Required for generating physiological ellipses and percent overlap calculations
library(vegan) #Required for NMDS
library(Hmisc) #Required for correlations
library(corrplot) #Required for correlation plot
library(tidyverse) #Required for data manipulation
library(Rmisc) #Required for data summary
library(ggpubr) #Required for stat_conf_ellipse
library(lme4) #Required for mixed effects models
library(lmerTest) #Required for significance in mixed effects models
library(DHARMa) #Required to check residuals of mixed effects models
library(emmeans) #Required for pairwise comparisons 

```


# Sample Data and Metadata

### Load and Organize Data
```{r}
##Load Data
#Note: Physiological metrics calculated in 02_PhysiologyMetrics.R file
Phys<-read.csv("Outputs/CoralData.csv", header=TRUE)

##Set factor variables
Phys$TimeP<-factor(Phys$TimeP, levels=c("TP1", "TP2", "TP3", "TP4"))
Phys$Site<-factor(Phys$Site, levels=c("KL", "SS"))
Phys$Genotype<-factor(Phys$Genotype, levels=c("AC8", "AC10", "AC12"))
Phys$Orig<-factor(Phys$Orig, levels=c("N", "T"))
Phys$Origin<-factor(Phys$Origin, levels=c("Native", "Transplant"))
Phys$Site.Orig<-factor(Phys$Site.Orig, levels=c("KL.N", "KL.T", "SS.N", "SS.T"))

##Set rownames
rownames(Phys)<-Phys$ID

##Load Environmental Data
#Note: Averaged environmental metrics calculated in 01_ExperimentalSetup.R file
Env<-read.csv("Outputs/AvgEnvData.csv", header=TRUE)
Env$TimeP<-factor(Env$TimeP, levels=c("TP1", "TP2", "TP3", "TP4"))
Env$Site<-factor(Env$Site, levels=c("KL", "SS"))
```


# SIBER Ellipses

### Prepare Input Data

#### NMDS Ordination
SIBER runs on a two dimensional data set. So first, we need to reduce the dimensionality of the physiology data. We're going to do this via Non-Metric Multidimensional Scaling (NMDS) on the scaled physiological metrics of interest.
```{r}
##Remove NA's
Phys.rm<-na.omit(Phys)

##Standardize Variables for similar scales
Phys.St <- Phys.rm %>% mutate_at(c("TP_ug.cm2", "AFDW_mg.cm2", "Chl_ug.cm2"), ~(scale(.) %>% as.vector))

#Select only the physiology columns of interest (without covariates) to put into the NMDS
names(Phys.St)
phys.num <- Phys.St[,c(11:13)]

#Run NMDS
phys.nmds <- metaMDS(phys.num, distance = "euclidean", k = 2, autotransform = F)

#Quick plot of NMDS output
plot(phys.nmds)
```


#### NMDS Model
```{r}
##Correlation with Physiology Metrics
phys.envi.fit <- envfit(phys.nmds, phys.num)
phys.envi.fit
phys.envi.vectors <- phys.envi.fit[["vectors"]][["arrows"]]

##Check Stress
stressplot(phys.nmds)
NMDS.stress<- sprintf("Stress = %.2f", phys.nmds$stress)
NMDS.stress

#Prepare for Plotting
phys.nmds.samples.plot.data <- cbind(phys.num, phys.nmds$points) 
phys.nmds.samples.plot.data$ID<-rownames(phys.nmds.samples.plot.data)
phys.nmds.samples.plot.data<-merge(phys.nmds.samples.plot.data, Phys.rm[,c(1:9)])
phys.nmds.samples.plot.data$Geno.Orig<-paste(phys.nmds.samples.plot.data$Genotype, phys.nmds.samples.plot.data$Orig, sep=".")
phys.nmds.samples.plot.data$Geno.Orig<-factor(phys.nmds.samples.plot.data$Geno.Orig, levels=c("AC8.N","AC8.T", "AC10.N", "AC10.T", "AC12.N", "AC12.T"), ordered=TRUE)
  
phys.nmds.covariates.plot.data <- data.frame(phys.envi.vectors)
phys.nmds.covariates.plot.data$Covariate<-c(rownames(phys.envi.vectors))
names(phys.nmds.covariates.plot.data)[1:2]<-c("MDS1", "MDS2")
phys.nmds.covariates.plot.data$Covariate
phys.nmds.covariates.plot.data$Metric<-c("Protein", "Biomass", "Chlorophyll")

```


#### Plot NMDS
```{r}
phys.nmds.plot<-ggplot(data = phys.nmds.samples.plot.data, aes(x = MDS1, y = MDS2)) + 
  geom_point(data = phys.nmds.samples.plot.data, aes(colour = Site.Orig, shape=Genotype), size = point.sz, alpha = 0.8) + 
  scale_colour_manual(values =Orig.colors.o)+
  scale_shape_manual(values=Geno.shapes.o)+
  theme_classic()+
  scale_x_continuous(limits = c(-4, 4))+
  scale_y_continuous(limits = c(-4, 4))+
   theme(axis.title.x = element_text(size = axis.title.sz), 
        axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"),
        legend.text=element_text(size=leg.txt.sz), 
        legend.title=element_text(size=leg.title.sz),
        legend.box.background = element_rect(color = "black"), 
         legend.position=c(0.5, .9), 
        legend.direction="horizontal", 
         legend.box = "horizontal")+
  geom_segment(aes(x = 0, y = 0, xend = MDS1, yend = MDS2), 
               data = phys.nmds.covariates.plot.data, 
               linewidth =bar.sz, alpha = 0.8, colour = "grey20",
               arrow=arrow(length=unit(0.01,"npc"))) +
  geom_text(data=phys.nmds.covariates.plot.data, 
            aes(x=MDS1, y=MDS2, label=Metric, fontface = "bold",
                hjust=0.15*(1-sign(MDS1)),vjust=0.6*(1-sign(MDS2))), 
            colour = "grey20", size=levels.sz-1)+
  labs(x="NMDS 1", y="NMDS 2", color="Site.Origin")+
  guides(color=guide_legend(nrow=2, byrow=TRUE), shape=guide_legend(nrow=2, byrow=TRUE)) +
  annotate("text", x = -4, y = -4, label = NMDS.stress, hjust=0, size=sig.sz); phys.nmds.plot
```



#### Data Prep for SIBER
```{r}
#Prepare data to be input into SIBER
phys.siber.input <- phys.nmds$points %>% #extract 2-D physiology axes from NMDS
  bind_cols(#add back covariates from original phys data
    select(Phys.log, ID:Site.Orig) #select those covariate columns
  ) %>% 
  #SIBER object needs four specific columns
  #iso1, iso2, group, and community
  #Need to manipulate the data to match that format
  #create  group variable that combines site, genotype and origin
  unite("group", c(Site, Genotype, Orig), sep = "_") %>% 
  #rename variables to match names for SIBER
  dplyr::rename(community = TimeP, iso1 = MDS1, iso2 = MDS2) %>% 
  #retain only columns for SIBER
  select(iso1, iso2, group, community)

```


### Run SIBER
Grouping by each genotype at each site, timepoint, and from each origin
```{r}
# create SIBER object
phys.siber.object <- createSiberObject(phys.siber.input)

# sample sizes by group
phys.siber.object[["sample.sizes"]]

# Calculate sumamry statistics for each group: TA, SEA and SEAc
group.ML <- groupMetricsML(phys.siber.object)
print(group.ML)

##Save SIBER Object in Outputs
save(phys.siber.object, file="Outputs/SIBER/Phys.SIBER.Object.RData")
```


# Ellipse Overlap

### Calculate Percent Overlap
Calculating the overlap of physiological niche space between corals of the same genotype, at the same site and timepoint. Using Percent of Niche Overlap as a more informative indicator of how similar the Transplant corals are to their Native clonemates. 

First, calculating the area of the 95% prediction ellipses for each group as well as the maximum likelihood of overlap between physiological niche ellipses.

#### Calculate Niche Areas
```{r}
##Load SIBER Object
load("Outputs/SIBER/Phys.SIBER.Object.RData")

##Dataframe of Niche Overlap
Niches<-data.frame(rbind(
  #KL AC10
  maxLikOverlap("TP1.KL_AC10_N", "TP1.KL_AC10_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP2.KL_AC10_N", "TP2.KL_AC10_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP3.KL_AC10_N", "TP3.KL_AC10_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP4.KL_AC10_N", "TP4.KL_AC10_T", phys.siber.object, p = 0.95, n = 360),
  #KL AC12
  maxLikOverlap("TP1.KL_AC12_N", "TP1.KL_AC12_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP2.KL_AC12_N", "TP2.KL_AC12_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP3.KL_AC12_N", "TP3.KL_AC12_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP4.KL_AC12_N", "TP4.KL_AC12_T", phys.siber.object, p = 0.95, n = 360),
  #KL AC8
  maxLikOverlap("TP1.KL_AC8_N", "TP1.KL_AC8_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP2.KL_AC8_N", "TP2.KL_AC8_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP3.KL_AC8_N", "TP3.KL_AC8_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP4.KL_AC8_N", "TP4.KL_AC8_T", phys.siber.object, p = 0.95, n = 360),
  #SS AC10
  maxLikOverlap("TP1.SS_AC10_N", "TP1.SS_AC10_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP2.SS_AC10_N", "TP2.SS_AC10_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP3.SS_AC10_N", "TP3.SS_AC10_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP4.SS_AC10_N", "TP4.SS_AC10_T", phys.siber.object, p = 0.95, n = 360),
  #SS AC12
  maxLikOverlap("TP1.SS_AC12_N", "TP1.SS_AC12_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP2.SS_AC12_N", "TP2.SS_AC12_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP3.SS_AC12_N", "TP3.SS_AC12_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP4.SS_AC12_N", "TP4.SS_AC12_T", phys.siber.object, p = 0.95, n = 360),
  #SS AC8
  maxLikOverlap("TP1.SS_AC8_N", "TP1.SS_AC8_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP2.SS_AC8_N", "TP2.SS_AC8_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP3.SS_AC8_N", "TP3.SS_AC8_T", phys.siber.object, p = 0.95, n = 360),
  maxLikOverlap("TP4.SS_AC8_N", "TP4.SS_AC8_T", phys.siber.object, p = 0.95, n = 360)
  ))

names(Niches)<-c("Native", "Transplant", "Overlap")
Niches$Site<-c(rep("KL", 12), rep("SS", 12))
Niches$Genotype<-c(rep(c(rep("AC10",4), rep("AC12",4), rep("AC8",4)),2))
Niches$TimeP<-rep(c("TP1", "TP2", "TP3", "TP4"),6)

```



#### Calculate Overlap
Proportion of niche overlap = Area of overlap between Native and Transplant physiological niches / Area of the Transplant niche

Proportion of the Transplant Niche area that is overlapping with that of the "expected" Native niche (if there was no lasting effect of Origin)
```{r}
##Proportion of Transplant Niche Area that is overlapping with Native niche
Niches$pOverlap<-Niches$Overlap/Niches$Transplant

##Write out dataframe
write.csv(Niches, "Outputs/SIBER/SIBER_p.Overlap.csv", row.names=FALSE)

##Read in
Niches<-read.csv("Outputs/SIBER/SIBER_p.Overlap.csv", header=TRUE)

##Set factor variables
Niches$Site<-factor(Niches$Site, levels=c("KL", "SS"))
Niches$Genotype<-factor(Niches$Genotype, levels=c("AC8", "AC10", "AC12"))
Niches$TimeP<-factor(Niches$TimeP)

##Add Months Post-Transplant
Niches$Time<-ifelse(Niches$TimeP=="TP1", 1, ifelse(Niches$TimeP=="TP2", 4, 
            ifelse(Niches$TimeP=="TP3", 8, ifelse(Niches$TimeP=="TP4", 12, NA))))
```


### Differences in Overlap

#### Across Time and Space
```{r}
##Check Normality
hist(Niches$pOverlap)
shapiro.test(Niches$pOverlap)
#Not Normal

hist(log(Niches$pOverlap+1))
shapiro.test(log(Niches$pOverlap+1))
#Normal

#Model with log(+1) transformation
##Model Overlap as a function of Site and Time, with Genotype as a Random effect
#Interaction between Site and Time
Overlap.lme<-lmer(log(pOverlap+1)~Time*Site+(1|Genotype), data=Niches)

##Check residuals
Overlap.lme_res <- simulateResiduals(fittedModel = Overlap.lme, plot = F)
plot(Overlap.lme_res)

##Model results
summary(Overlap.lme)

eta_squared(Overlap.lme)

##Save results
Overlap.res<-data.frame(summary(Overlap.lme)$coefficients[-1,])
Overlap.res$Predictor<-c("Time", "Site", "Time x Site")
Overlap.res$EtaSq<-c(eta_squared(Overlap.lme)$Eta2)
Overlap.res$Response<-rep("Niche Overlap", nrow(Overlap.res))
Overlap.res<-Overlap.res %>% dplyr::rename(SE = Std..Error, DF = df, t = t.value, p =Pr...t..)


```


### Plot Overlap over Time
```{r Plot Average Overlap over Time}
##Summary statistics by Site and Origin
pOverlap.orig.sum<-summarySE(Niches, measurevar="pOverlap", groupvars=c("Time"), na.rm=TRUE)

##Plot Average Niche Overlap across Timepoints
Overlap_Avg.plot<-ggplot(pOverlap.orig.sum, aes(x=Time, y=pOverlap, group=1)) + 
  scale_colour_manual(values=Site.colors.o)+
  geom_errorbar(aes(ymin=pOverlap-se, ymax=pOverlap+se), width=cap.sz, linewidth=bar.sz) +
  geom_path(linewidth=bar.sz)+
  geom_point(size=point.sz)+
  theme_classic()+
  theme(axis.title.x = element_text(size = axis.title.sz), 
        axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.position="none",
        strip.text=element_text(size=axis.title.sz))+
  scale_x_continuous(breaks = seq(0, 12, by = 2))+
  labs(x="Months Post-Transplant", y="Physiological Niche Overlap")+
  geom_text(aes(x=Time, y=pOverlap+se+0.03), label=c("TP1", "TP2", "TP3", "TP4"), size=sig.sz)+
  annotate("text", x=10, y=0.1, label="Time p < 0.05", size=sig.sz); Overlap_Avg.plot

```


# Environment and Physiology

### Merge Environmental and Physiology Data
```{r}
PhysEnv<-merge(Phys.rm, Env)

##Standardize Explanatory Environmental Variables for similar scales
PhysEnv.St <- PhysEnv %>% mutate_at(c("Temp.a", "Cond.a", "Temp.sd", "Cond.sd", "PAR.a", "Sat.hr", "TN_uM", "TP_uM", "N.P", "DIN_uM"), ~(scale(.) %>% as.vector))

##Set rownames
rownames(PhysEnv.St)<-PhysEnv.St$ID

```


### Distance Based RDA
```{r}
##dbRDA
PhysEnv.rda<-dbrda(vegdist(PhysEnv.St[, c(11:13)], "euclidean") ~ Temp.a + Cond.a + 
                         Temp.sd + Cond.sd + Sat.hr + N.P + DIN_uM, 
                       data=PhysEnv.St, dist="euclidean")

##Check model significance
anova(PhysEnv.rda)

##Check variance explained by model
(summary(PhysEnv.rda)$constr.chi/summary(PhysEnv.rda)$tot.chi)*100

RsquareAdj(PhysEnv.rda)$adj.r.squared*100

##Check variance explained by each Environmental metric
anova.cca(PhysEnv.rda, by="terms")

##Save Variance Results
PhysEnv.res<-data.frame(anova.cca(PhysEnv.rda, by="terms")[-8,])
PhysEnv.res$Predictor<-rownames(PhysEnv.res)
PhysEnv.res$Metric<-c("Average Temperature", "Average Conductivity", "Temperature Std. Dev.", 
                      "Conductivity Std. Dev.", "Hours Saturating Irradiance",
                      "Nitrogen:Phosphorus", "Dissolved Inorganic Nitrogen")
PhysEnv.res<-PhysEnv.res %>% dplyr::rename(DF = Df, p =Pr..F.)

##Quick dbRDA biplot
plot(PhysEnv.rda)
```


### RDA Plot

#### Prepare for Plotting
```{r}
##Samples
PhysEnv.rda.scores<- data.frame(scores(PhysEnv.rda)$sites)
PhysEnv.rda.scores$ID<- rownames(PhysEnv.rda.scores)
names(PhysEnv.St)
PhysEnv.rda.scores<-merge(PhysEnv.St[,c(1:9)], PhysEnv.rda.scores)

##Environmental Parameters
PhysEnv.rda.env<-data.frame(scores(PhysEnv.rda)$biplot[,1:2]) * ordiArrowMul(PhysEnv.rda, display="bp", fill=7)
PhysEnv.rda.env$Variable<-rownames(PhysEnv.rda.env)
PhysEnv.rda.env$Variable
PhysEnv.rda.env$Metric<-c("Temperature", "Conductivity", "Temp SD", "Cond SD", "Sat PAR", "N:P", "DIN")

##Adjusted R2
RDA.adjr2<- sprintf("%1.2f", (RsquareAdj(PhysEnv.rda)$adj.r.squared*100))
RDA.adjr2
```


#### Plot RDA
```{r}
phys.env.rda.plot<-ggplot(data = PhysEnv.rda.scores, aes(x = dbRDA1, y = dbRDA2)) + 
  geom_point(data = PhysEnv.rda.scores, size = point.sz, alpha = 0.8, aes(colour = Site.Orig, shape=Genotype)) + 
  theme_classic()+
  scale_colour_manual(values =Orig.colors.o)+
  scale_shape_manual(values=Geno.shapes.o)+
  theme(axis.title.x = element_text(size = axis.title.sz), 
        axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"),
        legend.text=element_text(size=leg.txt.sz), 
        legend.title=element_text(size=leg.title.sz),
        legend.box.background = element_rect(color = "black"), 
        legend.position=c(0.5, .9), 
        legend.direction="horizontal", 
         legend.box = "horizontal")+
  xlim(-18, 16)+
  ylim(-14, 20)+
  labs(x="RDA 1", y="RDA 2", color="Site.Origin")+
  geom_segment(aes(x = 0, y = 0, xend = dbRDA1, yend = dbRDA2), 
               data = PhysEnv.rda.env, linewidth =bar.sz, alpha = 0.5, colour = "grey30", 
               arrow=arrow(length=unit(0.01,"npc")))+
  geom_text(data = PhysEnv.rda.env[-c(3,6),], 
            aes(x=c(dbRDA1[1:2], dbRDA1[3]-3, dbRDA1[4:5]), 
                y = c(dbRDA2[1:2], dbRDA2[3]+1, dbRDA2[4:5]), 
                label = Metric, fontface = "bold",
                hjust=0.5*(1-sign(dbRDA1)),vjust=0.6*(1-sign(dbRDA2))), 
            colour = "grey20", size=levels.sz-1)+
  geom_text(data = PhysEnv.rda.env[c(3,6),], 
            aes(x=c(dbRDA1[1], dbRDA1[2]+2), y = c(dbRDA2[1]+2, dbRDA2[2]-1), 
                label = Metric, hjust=0.5*(1-sign(dbRDA1)),vjust=0.6*(1-sign(dbRDA2))), 
            colour = "grey20", size=levels.sz-1)+
  guides(color=guide_legend(nrow=2, byrow=TRUE), shape=guide_legend(nrow=2, byrow=TRUE)) +
  annotate("text", x = -18, y = -14, label = expression(paste('Adjusted R'^2*"= 27.10%")), hjust=0, size=sig.sz); phys.env.rda.plot
  
```


# Figures

### Figure 2
```{r}
##Create Panel
PhysEnvNiche_fig<-plot_grid(phys.env.rda.plot, phys.nmds.plot, Overlap_Avg.plot,
                    rel_widths=1, 
                    rel_heights = 1,
                    nrow=1, ncol=3, byrow=T, 
                    labels = c("A", "B", "C"), 
                    label_size=panel.lab.sz)

##Save Figure
ggsave(filename="Figures/02_Physiology/Fig2_PhysEnvNiches.png", plot=PhysEnvNiche_fig, dpi=300, width=20, height=6.5, units="in")

```


# Tables

### Table 1
```{r}
##Save Results
Table1_EnvPredictors<-PhysEnv.res

##Organize
Table1_EnvPredictors<-Table1_EnvPredictors[,c("Metric", "DF", "Variance", "F", "p")]

##Add Significance levels
Table1_EnvPredictors$Significance<-ifelse(Table1_EnvPredictors$p<0.001, "***",
                              ifelse(Table1_EnvPredictors$p<0.01, "**",
                                     ifelse(Table1_EnvPredictors$p<0.05, "*",
                                            ifelse(Table1_EnvPredictors$p<0.1, "-", ""))))
#Round to 3 digits
Table1_EnvPredictors$Variance<-round(Table1_EnvPredictors$Variance, 3)
Table1_EnvPredictors$F<-round(Table1_EnvPredictors$F, 3)

##Rename Predictor
Table1_EnvPredictors<-Table1_EnvPredictors %>% dplyr::rename(Predictor = Metric)

##Write Out Table
write.csv(Table1_EnvPredictors, "Tables/Table1_Physiology_Environment_RDA_Results.csv", row.names=FALSE)

```


### Table S3B Niche Overlap
```{r}
##Save Results
TableS3B_NicheOverlap<-Overlap.res

##Organize
TableS3B_NicheOverlap<-TableS3B_NicheOverlap[,c("Response", "Predictor", "Estimate", "SE", "DF", "t", "p", "EtaSq")]

#Round to 3 digits
TableS3B_NicheOverlap$Estimate<-round(TableS3B_NicheOverlap$Estimate, 3)
TableS3B_NicheOverlap$SE<-round(TableS3B_NicheOverlap$SE, 3)
TableS3B_NicheOverlap$t<-round(TableS3B_NicheOverlap$t, 3)
TableS3B_NicheOverlap$p<-round(TableS3B_NicheOverlap$p, 3)
TableS3B_NicheOverlap$EtaSq<-round(TableS3B_NicheOverlap$EtaSq, 3)

#Integer
TableS3B_NicheOverlap$DF<-round(TableS3B_NicheOverlap$DF, 0)

##Write Out Table
write.csv(TableS3B_NicheOverlap, "Tables/TableS3B_Niche_Overlap_Results.csv", row.names=FALSE)

```