Golshan is glimmering.Rmd

---
title: "lol final project"
output: pdf_document
date: "2024-11-14"
---

```{r, message=F, warning=F}
require(tidyverse)
require(deSolve)
require(vegan)
require(car)
require(rgl)
require(ggfortify)

options(rgl.useNULL = TRUE) 
rgl::setupKnitr(autoprint = TRUE)
```

```{r}
install.packages("ggfortify")
```

```{r}
install.packages("rgl")
```

```{r}
install.packages("fitdistrplus")
require(fitdistrplus)
```

## R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see <http://rmarkdown.rstudio.com>.

When you click the **Knit** button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

```{r}
read_csv("Arthropod_List.csv") -> Arthropod_List

read_csv("Plant_Traits.csv") -> Plant_Traits

read_csv("Soil_Traits.csv") -> Soil_Traits

read_csv("Stable_Isotopes.csv") -> Stable_Isotopes

read_csv("final_arthro.csv") -> Arthropod_List_Updated

```

```{r}
Arthropod_list[1,]
```


#Shannon Diversity Index
```{r}
# Sample data: abundance of species
abundance_comm1 <- c()
abundance_comm2 <- c()
abundance_comm3 <- c()
abundance_comm4 <- c()
abundance_comm5 <- c()

shannon_comm1 <- NULL
shannon_comm2 <- NULL
shannon_comm3 <- NULL
shannon_comm4 <- NULL
shannon_comm5 <- NULL

# For loop, First separate by community, and then sum the vector of the row of that species and put into an abundance vector

for (community in 1:5){
  if (community == 1){ #in OP
    for (species in 1:349){
      abundance_comm1[species] = sum(Arthropod_List[species,5:19])
    }
    shannon_comm1 <- diversity(abundance_comm1, index = "shannon")
  }
  if (community == 2){ #in OP
    for (species in 1:349){
      abundance_comm2[species] = sum(Arthropod_List[species,20:34])
    }
    shannon_comm2 <- diversity(abundance_comm2, index = "shannon")
  }
  if (community == 3){ #in OP
    for (species in 1:349){
      abundance_comm3[species] = sum(Arthropod_List[species,35:49])
    }
    shannon_comm3 <- diversity(abundance_comm3, index = "shannon")
  }
  if (community == 4){ #in OP
    for (species in 1:349){
      abundance_comm4[species] = sum(Arthropod_List[species,50:64])
    }
    shannon_comm4 <- diversity(abundance_comm4, index = "shannon")
  }
  if (community == 5){ #in OP
    for (species in 1:349){
      abundance_comm5[species] = sum(Arthropod_List[species,65:79])
    }
    shannon_comm5 <- diversity(abundance_comm5, index = "shannon")
  }
}
```


```{r}
shannon_frame = data.frame(Arthropod_List_Updated$Plot, "Shannon_Diversity")

for (i in 1:75){
  abundance_vec = Arthropod_List_Updated[i,2:250]
  diversity = diversity(abundance_vec, index = "shannon")
  shannon_frame[i,2] = diversity
}
```

```{r}
fixed_arthro = Arthropod_List_Updated[,2:350]

shannon_indices = diversity(fixed_arthro, index = "shannon")
view(shannon_indices)
```

```{r}
final_shannon = data.frame(Treatment = Arthropod_List_Updated$Plot,
                           Shannon = shannon_indices)
```

```{r}
ggplot(final_shannon, aes(x = Treatment, y = Shannon)) +
  geom_bar(stat = "identity", fill = "skyblue", color = "black") +
  theme_minimal() +
  labs(title = "Shannon Diversity for each Treatment", x = "Treatment", y = "Shannon Diversity")
```
```{r}
install.packages("lme4")    # For mixed effects models
install.packages("lmerTest")
```


```{r}
library(lme4)
library(lmerTest)
```


```{r}
soil_C = Soil_Traits$Soil_C
soil_N = Soil_Traits$Soil_N
soil_P = Soil_Traits$Soil_P
soil_pH = Soil_Traits$Soil_pH
soil_sal = Soil_Traits$Soil_Salinity
soil_wat = Soil_Traits$Soil_Water

biomass_above = Plant_Traits$Aboveground_Biomass
plant_dens = Plant_Traits$Plant_Density
leaf_N = Plant_Traits$Leaf_N
leaf_C = Plant_Traits$Leaf_C
leaf_P = Plant_Traits$Leaf_P
```


```{r}
final_shannon$soil_C <- soil_C
final_shannon$soil_N <- soil_N
final_shannon$soil_P <- soil_P
final_shannon$soil_pH <- soil_pH
final_shannon$soil_sal <- soil_sal
final_shannon$soil_wat <- soil_wat

final_shannon$biomass_above <- biomass_above
final_shannon$plant_dens <- plant_dens
final_shannon$leaf_N <- leaf_N
final_shannon$leaf_C <- leaf_C
final_shannon$leaf_P <- leaf_P

```

```{r}
write.csv(final_shannon, "final_shannon.csv")
```


```{r}
# Example: Group treatments into broader categories
final_shannon$Treatment_grouped <- ifelse(startsWith(final_shannon$Treatment, "OP"), "OP",
                                          ifelse(startsWith(final_shannon$Treatment, "RP"), "RP",
                                                 ifelse(startsWith(final_shannon$Treatment, "TP"), "TP",
                                                        ifelse(startsWith(final_shannon$Treatment, "IS"), "IS",
                                                               ifelse(startsWith(final_shannon$Treatment, "PS"), "PS", "Other")))))

# Fit the model with grouped Treatment
model2 <- lmer(Shannon ~ Treatment_grouped, data = final_shannon, family = Gamma(link = "log"))
summary(model2)
```


```{r}
model <- glmer(Shannon ~ soil_sal + biomass_above + plant_dens + leaf_N + (1 | Treatment), data = final_shannon)
summary(model)

run a PCA on all the variables - 
```

```{r}
final_shannon$soil_C <- as.numeric(as.character(final_shannon$soil_C))
final_shannon$soil_N <- as.numeric(as.character(final_shannon$soil_N))
final_shannon$soil_P <- as.numeric(as.character(final_shannon$soil_P))
final_shannon$soil_sal <- as.numeric(as.character(final_shannon$soil_sal))
final_shannon$soil_wat <- as.numeric(as.character(final_shannon$soil_wat))
final_shannon$biomass_above <- as.numeric(as.character(final_shannon$biomass_above))
final_shannon$plant_dens <- as.numeric(as.character(final_shannon$plant_dens))
final_shannon$leaf_N <- as.numeric(as.character(final_shannon$leaf_N))
final_shannon$leaf_C <- as.numeric(as.character(final_shannon$leaf_C))
final_shannon$leaf_P <- as.numeric(as.character(final_shannon$leaf_P))
```

```{r}
is.na(final_shannon$soil_C)
is.na(final_shannon$soil_N)
is.na(final_shannon$soil_P)
is.na(final_shannon$soil_sal)
is.na(final_shannon$soil_wat)
is.na(final_shannon$biomass_above)
is.na(final_shannon$plant_dens)
is.na(final_shannon$leaf_N)
is.na(final_shannon$leaf_C)
is.na(final_shannon$leaf_P)
```



```{r}
noNA_final_shannon <- filter(final_shannon, !is.na(final_shannon$leaf_N))
```


```{r}
mean(noNA_final_shannon$leaf_N)
```

```{r}
read.csv("final_shannon_noNA.csv") -> final_shannon_UPDATED
```


```{r}
par(mfrow=c(1,1))

# can look at all variables together
# can see that all morphological traits are highly linearly and positively correlated
# ie species with large masses ALSO have high wing length and high tarsus length (lower leg bone, the bone below what looks like a backwards facing knee)
pairs(final_shannon_UPDATED[, c("soil_C", "soil_N", "soil_P", "soil_pH", 
                                    "soil_sal", "soil_wat", "biomass_above", 
                                    "plant_dens", "leaf_N", "leaf_C", "leaf_P")])
```

```{r}
pca.out <- prcomp(final_shannon_UPDATED[, c("soil_C", "soil_N", "soil_P", "soil_pH", 
                                    "soil_sal", "soil_wat", "biomass_above", 
                                    "plant_dens", "leaf_N", "leaf_C", "leaf_P")], scale.=T)
```
```{r}
pca.out
```


```{r}
summary(pca.out)
```


```{r}
head(pca.out$rotation)
```
```{r}
#soil_n, water, ph, c, aboveground biomass
#
```

```{r}

```


```{r}
autoplot(pca.out, x=1,y=2,
         data=final_shannon_UPDATED, 
         colour="Treatment_grouped",
         frame=TRUE) +
  theme_classic()
```


```{r}
cor_matrix <- cor(final_shannon[, c("soil_C", "soil_N", "soil_P", "soil_pH", 
                                    "soil_sal", "soil_wat", "biomass_above", 
                                    "plant_dens", "leaf_N", "leaf_C", "leaf_P")], 
                  use = "pairwise.complete.obs")
print(cor_matrix)
```