#Analysis of Thule seabird data (2010-2023) ----
#Kimberlee Whitmore (whitmorekimberlee@gmail.com) Texas Christian University 2024
setwd("C:/LearningR/R code")
library(tidyr)   # Formatting data for analysis
library (dplyr)  # Manipulating data
library (ggplot2)#Visualizing results
library (readr)  #Manipulating data

#Plot formatting----
theme_Publication <- function(base_size=14, base_family="helvetica") {
  library(grid)
  library(ggthemes)
  (theme_foundation(base_size=base_size, base_family=base_family)
    + theme(plot.title = element_text(face = "bold",
                                      size = rel(1.2), hjust = 0.5),
            text = element_text(),
            panel.background = element_rect(colour = NA),
            plot.background = element_rect(colour = NA),
            panel.border = element_rect(colour = NA),
            axis.title = element_text(face = "bold",size = rel(1)),
            axis.title.y = element_text(angle=90,vjust =2),
            axis.title.x = element_text(vjust = -0.2),
            axis.text = element_text(), 
            axis.line = element_line(colour="black"),
            axis.ticks = element_line(),
            panel.grid.major = element_line(colour="#f0f0f0"),
            panel.grid.minor = element_blank(),
            legend.key = element_rect(colour = NA),
            legend.position = "bottom",
            legend.direction = "horizontal",
            legend.key.size= unit(0.2, "cm"),
            legend.margin = unit(0, "cm"),
            legend.title = element_text(face="italic"),
            plot.margin=unit(c(10,5,5,5),"mm"),
            strip.background=element_rect(colour="#f0f0f0",fill="#f0f0f0"),
            strip.text = element_text(face="bold")
    ))
  
}
thule_birds <- read.csv("C:/LearningR/R code/thule_birds_r.csv")
#filtering out 2015
thule_birds <- subset(thule_birds, year !=2015)

#Data prep ----

str(thule_birds) #columns with blanks need to be reclassed as numeric

thule_birds$Weight <- as.numeric(thule_birds$Weight)
thule_birds$Beak <- as.numeric(thule_birds$Beak)
thule_birds$Tarsus <- as.numeric(thule_birds$Tarsus)
thule_birds$Tail <- as.numeric(thule_birds$Tail)
thule_birds$Wing <- as.numeric(thule_birds$Wing)
thule_birds$Head <- as.numeric(thule_birds$Head)
str(thule_birds)
thule_birds$year1 <- thule_birds$year - min(thule_birds$year) + 1
thule_birds$lnHg <- log(thule_birds$Hg) #natural log
sp_order <- c("DOVE", "BLKI", "ATPU", "BLGU", "TBMU") #change order of graphs
thule_birds$sp_code <- factor(thule_birds$sp_code, levels = sp_order)
thule_birds2 <- thule_birds[!is.na(thule_birds$carbon), ] #only samples with isotope data


#Data frames for each species ----

ATPU <- filter(thule_birds, sp_code == "ATPU")
BLGU <- filter(thule_birds, sp_code == "BLGU")
BLKI <- filter(thule_birds, sp_code == "BLKI")
DOVE <- filter(thule_birds, sp_code == "DOVE")
TBMU <- filter(thule_birds, sp_code == "TBMU")

#Data summary ------
data_summary <- thule_birds %>%
  group_by(sp_code) %>%
  summarise(
    Hg_Average = mean(Hg),
    Hg_StandardDeviation = sd(Hg),
    StandardError = sd(Hg) / sqrt(n()),
    Minimum = min(Hg),
    Maximum = max(Hg),
    N = n()
  ) %>%
  arrange(sp_code)

SI_av_data <- thule_birds2 %>%
  group_by(sp_code) %>%
  summarise(
    TotalCount = n(),        # Count of occurrences
    MeanlnHg = mean(lnHg), 
    SDlnHg = sd(lnHg),
    MeanC = mean(carbon),
    SDC = sd(carbon),
    MeanN = mean(nitrogen),
    SDN = sd(nitrogen)
  )
#Data analysis- Hg vs nitrogen----

#linear regression for ATPU Hg vs nitrogen
ATPUlinreg = lm(formula = lnHg ~ nitrogen,
                data = ATPU)
summary(ATPUlinreg)
#create scatterplot

ggplot(ATPU, aes(x = nitrogen, y = logHg, colour = as.factor(year))) +
  geom_point() +
  geom_line(aes(x = nitrogen, y = predict(ATPUlinreg, newdata = ATPU)), colour = 'blue') +
  ggtitle('ATPU logHg vs N') +
  xlab('Nitrogen') +
  ylab('logHg') +
  theme_classic()
 #separate linear regressions (no formula)
ggplot(ATPU, aes(x = nitrogen, y = lnHg, colour = as.factor(year))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  ggtitle('ATPU logHg vs N') +
  xlab('Nitrogen') +
  ylab('logHg') +
  theme_classic()

#linear regression for BLGU Hg vs nitrogen
BLGUlinreg = lm(formula = lnHg ~ nitrogen,
                data = BLGU)
summary(BLGUlinreg)
#create scatterplot
ggplot(BLGU, aes(x = nitrogen, y = logHg, colour = as.factor(year))) +
  geom_point() +
  geom_line(aes(x = nitrogen, y = predict(BLGUlinreg, newdata = BLGU)), colour = 'blue') +
  ggtitle('BLGU logHg vs N') +
  xlab('Nitrogen') +
  ylab('logHg') +
  theme_classic()
#separate linear regressions (no formula)
ggplot(BLGU, aes(x = nitrogen, y = logHg, colour = as.factor(year))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  ggtitle('BLGU logHg vs N') +
  xlab('Nitrogen') +
  ylab('logHg') +
  theme_classic()

#linear regression for BLKI Hg vs nitrogen
BLKIlinreg = lm(formula = logHg ~ nitrogen,
         data = BLKI)
summary(BLKIlinreg)
#create scatterplot
ggplot(BLKI, aes(x = nitrogen, y = logHg, colour = as.factor(year))) +
  geom_point() +
  geom_line(aes(x = nitrogen, y = predict(BLKIlinreg, newdata = BLKI)), colour = 'blue') +
  ggtitle('BLKI logHg vs N') +
  xlab('Nitrogen') +
  ylab('logHg') +
  theme_classic()
#separate linear regressions (no formula)
ggplot(BLKI, aes(x = nitrogen, y = logHg, colour = as.factor(year))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  ggtitle('BLKI logHg vs N') +
  xlab('Nitrogen') +
  ylab('logHg') +
  theme_classic()

#linear regression for DOVE Hg vs nitrogen
DOVElinreg = lm(formula = logHg ~ nitrogen,
                data = DOVE)
summary(DOVElinreg)
#create scatterplot
ggplot(DOVE, aes(x = nitrogen, y = logHg, colour = as.factor(year))) +
  geom_point() +
  geom_line(aes(x = nitrogen, y = predict(DOVElinreg, newdata = DOVE)), colour = 'blue') +
  ggtitle('DOVE logHg vs N') +
  xlab('Nitrogen') +
  ylab('logHg') +
  theme_classic()
#separate linear regressions (no formula)
ggplot(DOVE, aes(x = nitrogen, y = logHg, colour = as.factor(year))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  ggtitle('DOVE logHg vs N') +
  xlab('Nitrogen') +
  ylab('logHg') +
  theme_classic()

#linear regression for TBMU Hg vs nitrogen
TBMUlinreg = lm(formula = logHg ~ nitrogen,
                data = TBMU)
summary(TBMUlinreg)
#create scatterplot
ggplot(TBMU, aes(x = nitrogen, y = logHg, colour = as.factor(year))) +
  geom_point() +
  geom_line(aes(x = nitrogen, y = predict(TBMUlinreg, newdata = TBMU)), colour = 'blue') +
  ggtitle('TBMU logHg vs N') +
  xlab('Nitrogen') +
  ylab('logHg') +
  theme_classic()
#separate linear regressions (no formula)
ggplot(TBMU, aes(x = nitrogen, y = logHg, colour = as.factor(year))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  ggtitle('TBMU logHg vs N') +
  xlab('Nitrogen') +
  ylab('logHg') +
  theme_classic()


#HgN Analysis across species
ALLlinreg = lm(formula = lnHg ~ nitrogen,
                data = thule_birds)
summary(ALLlinreg)

ggplot(thule_birds, aes(x = nitrogen, y = lnHg, colour = as.factor(sp_number))) +
  geom_point() +
  geom_smooth(method = 'lm', aes(group = 1), se = FALSE, color = "black") +
  ggtitle('All species logHg vs N') +
  xlab('Nitrogen') +
  ylab('logHg') +
  theme_classic() 
#species averages
ALLlinreg = lm(formula = MeanlnHg ~ MeanN,
               data = SI_av_data)
summary(ALLlinreg)

library(ggplot2)

# Your existing ggplot code
ggplot(SI_av_data, aes(x = MeanN, y = MeanlnHg, colour = as.factor(sp_code))) +
  geom_abline(intercept = 1.59, slope = 0.315, color = "black", size = 1) +  # Adding regression of all points  
  geom_point() +
  geom_errorbar(aes(x = MeanN, ymin = MeanlnHg - SDlnHg, ymax = MeanlnHg + SDlnHg), width = 0.01) +
  geom_errorbarh(aes(y = MeanlnHg, xmin = MeanN - SDN, xmax = MeanN + SDN), height = 0.1) +
  ggtitle('All species lnHg vs N') +
  xlab('δ15N') +
  ylab('lnHg') +
  theme_classic()

#HgC Analysis accross species
ALLlinregC = lm(formula = logHg ~ carbon,
               data = thule_birds)
summary(ALLlinregC)

ggplot(thule_birds, aes(x = carbon, y = logHg, colour = as.factor(sp_number))) +
  geom_point() +
  geom_smooth(method = 'lm', aes(group = 1), se = FALSE, color = "black") +
  ggtitle('All species lnHg vs C') +
  xlab('Carbon') +
  ylab('logHg') +
  theme_classic() 
#species averages
ALLlinregC = lm(formula = MeanlnHg ~ MeanC,
               data = SI_av_data)
summary(ALLlinregC)

ggplot(SI_av_data, aes(x = MeanC, y = MeanlnHg, colour = as.factor(sp_code))) +
  geom_point() +
  geom_smooth(method = 'lm', aes(group = 1), se = FALSE, color = "black") +
  geom_errorbar(aes(x = MeanC, ymin = MeanlnHg - SDlnHg, ymax = MeanlnHg + SDlnHg), width = 0.01) +
  geom_errorbarh(aes(y = MeanlnHg, xmin = MeanC - SDC, xmax = MeanC + SDC), height = 0.1)+
  ggtitle('All species lnHg vs C') +
  xlab(' δ13C') +
  ylab('lnHg') +
  theme_classic() 
#Data analysis - Nitrogen and carbon ----
#create scatterplot for all species
#I had to use chatGPT to correct this code, on a side note it was wildly helpful
#When you put the data frame first you can name columns of that frame without frame$column
#with ellipses
ggplot(thule_birds, aes(x = carbon, y = nitrogen, colour = as.factor(sp_code))) +
  geom_point() +
  stat_ellipse(aes(group = as.factor(sp_number)), alpha = 0.7) +
  ggtitle('N vs C') +
  xlab('Carbon') +
  ylab('δ15N') +
  theme_classic()
#with markers for average
ggplot(thule_birds, aes(x = carbon, y = nitrogen, colour = as.factor(sp_code))) +
  geom_point() +
  stat_ellipse(aes(group = as.factor(sp_code)), alpha = 0.7, size = 1) +
  geom_point(data = SI_av_data, aes(x = MeanC, y = MeanN), color = "black", shape = 15, size = 3) +
  ggtitle('N vs C') +
  xlab('δ13C') +
  ylab('δ15N ') +
  theme_classic()

#with linear regressions
ggplot(thule_birds, aes(x = carbon, y = nitrogen, colour = as.factor(sp_code))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +  # Add linear regression lines without confidence intervals
  ggtitle('N vs C') +
  xlab('Carbon') +
  ylab('Nitrogen') +
  theme_classic()

#Kruskal-Wallis differences between species N

kruskal.test(nitrogen ~ sp_code, data = thule_birds)

#group comparisons

pairwise.wilcox.test(thule_birds$nitrogen, thule_birds$sp_code,
                     p.adjust.method = "bonferroni") 

#Kruskal-Wallis differences between species C

kruskal.test(carbon ~ sp_code, data = thule_birds)

#group comparisons

pairwise.wilcox.test(thule_birds$carbon, thule_birds$sp_code,
                     p.adjust.method = "bonferroni")  

#Puffin CN comparison
ggplot(ATPU, aes(x=carbon, y=nitrogen, colour = as.factor(year))) +
  geom_point() +
  stat_ellipse(aes(group = as.factor(year)), alpha = 0.7) +
  ggtitle('ATPU N vs C') +
  xlab('δ13C') +
  ylab('Nitrogen') +
  theme_classic()

#Guillemot CN comparison
ggplot(BLGU, aes(x=carbon, y=nitrogen, colour = as.factor(year))) +
  geom_point() +
  stat_ellipse(aes(group = as.factor(year)), alpha = 0.7) +
  ggtitle('BLGU N vs C') +
  xlab('Carbon') +
  ylab('Nitrogen') +
  theme_classic()

#Kittiwake CN comparison
ggplot(BLKI, aes(x=carbon, y=nitrogen, colour = as.factor(year))) +
  geom_point() +
  stat_ellipse(aes(group = as.factor(year)), alpha = 0.7) +
  ggtitle('BLKI N vs C') +
  xlab('Carbon') +
  ylab('Nitrogen') +
  theme_classic()

#Dovekie CN comparison
ggplot(DOVE, aes(x = carbon, y = nitrogen, colour = as.factor(year))) +
  geom_point() +
  stat_ellipse(aes(group = as.factor(year)), alpha = 0.7) +
  ggtitle('DOVE N vs C') +
  xlab('Carbon') +
  ylab('Nitrogen') +
  theme_classic()

#Murre CN comparison
ggplot(TBMU, aes(x=carbon, y=nitrogen, colour = as.factor(year))) +
  geom_point() +
  stat_ellipse(aes(group = as.factor(year)), alpha = 0.7) +
  ggtitle('TBMU N vs C') +
  xlab('Carbon') +
  ylab('Nitrogen') +
  theme_classic()



#Data analysis - Hg across species ----

#box plot

boxplot(Hg~sp_code,data=thule_birds, main="Hg by species",
        xlab="Species", ylab="Hg (mg/kg)")

library("ggpubr")
ggboxplot(thule_birds,
          x = "sp_code",
          y = "Hg",
          color = "sp_code",
          ylab = "Total Hg (ng/g ww)",
          xlab = "Species",
          legend = "none")
ggboxplot(thule_birds,
          x = "sp_code",
          y = "Hg",
          fill = "sp_code",  # Use fill instead of color
          ylab = "Total Hg (ng/g ww)",
          xlab = "Species",
          legend = "none") +
  geom_hline(yintercept = c(200, 1000), linetype = "dashed", color = "red")


#Kruskal-Wallis test

kruskal.test(Hg ~ sp_code, data = thule_birds)

#group comparisons

pairwise.wilcox.test(thule_birds$Hg, thule_birds$sp_code,
                     p.adjust.method = "bonferroni")  


#Data analysis - Hg over time ----
#Temporal mercury all species
ggplot(thule_birds, aes(x = year, y = Hg, colour = as.factor(sp_code))) +
  geom_point(position = position_jitter(width = 0.5, height = 0), alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  ggtitle('Temporal Hg') +
  xlab('Year') +
  ylab('Hg ppb') +
  theme_classic()+
  facet_wrap(~sp_code, scales = 'free_y')
#temp Hg with dot for mean of each year ##add CI
ggplot(thule_birds, aes(x = year, y = logHg, colour = as.factor(sp_code))) +
  geom_point(position = position_jitter(width = 0.5, height = 0), alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  stat_summary(fun = mean, geom = "point", shape = 8, size = 3, colour = "black", position = position_dodge(width = 0.5)) +
  ggtitle('Temporal Hg') +
  xlab('Year') +
  ylab('logHg ppb') +
  theme_classic() +
  facet_wrap(~sp_code, scales = 'free_y')

#average Hg linear regression
ggplot(year_av_data, aes(x = year, y = MeanHg, colour = as.factor(sp_code))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  geom_errorbar(aes(ymin = MeanHg - SDHg, ymax = MeanHg + SDHg), 
                position = position_dodge(width = 0.5), width = 0.2) +  # Add error bars
  ggtitle('Temporal Hg') +
  xlab('Year') +
  ylab('Hg ppb') +
  theme_classic() +
  facet_wrap(~sp_code, scales = 'free_y')
TBMUlr <- lm(lnHg ~ year, data = TBMU)
summary (TBMUlr)
BLGUlr <- lm(lnHg ~ year, data = BLGU)
summary (BLGUlr)
#Puffin temporal trend
#linear regression
shapiro.test(ATPU$Hg) #non normal
shapiro.test(ATPU$logHg) #normal!
ATPUHglinreg = lm(formula = Hg ~ year1,
                data = ATPU)
summary(ATPUHglinreg)


# Scatterplot with linear regression line
ggplot(ATPU, aes(x = year, y = Hg)) +
  geom_point(position = position_jitter(width = 0.3), color = 'blue', alpha = 0.7) +
  geom_line(aes(x = year, y = predict(ATPUHglinreg, newdata = ATPU)), colour = 'black') +
  ggtitle('ATPU Temporal Hg') +
  xlab('Year') +
  ylab('Hg ppb') +
  theme_classic() +
  scale_x_continuous(breaks = c(2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019))


  #Guillemot temporal trend
#linear regression
shapiro.test(BLGU$Hg) #non normal
shapiro.test(BLGU$logHg) #normal!
BLGUHglinreg = lm(formula = lnHg ~ year1,
                  data = BLGU)
summary(BLGUHglinreg)
BLGUHglinreg2 = lmer(formula = logHg ~ year + (1|year),
                  data = BLGU)
summary(BLGUHglinreg2)
BLGUHglinreg3 = lmer(formula = logHg ~ 1 + (1|year),
                     data = BLGU)
summary(BLGUHglinreg3)
AIC(BLGUHglinreg3)
AIC(BLGUHglinreg2)
AIC(BLGUHglinreg)
# Scatterplot with linear regression line
ggplot(BLGU, aes(x = year, y = logHg)) +
  geom_point(position = position_jitter(width = 0.3), color = 'blue', alpha = 0.7) +
  geom_line(aes(x = year, y = predict(BLGUHglinreg, newdata = BLGU)), colour = 'red') +
  geom_line(aes(x = year, y = predict(BLGUHglinreg2, newdata = BLGU)), colour = 'black') +
  geom_line(aes(x = year, y = predict(BLGUHglinreg3, newdata = BLGU)), colour = 'green') +
  ggtitle('BLGU Temporal Hg') +
  xlab('Year') +
  ylab('Hg ppb') +
  theme_classic() +
  scale_x_continuous(breaks = c(2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020,2021,2022,2023))

#Kittiwake temporal trend
#linear regression
shapiro.test(BLKI$Hg) #non normal
shapiro.test(BLKI$logHg) #normal!
BLKIHglinreg = lm(formula = logHg ~ year1,
                  data = BLKI)
summary(BLKIHglinreg)

# Scatterplot with linear regression line
ggplot(BLKI, aes(x = year, y = Hg)) +
  geom_point(position = position_jitter(width = 0.3), color = 'blue', alpha = 0.7) +
  geom_line(aes(x = year, y = predict(BLKIHglinreg, newdata = BLKI)), colour = 'black') +
  ggtitle('BLKI Temporal Hg') +
  xlab('Year') +
  ylab('Hg ppb') +
  theme_classic() +
  scale_x_continuous(breaks = c(2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020,2021,2022,2023))

#DOVE temporal trend
#linear regression
shapiro.test(DOVE$Hg) #non normal
shapiro.test(DOVE$logHg) #normal!
DOVEHglinreg = lm(formula = logHg ~ year1,
                  data = DOVE)
summary(DOVEHglinreg)

# Scatterplot with linear regression line
ggplot(DOVE, aes(x = year, y = Hg)) +
  geom_point(position = position_jitter(width = 0.3), color = 'blue', alpha = 0.7) +
  geom_line(aes(x = year, y = predict(DOVEHglinreg, newdata = DOVE)), colour = 'black') +
  ggtitle('DOVE Temporal Hg') +
  xlab('Year') +
  ylab('Hg ppb') +
  theme_classic() +
  scale_x_continuous(breaks = c(2010, 2019,2022,2023))

#Murre temporal trend
#linear regression
shapiro.test(TBMU$Hg) #non normal
shapiro.test(TBMU$logHg) #normal!
TBMUHglinreg = lm(formula = logHg ~ year1,
                  data = TBMU)
summary(TBMUHglinreg)

# Scatterplot with linear regression line
ggplot(TBMU, aes(x = year, y = Hg)) +
  geom_point(position = position_jitter(width = 0.3), color = 'blue', alpha = 0.7) +
  geom_line(aes(x = year, y = predict(TBMUHglinreg, newdata = TBMU)), colour = 'black') +
  ggtitle('TBMU Temporal Hg') +
  xlab('Year') +
  ylab('Hg ppb') +
  theme_classic() +
  scale_x_continuous(breaks = c(2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020,2021,2022,2023))

#linear-mixed effects model to assess temporal mercury with species as random effect

Allsp_linreg = lmer(logHg ~ year1 + (1|sp_code), data=thule_birds)
summary(Allsp_linreg)
report (Allsp_linreg)
# Temporal assessment of stable isotopes----

#Carbon
tempC <- lm(formula = carbon ~ year1, data = thule_birds)
summary (tempC)
ggplot(thule_birds, aes(x = year, y = carbon, colour = as.factor(sp_code))) +
  geom_point(position = position_jitter(width = 0.5, height = 0), alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE) +
  ggtitle('Temporal C') +
  xlab('Year') +
  ylab('carbon') +
  theme_classic()
#Carbon one regression
ggplot(thule_birds, aes(x = year, y = carbon, colour = as.factor(sp_code))) +
  geom_point(position = position_jitter(width = 0.5, height = 0), alpha = 0.7) +
  geom_smooth (method = 'lm', aes(group = 1), se = FALSE, color = "black") +
  ggtitle('Temporal C') +
  xlab('Year') +
  ylab('Carbon') +
  theme_classic()

ggplot(thule_birds, aes(x = year, y = carbon, colour = as.factor(sp_code))) +
  geom_point(position = position_jitter(width = 0.5, height = 0), alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  stat_summary(fun = mean, geom = "point", shape = 8, size = 3, colour = "black", position = position_dodge(width = 0.5)) +
  ggtitle('Temporal C') +
  xlab('Year') +
  ylab('C') +
  theme_classic() +
  facet_wrap(~sp_code, scales = 'free_y')

tempN <- lm(formula = nitrogen ~ year1, data = thule_birds)

tempN <- lmer(nitrogen ~ year1 + (1|sp_code), data = thule_birds)
summary (tempN)
#carbon linear regressions
shapiro.test(ATPU$carbon) #normal
ATPUClinreg = lm(formula = carbon ~ year1,
                  data = ATPU)
summary(ATPUClinreg)
shapiro.test(BLGU$carbon) #normal
BLGUClinreg = lm(formula = carbon ~ year1,
                 data = BLGU)
summary(BLGUClinreg)
shapiro.test(BLKI$carbon) #normal
BLKIClinreg = lm(formula = carbon ~ year1,
                 data = BLKI)
summary(BLKIClinreg)
shapiro.test(DOVE$carbon) 
DOVEClinreg = lm(formula = carbon ~ year1,
                 data = DOVE)
summary(DOVEClinreg)
shapiro.test(TBMU$carbon)
TBMUClinreg = lm(formula = carbon ~ year1,
                 data = TBMU)
summary(TBMUClinreg)
#nitrogen
ggplot(thule_birds, aes(x = year, y = nitrogen, colour = as.factor(sp_code))) +
  geom_point(position = position_jitter(width = 0.5, height = 0), alpha = 0.7) +
  geom_smooth(method = "lm", se = TRUE) +
  ggtitle('Temporal N') +
  xlab('Year') +
  ylab('N') +
  theme_classic()
#Carbon one regression
ggplot(thule_birds, aes(x = year, y = nitrogen, colour = as.factor(sp_code))) +
  geom_point(position = position_jitter(width = 0.5, height = 0), alpha = 0.7) +
  geom_smooth (method = 'lm', aes(group = 1), se = FALSE, color = "black") +
  ggtitle('Temporal N') +
  xlab('Year') +
  ylab('Nitrogen') +
  theme_classic()

ggplot(thule_birds, aes(x = year, y = nitrogen, colour = as.factor(sp_code))) +
  geom_point(position = position_jitter(width = 0.5, height = 0), alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  stat_summary(fun = mean, geom = "point", shape = 8, size = 3, colour = "black", position = position_dodge(width = 0.5)) +
  ggtitle('Temporal N') +
  xlab('Year') +
  ylab('N') +
  theme_classic() +
  facet_wrap(~sp_code, scales = 'free_y')
#nitrogen smoothed
ggplot(thule_birds, aes(x = year, y = nitrogen, colour = as.factor(sp_code))) +
  geom_point(position = position_jitter(width = 0.5, height = 0), alpha = 0.7) +
  geom_smooth(method = 'lm', aes(group = 1), se = TRUE) +
  ggtitle('Temporal N') +
  xlab('Year') +
  ylab('Nitrogen') +
  theme_classic()
#nitrogen linear regressions
shapiro.test(ATPU$nitrogen) #not normal
ATPUNlinreg = lm(formula = nitrogen ~ year1,
                 data = ATPU)
summary(ATPUNlinreg)
AIC(ATPUNlinreg)
shapiro.test(BLGU$nitrogen) #normal
BLGUNlinreg = lm(formula = nitrogen ~ year1,
                 data = BLGU)
summary(BLGUNlinreg)
AIC(BLGUNlinreg)
shapiro.test(BLKI$nitrogen) #normal
BLKINlinreg = lm(formula = nitrogen ~ year1,
                 data = BLKI)
summary(BLKINlinreg)
AIC(BLKINlinreg)
shapiro.test(DOVE$nitrogen) #not normal
DOVENlinreg = lm(formula = nitrogen ~ year1,
                 data = DOVE)
summary(DOVENlinreg)
AIC(DOVENlinreg)
shapiro.test(TBMU$nitrogen) #barely normal
TBMUNlinreg = lm(formula = nitrogen ~ year1,
                 data = TBMU)
summary(TBMUNlinreg)
AIC(TBMUNlinreg)