Pastos marinos

library(readxl)

p.pastos <- read_excel("/media/m/Margarita Loría Naranjo/CIMAR/CARICOMP/caricomp pastos/Datos completos para análisis.xlsx", sheet = "Parametros pastos 2sp")

p.pastos <- p.pastos[,!is.na(names(p.pastos))]

p.pastos$Date <- as.Date(p.pastos$Date, origin = "1899-12-30")

p.pastos$year <- sapply(strsplit(as.character(p.pastos$Date), "-",fixed=T), "[[", 1)

p.pastos$month <- sapply(strsplit(as.character(p.pastos$Date), "-",fixed=T), "[[", 2)


p.pastos$year <- as.factor(p.pastos$year)
p.pastos$month <- as.factor(p.pastos$month)

pSy <- p.pastos[p.pastos$Specie == "Syringodium filiforme",]
pSy <- pSy[pSy$Biomass < 400,]

pTh <- p.pastos[p.pastos$Specie != "Syringodium filiforme",]


library(nlme)

Biomass Syringodium filiforme

linear mixed models with month as predictor and year as random effect

m1 <- lme(Biomass ~ month, random = ~1|year, data = pSy)
anova(m1)

##             numDF denDF  F-value p-value
## (Intercept)     1   115 5.782544  0.0178
## month          11   115 0.793912  0.6457

plot(m1)

agbiom <- aggregate(pSy$Biomass, by = list(pSy$month), mean)
agbiomsd <- aggregate(pSy$Biomass, by = list(pSy$month), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agbiom <- data.frame(agbiom, agbiomsd[,2])

names(agbiom) <- c("Month", "Biomass", "sd")

library(ggplot2)

ggplot(agbiom, aes(x=Month, y= Biomass))  +
  geom_errorbar(width=.2, aes(ymin= Biomass - sd, ymax = Biomass + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Biomass (g/m2)") +
  xlab("Month") 

linear mixed models with year as predictor and month as random effect

m2 <- lme(Biomass ~ year, random = ~1|month, data = pSy)

anova(m2)

##             numDF denDF   F-value p-value
## (Intercept)     1   115 19.566193  <.0001
## year            5   115  1.807969  0.1167

plot(m2)

agbiom <- aggregate(pSy$Biomass, by = list(pSy$year), mean)
agbiomsd <- aggregate(pSy$Biomass, by = list(pSy$year), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agbiom <- data.frame(agbiom, agbiomsd[,2])

names(agbiom) <- c("Month", "Biomass", "sd")

library(ggplot2)

ggplot(agbiom, aes(x=Month, y= Biomass))  +
  geom_errorbar(width=.2, aes(ymin= Biomass - sd, ymax = Biomass + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Biomass (g/m2)") +
  xlab("Year") 

Biomass Thalassia testudinum

linear mixed models with month as predictor and year as random effect

m1 <- lme(Biomass ~ month, random = ~1|year, data = pTh)
anova(m1)

##             numDF denDF   F-value p-value
## (Intercept)     1   151 1408.8705  <.0001
## month          11   151    2.6399  0.0041

plot(m1)

agbiom <- aggregate(pTh$Biomass, by = list(pTh$month), mean)
agbiomsd <- aggregate(pTh$Biomass, by = list(pTh$month), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agbiom <- data.frame(agbiom, agbiomsd[,2])

names(agbiom) <- c("Month", "Biomass", "sd")

library(ggplot2)

ggplot(agbiom, aes(x=Month, y= Biomass))  +
  geom_errorbar(width=.2, aes(ymin= Biomass - sd, ymax = Biomass + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Biomass (g/m2)") +
  xlab("Month") 

linear mixed models with year as predictor and month as random effect

m2 <- lme(Biomass ~ year, random = ~1|month, data = pTh)

anova(m2)

##             numDF denDF  F-value p-value
## (Intercept)     1   151 904.8750  <.0001
## year            6   151   1.8085   0.101

plot(m2)

agbiom <- aggregate(pTh$Biomass, by = list(pTh$year), mean)
agbiomsd <- aggregate(pTh$Biomass, by = list(pTh$year), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agbiom <- data.frame(agbiom, agbiomsd[,2])

names(agbiom) <- c("Month", "Biomass", "sd")

library(ggplot2)

ggplot(agbiom, aes(x=Month, y= Biomass))  +
  geom_errorbar(width=.2, aes(ymin= Biomass - sd, ymax = Biomass + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Biomass (g/m2)") +
  xlab("Year") 

Productivity

linear mixed models with year as predictor and month as random effect

pro <- read_excel("/media/m/Margarita Loría Naranjo/CIMAR/CARICOMP/caricomp pastos/Datos completos para análisis.xlsx", sheet = "Productividad")

names(pro)[4] <- "Prod"

pro$Date <- as.Date(pro$Date, origin = "1899-12-30")

pro$year <- sapply(strsplit(as.character(pro$Date), "-",fixed=T), "[[", 1)

pro$month <- sapply(strsplit(as.character(pro$Date), "-",fixed=T), "[[", 2)

m4 <- lme(Prod ~ year, random = ~1|month, data = pro)

anova(m4)

##             numDF denDF   F-value p-value
## (Intercept)     1   112 137.71083  <.0001
## year            6   112   9.05566  <.0001

plot(m4)

qqnorm(resid(m4))
qqline(resid(m4))

agpro <- aggregate(pro$Prod, by = list(pro$year), mean)
agprosd <- aggregate(pro$Prod, by = list(pro$year), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agpro <- data.frame(agpro, agprosd[,2])

names(agpro) <- c("Month", "Prod", "sd")


ggplot(agpro, aes(x=Month, y= Prod))  +
  geom_errorbar(width=.2, aes(ymin= Prod - sd, ymax = Prod + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Productivity (g/m2/day)") +
  xlab("Year") 

linear mixed models with month as predictor and year as random effect

pro <- read_excel("/media/m/Margarita Loría Naranjo/CIMAR/CARICOMP/caricomp pastos/Datos completos para análisis.xlsx", sheet = "Productividad")

names(pro)[4] <- "Prod"

pro$Date <- as.Date(pro$Date, origin = "1899-12-30")

pro$year <- sapply(strsplit(as.character(pro$Date), "-",fixed=T), "[[", 1)

pro$month <- sapply(strsplit(as.character(pro$Date), "-",fixed=T), "[[", 2)

m5 <- lme(Prod ~ month, random = ~1|year, data = pro)

anova(m5)

##             numDF denDF  F-value p-value
## (Intercept)     1   112 52.26354  <.0001
## month          10   112  4.07697   1e-04

plot(m5)

qqnorm(resid(m5))
qqline(resid(m5))

agpro <- aggregate(pro$Prod, by = list(pro$month), mean)
agprosd <- aggregate(pro$Prod, by = list(pro$month), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agpro <- data.frame(agpro, agprosd[,2])

names(agpro) <- c("Month", "Prod", "sd")


ggplot(agpro, aes(x=Month, y= Prod))  +
  geom_errorbar(width=.2, aes(ymin= Prod - sd, ymax = Prod + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Productivity (g/m2/day)") +
  xlab("Month") 

Change rate

linear mixed models with year as predictor and month as random effect

rec <- read_excel("/media/m/Margarita Loría Naranjo/CIMAR/CARICOMP/caricomp pastos/Datos completos para análisis.xlsx", sheet = "Tasa de recambio")

names(rec)[4] <- "Change"

rec$Date <- as.Date(rec$Date, origin = "1899-12-30")

rec$year <- sapply(strsplit(as.character(rec$Date), "-",fixed=T), "[[", 1)

rec$month <- sapply(strsplit(as.character(rec$Date), "-",fixed=T), "[[", 2)

m6 <- lme(Change ~ year, random = ~1|month, data = rec)

anova(m6)

##             numDF denDF   F-value p-value
## (Intercept)     1   112 188.34833  <.0001
## year            6   112  10.10165  <.0001

plot(m6)

qqnorm(resid(m6))
qqline(resid(m6))

agrec <- aggregate(rec$Change, by = list(rec$year), mean)
agrecsd <- aggregate(rec$Change, by = list(rec$year), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agrec <- data.frame(agrec, agrecsd[,2])

names(agrec) <- c("Year", "Change", "sd")


ggplot(agrec, aes(x=Year, y= Change))  +
  geom_errorbar(width=.2, aes(ymin= Change - sd, ymax = Change + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Change rate (%/day)") +
  xlab("Year") 

linear mixed models with month as predictor and year as random effect

m7 <- lme(Change ~ month, random = ~1|year, data = rec)

anova(m7)

##             numDF denDF   F-value p-value
## (Intercept)     1   112 164.53607  <.0001
## month          10   112   6.07953  <.0001

plot(m7)

qqnorm(resid(m7))
qqline(resid(m7))

agrec <- aggregate(rec$Change, by = list(rec$month), mean)
agrecsd <- aggregate(rec$Change, by = list(rec$month), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agrec <- data.frame(agrec, agrecsd[,2])

names(agrec) <- c("Month", "Change", "sd")


ggplot(agrec, aes(x=Month, y= Change))  +
  geom_errorbar(width=.2, aes(ymin= Change - sd, ymax = Change + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Change rate (%/day)") +
  xlab("Year") 

Density

linear mixed models with year as predictor and month as random effect

den <- read_excel("/media/m/Margarita Loría Naranjo/CIMAR/CARICOMP/caricomp pastos/Datos completos para análisis.xlsx", sheet = "Densidad")

names(den)[4] <- "Den"

den$Date <- as.Date(den$Date, origin = "1899-12-30")

den$year <- sapply(strsplit(as.character(den$Date), "-",fixed=T), "[[", 1)

den$month <- sapply(strsplit(as.character(den$Date), "-",fixed=T), "[[", 2)

m8 <- lme(Den ~ year, random = ~1|month, data = den)

anova(m8)

##             numDF denDF   F-value p-value
## (Intercept)     1   112 129.64789  <.0001
## year            6   112   8.45669  <.0001

plot(m8)

qqnorm(resid(m8))
qqline(resid(m8))

agden <- aggregate(den$Den, by = list(den$year), mean)
agdensd <- aggregate(den$Den, by = list(den$year), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agden <- data.frame(agden, agdensd[,2])

names(agden) <- c("Year", "Den", "sd")

ggplot(agden, aes(x=Year, y= Den))  +
  geom_errorbar(width=.2, aes(ymin= Den - sd, ymax = Den + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Density (shoot/m2)") +
  xlab("Year") 

linear mixed models with month as predictor and year as random effect

m9 <- lme(Den ~ month, random = ~1|year, data = den)

anova(m9)

##             numDF denDF  F-value p-value
## (Intercept)     1   112 49.28314  <.0001
## month          10   112  6.06628  <.0001

plot(m9)

qqnorm(resid(m9))
qqline(resid(m9))

agden <- aggregate(den$Den, by = list(den$month), mean)
agdensd <- aggregate(den$Den, by = list(den$month), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agden <- data.frame(agden, agdensd[,2])

names(agden) <- c("Month", "Den", "sd")

ggplot(agden, aes(x=Month, y= Den))  +
  geom_errorbar(width=.2, aes(ymin= Den - sd, ymax = Den + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Density (shoot/m2)") +
  xlab("Month") 

Leaf Area Index

linear mixed models with year as predictor and month as random effect

lai <- read_excel("/media/m/Margarita Loría Naranjo/CIMAR/CARICOMP/caricomp pastos/Datos completos para análisis.xlsx", sheet = "Indice area foliar")

names(lai)[4] <- "lai"

lai$Date <- as.Date(lai$Date, origin = "1899-12-30")

lai$year <- sapply(strsplit(as.character(lai$Date), "-",fixed=T), "[[", 1)

lai$month <- sapply(strsplit(as.character(lai$Date), "-",fixed=T), "[[", 2)

m10 <- lme(lai ~ year, random = ~1|month, data = lai)

anova(m10)

##             numDF denDF  F-value p-value
## (Intercept)     1    89 35.74425  <.0001
## year            5    89  2.81058   0.021

plot(m10)

qqnorm(resid(m10))
qqline(resid(m10))

aglai <- aggregate(lai$lai, by = list(lai$year), mean)
aglaisd <- aggregate(lai$lai, by = list(lai$year), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

aglai <- data.frame(aglai, aglaisd[,2])

names(aglai) <- c("Year", "lai", "sd")

ggplot(aglai, aes(x=Year, y= lai))  +
  geom_errorbar(width=.2, aes(ymin= lai - sd, ymax = lai + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Leaf Area Index") +
  xlab("Year") 

linear mixed models with month as predictor and year as random effect

m11 <- lme(lai ~ month, random = ~1|year, data = lai)

anova(m11)

##             numDF denDF   F-value p-value
## (Intercept)     1    89 26.140627  <.0001
## month          10    89  4.792835  <.0001

plot(m11)

qqnorm(resid(m11))
qqline(resid(m11))

aglai <- aggregate(lai$lai, by = list(lai$month), mean)
aglaisd <- aggregate(lai$lai, by = list(lai$month), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

aglai <- data.frame(aglai, aglaisd[,2])

names(aglai) <- c("Month", "lai", "sd")

ggplot(aglai, aes(x=Month, y= lai))  +
  geom_errorbar(width=.2, aes(ymin= lai - sd, ymax = lai + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Leaf Area Index") +
  xlab("Month") 

Leaf Area

linear mixed models with year as predictor and month as random effect

la <- read_excel("/media/m/Margarita Loría Naranjo/CIMAR/CARICOMP/caricomp pastos/Datos completos para análisis.xlsx", sheet = "Area foliar")

names(la)[5] <- "la"

la$Date <- as.Date(la$Date, origin = "1899-12-30")

la$year <- sapply(strsplit(as.character(la$Date), "-",fixed=T), "[[", 1)

la$month <- sapply(strsplit(as.character(la$Date), "-",fixed=T), "[[", 2)

m12 <- lme(la ~ year, random = ~1|month, data = la)

anova(m12)

##             numDF denDF  F-value p-value
## (Intercept)     1   304 95.15876  <.0001
## year            5   304  5.18906   1e-04

plot(m12)

qqnorm(resid(m12))
qqline(resid(m12))

agla <- aggregate(la$la, by = list(la$year), mean)
aglasd <- aggregate(la$la, by = list(la$year), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agla <- data.frame(agla, aglasd[,2])

names(agla) <- c("Year", "la", "sd")

ggplot(agla, aes(x=Year, y= la))  +
  geom_errorbar(width=.2, aes(ymin= la - sd, ymax = la + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Leaf Area") +
  xlab("Year") 

linear mixed models with month as predictor and year as random effect

m11 <- lme(la ~ month, random = ~1|year, data = la)

anova(m11)

##             numDF denDF  F-value p-value
## (Intercept)     1   304 33.10026  <.0001
## month           9   304  4.24415  <.0001

plot(m11)

qqnorm(resid(m11))
qqline(resid(m11))

agla <- aggregate(la$la, by = list(la$month), mean)
aglasd <- aggregate(la$la, by = list(la$month), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agla <- data.frame(agla, aglasd[,2])

names(agla) <- c("Month", "la", "sd")

ggplot(agla, aes(x=Month, y= la))  +
  geom_errorbar(width=.2, aes(ymin= la - sd, ymax = la + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Leaf Area") +
  xlab("Month") 

Parámetros ambientales

pa <- read_excel("/media/m/Margarita Loría Naranjo/CIMAR/CARICOMP/caricomp pastos/Datos completos para análisis.xlsx", sheet = "Parámetros ambientales")

rain <- pa[,1:2]

names(rain)[2] <- "rain"

rain$Date <- as.Date(rain$Date, origin = "1899-12-30")

rain <- rain[!is.na(rain$rain),]

rain$year <- sapply(strsplit(as.character(rain$Date), "-",fixed=T), "[[", 1)

rain$month <- sapply(strsplit(as.character(rain$Date), "-",fixed=T), "[[", 2)

agrain <- aggregate(rain$rain, by = list(rain$month), mean)
agrainsd <- aggregate(rain$rain, by = list(rain$month), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agrain <- data.frame(agrain, agrainsd[,2])

names(agrain) <- c("Month", "rain", "sd")

ggplot(agrain, aes(x=Month, y= rain))  +
  geom_errorbar(width=.2, aes(ymin= rain - sd, ymax = rain + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Rain fall (mm)") +
  xlab("Month") 

Temperature

temp <- pa[,3:4]

names(temp)[2] <- "temp"

temp$Date <- as.Date(temp$Date, origin = "1899-12-30")

temp <- temp[!is.na(temp$temp),]

temp$year <- sapply(strsplit(as.character(temp$Date), "-",fixed=T), "[[", 1)

temp$month <- sapply(strsplit(as.character(temp$Date), "-",fixed=T), "[[", 2)

agtemp <- aggregate(temp$temp, by = list(temp$month), mean)
agtempsd <- aggregate(temp$temp, by = list(temp$month), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agtemp <- data.frame(agtemp, agtempsd[,2])

names(agtemp) <- c("Month", "temp", "sd")

ggplot(agtemp, aes(x=Month, y= temp))  +
  geom_errorbar(width=.2, aes(ymin= temp - sd, ymax = temp + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Temperature (C)") +
  xlab("Month") 

Salinity

sal <- pa[,5:8]

names(sal)[4] <- "sal"

sal$Date <- as.Date(sal$Date, origin = "1899-12-30")

sal <- sal[!is.na(sal$sal),]

sal$year <- sapply(strsplit(as.character(sal$Date), "-",fixed=T), "[[", 1)

sal$month <- sapply(strsplit(as.character(sal$Date), "-",fixed=T), "[[", 2)

agsal <- aggregate(sal$sal, by = list(sal$month), mean)
agsalsd <- aggregate(sal$sal, by = list(sal$month), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agsal <- data.frame(agsal, agsalsd[,2])

names(agsal) <- c("Month", "sal", "sd")

ggplot(agsal, aes(x=Month, y= sal))  +
  geom_errorbar(width=.2, aes(ymin= sal - sd, ymax = sal + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Salinity (ppt)") +
  xlab("Month") 

Secchi

sec <- pa[,9:12]

names(sec)[4] <- "sec"

sec$Date <- as.Date(sec$Date, origin = "1899-12-30")

sec <- sec[!is.na(sec$sec),]

sec$year <- sapply(strsplit(as.character(sec$Date), "-",fixed=T), "[[", 1)

sec$month <- sapply(strsplit(as.character(sec$Date), "-",fixed=T), "[[", 2)

agsec <- aggregate(sec$sec, by = list(sec$month), mean)
agsecsd <- aggregate(sec$sec, by = list(sec$month), function(x) sd(x)/sqrt(length(x)))
# agbiomsd <- aggregate(p.pastos$Biomass, by = list(p.pastos$month), sd)

agsec <- data.frame(agsec, agsecsd[,2])

names(agsec) <- c("Month", "sec", "sd")

ggplot(agsec, aes(x=Month, y= sec))  +
  geom_errorbar(width=.2, aes(ymin= sec - sd, ymax = sec + sd)) +
  geom_point( shape = 21, size= 4, fill="white") + ylab("Secchi (m)") +
  xlab("Month")