library(vegan)
library(dave)
library(ggplot2)
library(lme4)
library(gridExtra)
library(devtools)
library(multcomp)
require(ggbiplot)

data1 <- read.table("r_vegdatanew.csv", header = TRUE, sep = ",")

veg.pca <- princomp(data1[,3:19], cor = TRUE)
print(veg.pca)
summary(veg.pca, loadings=TRUE)

biplot = ggbiplot(pcobj = veg.pca, choices = c(1,2), obs.scale = 1, 
  var.scale = 1, ellipse = TRUE, groups = data1$Habitat) +
  labs(colour = "Habitat") +
  theme_classic() +
  theme(axis.text = element_text(size = 11))  # Increase font size for axes

print(biplot)
# To print or save the combined plot
ggsave("Figure4.pdf", biplot, width = 3000, units = "px")

cor(data1[,3:19])

scores <- veg.pca$scores
loadings <- veg.pca$loadings

library(writexl)
write.csv(scores, "scores.csv")
write.csv(loadings, "loadings.csv")

lr.data <- read.csv("pc_outnew.csv")

plot1TT <- ggplot(lr.data, aes(x = PC1, y = tthet_cor)) +
  theme_classic() +
  geom_point(shape = 1) +
  geom_smooth(method = lm) +
  labs(y = "T. thetis Detections")
plot2TT <- ggplot(lr.data, aes(x = PC2, y = tthet_cor)) + 
  theme_classic() +
  geom_point(shape = 1) +
  geom_smooth(method = lm) +
  labs(y = "T. thetis Detections")



# Linear regression graph #
grid.arrange(plot1TT, plot2TT, ncol = 1, nrow = 2)

# Linear regression code #
tthet.pc <- glm(tthet_cor ~ PC1 + PC2 + PC3, data=lr.data)
summary(tthet.pc)

plot1TS <- ggplot(lr.data, aes(x = PC1, y = thstig_cor)) +
  theme_classic() +
  geom_point(shape = 1) +
  geom_smooth(method = glm) +
  labs(y = "T. stigmatica Detections")
plot2TS <- ggplot(lr.data, aes(x = PC2, y = thstig_cor)) + 
  theme_classic() +
  geom_point(shape = 1) +
  geom_smooth(method = glm) +
  labs(y = "T. stigmatica Detections")

# Linear regression graph #
grid.arrange(plot1TS, plot2TS, ncol = 1, nrow = 2)

tstig.pc <- lm(thstig_cor ~ PC1 + PC2 + PC3, data=lr.data)
summary(tstig.pc)

################# ANOVA

pade_anova <- read.table("eventsbyforest_animalsperday2.csv", header = TRUE, sep = ",")

pade_aov <- aov(logevents_corr ~ forest * species, data=pade_anova)
summary(pade_aov)

# Tukey's HSD test for post hoc analysis
posthoc_results <- TukeyHSD(pade_aov, "forest:species")
print(posthoc_results)

# Tukey post-hoc test for the interaction between forest and species
tukey_interaction <- glht(pade_aov, linfct = mcp(forest:species = "Tukey"))
summary(tukey_interaction)

residuals <- residuals(pade_aov)
hist(residuals, main = "Histogram of Residuals", xlab = "Residuals")

qqnorm(residuals)
qqline(residuals, col = "red")

shapiro.test(residuals)

ws_anova <- subset(pade_anova, forest == "WS")

#### charts

sp_means <- aggregate(pade_anova[,5], by=list(pade_anova$species), mean,na.rm=TRUE) # get means
ws_means <- aggregate(ws_anova[,5], by=list(ws_anova$species), mean,na.rm=TRUE) # get means
sp_hab_means <- aggregate(logevents_corr ~ species + forest, data = pade_anova, mean, na.rm = TRUE)

boxplot(logevents ~ species, xlab = "Species", ylab = "Independent events/camera", main = "(a)",  
        par(las=1), frame.plot = FALSE, whiskcol = "black", whisklty = "solid", boxwex = 0.7, 
        col = "grey90", outline = FALSE, varwidth = FALSE, data = pade_anova)

# Add points with stripchart
stripchart(logevents ~ species, data = pade_anova, 
           method = "jitter", # Adds a small amount of random variation to the position of the points
           pch = 1,         # Type of point. 20 is a solid circle
           col = "black",    # Color of points
           add = TRUE,       # Adds the stripchart to the existing plot
           vertical = TRUE)  # Aligns points in the same direction as the boxplot
points(sp_means, col="black", lwd=1, pch = 21, bg = "black", cex = 2)


boxplot(logevents ~ species, xlab = "Species", ylab = "Independent events/camera", main = "(b)",  
        par(las=1), frame.plot = FALSE, whiskcol = "black", whisklty = "solid", boxwex = 0.7, 
        col = "grey90", outline = FALSE, varwidth = FALSE, data = ws_anova)

stripchart(logevents ~ species, data = ws_anova, 
           method = "jitter", # Adds a small amount of random variation to the position of the points
           pch = 1,         # Type of point. 20 is a solid circle
           col = "black",    # Color of points
           add = TRUE,       # Adds the stripchart to the existing plot
           vertical = TRUE)  # Aligns points in the same direction as the boxplot
points(ws_means, col="black", lwd=1, pch = 21, bg = "black", cex = 2)


boxplot(logevents ~ species*forest, xlab = "Species", ylab = "Independent events/camera", main = "",  
        par(las=1), frame.plot = FALSE, whiskcol = "black", whisklty = "solid", boxwex = 0.7, 
        col = "grey90", outline = FALSE, varwidth = FALSE, data = pade_anova)

stripchart(logevents ~ species*forest, data = pade_anova, 
           method = "jitter", # Adds a small amount of random variation to the position of the points
           pch = 1,         # Type of point. 20 is a solid circle
           col = "black",    # Color of points
           add = TRUE,       # Adds the stripchart to the existing plot
           vertical = TRUE)  # Aligns points in the same direction as the boxplot
with(sp_hab_means, points(x = as.numeric(interaction(species, forest)), 
                          y = logevents, lwd=1, pch = 21, bg = "black", cex = 2))



##### test

# Create the boxplot without x-axis labels
boxplot(logevents_corr ~ species*forest, 
        xlab = "Species by Habitat", 
        ylab = "Independent events/camera", 
        main = "",  
        par(las=1), 
        frame.plot = FALSE, 
        whiskcol = "black", 
        whisklty = "solid", 
        boxwex = 0.7, 
        col = "grey90", 
        outline = FALSE, 
        varwidth = FALSE, 
        axes = FALSE, # Prevents drawing of axes
        data = pade_anova)

# Add stripchart
stripchart(logevents_corr ~ species*forest, 
           data = pade_anova, 
           method = "jitter",
           pch = 1, 
           col = "black",
           add = TRUE, 
           vertical = TRUE)

# Add points
with(sp_hab_means, points(x = as.numeric(interaction(species, forest)), 
                          y = logevents_corr, 
                          lwd = 1, 
                          pch = 21, 
                          bg = "black", 
                          cex = 2))

# Manually add x-axis labels
labels <- c(expression(italic("T. stigmatica") ~ "Rainforest"), expression(italic("T. thetis") ~ "Rainforest"), expression(italic("T. stigmatica") ~ "Sclerophyll"), expression(italic("T. thetis") ~ "Sclerophyll"))  # Replace with your actual labels
axis(1, at = 1:length(labels), labels = labels, las = 1)  # Adjusts the label angles with 'las'

# Add y-axis
axis(2)



########## end of anova

thetis <- subset(pade_anova, species == "TTHET")
thetis_aov <- ?aov(logcount ~ habitat, data=thetis)
summary(thetis_aov)

####tests for normality:
residuals2 <- residuals(thetis_aov)
hist(residuals2, main = "Histogram of Residuals", xlab = "Residuals")

shapiro.test(residuals2)
boxplot(logcount ~ habitat, data = thetis, main = "Boxplot of Count by Habitat")

stigmatica <- subset(pade_anova, species == "TSTIG")
stigmatica_aov <- aov(logcount ~ habitat, data=stigmatica)
summary(stigmatica_aov)

####tests for normality:
residuals2 <- residuals(stigmatica_aov)
hist(residuals2, main = "Histogram of Residuals", xlab = "Residuals")

qqnorm(residuals2)
qqline(residuals2, col = "red")

shapiro.test(residuals2)
boxplot(count ~ habitat, data = thetis, main = "Boxplot of Count by Habitat")



####

boxplot(TTHET ~ Type, xlab = "Forest Type", ylab = "Independent events/camera", main = "(a)",  par(las=1), frame.plot = FALSE, 
        whiskcol = "black", whisklty = "solid", boxwex = 0.7, col = "grey90", outline = FALSE, 
        varwidth = FALSE, data=lr.data)
points(TTHET_means, col="black", lwd=1, pch = 21, bg = "black", cex = 2)


boxplot(TSTIG ~ Type, xlab = "Forest Type", ylab = "Independent events/camera", main = "(b)", par(las=1), frame.plot = FALSE, 
        whiskcol = "black", whisklty = "solid", boxwex = 0.7, col = "grey90", outline = FALSE, 
        varwidth = FALSE, data=lr.data)
points(TSTIG_means, col="black", lwd=1, pch = 21, bg = "black", cex = 2)


### individual covariates - TTHET:

# Make sure necessary packages are loaded
if (!requireNamespace("patchwork", quietly = TRUE)) install.packages("patchwork")
library(patchwork)
library(ggplot2)

# Assuming 'lr.data' and 'data1' are correctly set up, and variables_to_plot is defined.

# Initialize an empty list to store plots
plot_list <- list()

# Loop through each variable, calculate correlation, and generate a plot
for (var in variables_to_plot) {
  # Create a temporary dataframe
  temp_df <- data.frame(tthet_cor = lr.data$tthet_cor, variable_value = data1[[var]])
  
  # Calculate correlation and p-value
  correlation_test <- cor.test(temp_df$variable_value, temp_df$tthet_cor, method = "pearson")
  p_value <- correlation_test$p.value
  
  # Generate the plot
  p <- ggplot(temp_df, aes(x = variable_value, y = tthet_cor)) +
    theme_classic() +
    geom_point(shape = 1) +
    geom_smooth(method = "lm", se = FALSE) +
    labs(x = var, y = "T. thetis Detections", title = paste("T. thetis Detections vs", var, "\nP-value:", format(p_value, digits = 3))) +
    theme(plot.title = element_text(size = 10)) # Adjust title size as needed
  
  # Alternatively, add p-value as an annotation
  # p <- p + annotate("text", x = Inf, y = Inf, label = paste("P-value:", format(p_value, digits = 3)), 
  #                   hjust = 1.1, vjust = 1.1, size = 3)
  
  # Add the plot to the list
  plot_list[[var]] <- p
}

# Combine all plots using patchwork
combined_plot <- wrap_plots(plot_list, ncol = 3) # Adjust 'ncol' as needed

# Print the combined plot
print(combined_plot)

### individual covariates - TSTIG:

# Make sure necessary packages are loaded
library(patchwork)
library(ggplot2)

# Assuming 'lr.data' and 'data1' are correctly set up, and variables_to_plot is defined.

# Initialize an empty list to store plots
plot_list <- list()

# Loop through each variable, calculate correlation, and generate a plot
for (var in variables_to_plot) {
  # Create a temporary dataframe
  temp_df2 <- data.frame(thstig_cor = lr.data$thstig_cor, variable_value = data1[[var]])
  
  # Calculate correlation and p-value
  correlation_test <- cor.test(temp_df2$variable_value, temp_df2$thstig_cor, method = "spearman")
  p_value <- correlation_test$p.value
  
  # Assuming variables_to_plot is your list of variables
  num_tests <- length(variables_to_plot)
  alpha <- 0.05 / num_tests  # Bonferroni correction
  
  # Then, when you're evaluating the p-values:
  # Suppose p_value is the p-value obtained from cor.test()
  is_significant <- p_value < alpha  # This checks if the p-value is significant after correction
  
  # Generate the plot
  p <- ggplot(temp_df2, aes(x = variable_value, y = thstig_cor)) +
    theme_classic() +
    geom_point(shape = 1) +
    geom_smooth(method = "lm", se = FALSE) +
    labs(x = var, y = "T. stigmatica Detections", title = paste("T. stigmatica Detections vs", var, "\nP-value:", format(p_value, digits = 3))) +
    theme(plot.title = element_text(size = 10)) # Adjust title size as needed
  
  # Alternatively, add p-value as an annotation
  # p <- p + annotate("text", x = Inf, y = Inf, label = paste("P-value:", format(p_value, digits = 3)), 
  #                   hjust = 1.1, vjust = 1.1, size = 3)
  
  # Add the plot to the list
  plot_list[[var]] <- p
}

# Combine all plots using patchwork
combined_plot <- wrap_plots(plot_list, ncol = 3) # Adjust 'ncol' as needed

# Print the combined plot
print(combined_plot)