rm(list=ls())
setwd('D:/Users/m.bartrons/Dropbox/Madison_Postdoc/Spatial_gradient/ZOOPLANKTON/')
T1=read.table("Transect1_data.txt", sep="\t")

#########################################################################
### SPATIAL MAP                                                      
#########################################################################

library(sp)
library(rgdal)
library(maptools)
library(mapplots)

Myvatn <- readShapePoly("D:/Users/m.bartrons/Dropbox/Madison_Postdoc/Iceland/Myvatn_map/Myvatn_WSGUTM28.shp", CRS("+proj=utm +zone=28"),IDvar=NULL)

coordsT1<- cbind(T1$longitude, T1$latitude)
colnames(coordsT1)<- c("long", "lat")
WGScoordsT1 <- SpatialPoints(coordsT1, proj4string=CRS("+proj=longlat +datum=WGS84"))
UTM28NCoordsT1 <- spTransform(WGScoordsT1, CRS("+proj=utm +zone=28"))
newCoordsT1<- UTM28NCoordsT1@coords
PT1<-cbind(newCoordsT1,T1)
ZooT1=subset(PT1, select=c("long","lat","sampleid","dalo_nL","cycl_nL","naup_nL","chsp_nL"))

library(reshape)
ZooT1b <- melt(ZooT1, id=(c("long","lat","sampleid")))
ZooT1xyz <- make.xyz(x=ZooT1b$long,y=ZooT1b$lat,z=ZooT1b$value,group=ZooT1b$variable)

plot(Myvatn)
draw.pie(z=ZooT1xyz$z,x=ZooT1xyz$x, y=ZooT1xyz$y, scale=T, radius=450,col=c("palegreen2","red2","sienna1","slateblue2"))
legend.pie("top",labels=c("Daphnia","Cyclops","Nauplii","Chydorus"),col=c("palegreen2","red2","sienna1","slateblue2"), radius=400, bty='n',cex = 0.9)
title(main="Transect 1, 02-07-2012",cex= 0.2)

################################################################################################################################################
# analyze data: lmer (MLM)
################################################################################################################################################

T1z <- subset(PT1, select=c("long","lat","sampleid","sampledate","do","k","ph","sampledepth","wtemp","exp.wtemp","chl_ugl","log.chl_ugl","PCYV","log.PCYV","log.PCYVChl","logit.Oocystis","Oocystis","logit.Anabaena","logit.AnabaenaAO","logit.AnabaenaAOr","logit.AnabaenaOAC","logit.Diatoms","logit.OtherChlorophyta","logit.Chrysophyceae","log.Oocystisr","log.Anabaenar","log.Diatomsr","log.OtherChlorophytar","log.Chrysophyceaer","TN","NH4","Phosphate","Silica","SpCond","dist","horn","aspl_nL","naup_nL","cycl_nL","dalo_nL","alona_nL","eula_nL","chsp_nL","acha_nL","mahi_nL","sive_nL"))

#Normalize predictor variables:
for(i in 5:36){T1z[,i] <- (T1z[,i] - mean(T1z[,i]))/sd(T1z[,i])}  #T1
#for(i in 5:36){T2z[,i] <- (T2z[,i] - mean(T2z[,i]))/sd(T2z[,i])}  #T2 
#for(i in 5:36){T3z[,i] <- (T3z[,i] - mean(T3z[,i]))/sd(T3z[,i])}  #T3

#TT1z
T1zsp <- subset(T1z, select=c("long","lat","sampleid","sampledate","aspl_nL","naup_nL","cycl_nL","dalo_nL","alona_nL","eula_nL","chsp_nL","acha_nL","mahi_nL","sive_nL"))
T1zenv <- subset(T1z, select=c("long","lat","sampleid","sampledate","do","k","ph","sampledepth","wtemp","exp.wtemp","chl_ugl","log.chl_ugl","PCYV","log.PCYV","log.PCYVChl","logit.Oocystis","Oocystis","logit.Anabaena","logit.AnabaenaAOr","logit.AnabaenaAO","logit.AnabaenaOAC","logit.Diatoms","logit.OtherChlorophyta","logit.Chrysophyceae","log.Oocystisr","log.Anabaenar","log.Diatomsr","log.OtherChlorophytar","log.Chrysophyceaer","TN","NH4","Phosphate","Silica","SpCond","dist","horn"))
TT1zmelt <- melt(T1zsp, id=(c("long","lat","sampleid","sampledate")))
TT1z <- merge(T1zenv,TT1zmelt, by = c("long","lat","sampleid","sampledate"))
names(TT1z)[37] <- 'SPP' 
names(TT1z)[38] <- 'ABUNDANCE' 

#TT2z
T2zsp <- subset(T2z, select=c("long","lat","sampleid","sampledate","aspl_nL","kera_nL","naup_nL","cycl_nL","dalo_nL","alona_nL","eula_nL","chsp_nL","acha_nL","mahi_nL"))
T2zenv <- subset(T2z, select=c("long","lat","sampleid","sampledate","do","k","ph","sampledepth","wtemp","exp.wtemp","chl_ugl","log.chl_ugl","PCYV","log.PCYV","log.PCYVChl","logit.Oocystis","Oocystis","logit.Anabaena","logit.AnabaenaAOr","logit.AnabaenaAO","logit.AnabaenaOAC","logit.Diatoms","logit.OtherChlorophyta","logit.Chrysophyceae","log.Oocystisr","log.Anabaenar","log.Diatomsr","log.OtherChlorophytar","log.Chrysophyceaer","TN","NH4","Phosphate","Silica","SpCond","dist","horn"))
TT2zmelt <- melt(T2zsp, id=(c("long","lat","sampleid","sampledate")))
TT2z <- merge(T2zenv,TT2zmelt, by = c("long","lat","sampleid","sampledate"))
names(TT2z)[37] <- 'SPP' 
names(TT2z)[38] <- 'ABUNDANCE' 

#TT3z
T3zsp <- subset(T3z, select=c("long","lat","sampleid","sampledate","aspl_nL","kera_nL","naup_nL","cycl_nL","dalo_nL","alona_nL","eula_nL","chsp_nL","acha_nL","mahi_nL","sive_nL"))
T3zenv <- subset(T3z, select=c("long","lat","sampleid","sampledate","do","k","ph","sampledepth","wtemp","exp.wtemp","chl_ugl","log.chl_ugl","PCYV","log.PCYV","log.PCYVChl","logit.Oocystis","Oocystis","logit.Anabaena","logit.AnabaenaAOr","logit.AnabaenaAO","logit.AnabaenaOAC","logit.Diatoms","logit.OtherChlorophyta","logit.Chrysophyceae","log.Oocystisr","log.Anabaenar","log.Diatomsr","log.OtherChlorophytar","log.Chrysophyceaer","TN","NH4","Phosphate","Silica","SpCond","dist","horn"))
TT3zmelt <- melt(T3zsp, id=(c("long","lat","sampleid","sampledate")))
TT3z <- merge(T3zenv,TT3zmelt, by = c("long","lat","sampleid","sampledate"))
names(TT3z)[37] <- 'SPP' 
names(TT3z)[38] <- 'ABUNDANCE' 

#TT2z9
T2z9 <- T2z[c(1,3,5,7,9,11,13,15,17),]
T2z9sp <- subset(T2z9, select=c("long","lat","sampleid","sampledate","aspl_nL","kera_nL","naup_nL","cycl_nL","dalo_nL","alona_nL","eula_nL","chsp_nL","acha_nL","mahi_nL"))
T2z9env <- subset(T2z9, select=c("long","lat","sampleid","sampledate","do","k","ph","sampledepth","wtemp","exp.wtemp","chl_ugl","log.chl_ugl","PCYV","log.PCYV","log.PCYVChl","logit.Oocystis","Oocystis","logit.Anabaena","logit.AnabaenaAOr","logit.AnabaenaAO","logit.AnabaenaOAC","logit.Diatoms","logit.OtherChlorophyta","logit.Chrysophyceae","log.Oocystisr","log.Anabaenar","log.Diatomsr","log.OtherChlorophytar","log.Chrysophyceaer","TN","NH4","Phosphate","Silica","SpCond","dist","horn"))
TT2z9melt <- melt(T2z9sp, id=(c("long","lat","sampleid","sampledate")))
TT2z9 <- merge(T2z9env,TT2z9melt, by = c("long","lat","sampleid","sampledate"))
names(TT2z9)[37] <- 'SPP' 
names(TT2z9)[38] <- 'ABUNDANCE' 

TT1z$sampledepth2 <- TT1z$sampledepth^2
TT1z$log.ABUNDANCE <- log(TT1z$ABUNDANCE+0.5)

TT2z$sampledepth2 <- TT2z$sampledepth^2
TT2z$log.ABUNDANCE <- log(TT2z$ABUNDANCE+0.5)

TT3z$sampledepth2 <- TT3z$sampledepth^2
TT3z$log.ABUNDANCE <- log(TT3z$ABUNDANCE+0.5)

TT2z9$sampledepth2 <- TT2z9$sampledepth^2
TT2z9$log.ABUNDANCE <- log(TT2z9$ABUNDANCE+0.5)
TT2z9$log.ABUNDANCEstd <- TT2z9$log.ABUNDANCE - mean(TT2z9$log.ABUNDANCE)/sd(TT2z9$log.ABUNDANCE) 

#lme() because you can't include a correlation structure using lmer()
library(nlme)
library(arm)
library(car)
library(MASS)

mT1 <- lme(fixed = log.ABUNDANCE ~ sampledepth + k + log.chl_ugl + logit.Anabaena + logit.Oocystis, 
           random = list(~1|SPP,~ 0+log.chl_ugl|SPP,~ 0+logit.Anabaena|SPP),
           correlation = corExp(form = ~ long + lat, nugget=T), method = "ML",control=lmeControl(opt = "optim"), data=TT1z)

mT2 <- lme(fixed = log.ABUNDANCE ~ exp.wtemp +sampledepth + k + log.chl_ugl + logit.Anabaena,
            random = list(~1|SPP,~ 0+sampledepth|SPP,~ 0+log.chl_ugl|SPP,~ 0+logit.Anabaena|SPP),
            correlation = corExp(form = ~ long + lat, nugget=T), method = "ML",control=lmeControl(opt = "optim"), data=TT2z)

mT29 <- lme(fixed = log.ABUNDANCEstd ~ log.chl_ugl + logit.Anabaena,
           random = list(~1|SPP,~ 0+log.chl_ugl|SPP),
           correlation = corExp(form = ~ long + lat, nugget=T), method = "ML",control=lmeControl(opt = "optim"), data=TT2z9)

mT3 <- lme(fixed = log.ABUNDANCE ~ exp.wtemp + k + log.chl_ugl + logit.Oocystis, 
           random = list(~1|SPP,~ 0+log.chl_ugl|SPP,~ 0+logit.Oocystis|SPP),
           correlation = corExp(form = ~ long + lat, nugget=T), method = "ML",control=lmeControl(opt = "optim"), data=TT3z)

summary(mT1) #fixed+random
ranef(mT1) #fixed+random
intervals(mT1)

################################################################################################################################################
# analyze data: lmer (MLM)
# compute random effects pvalues 
################################################################################################################################################

##### T1
#full
mT1 <- lme(fixed = log.ABUNDANCE ~ sampledepth + k + log.chl_ugl + logit.Anabaena + logit.Oocystis, 
           random = list(~1|SPP,~ 0+log.chl_ugl|SPP,~ 0+logit.Anabaena|SPP),
           correlation = corExp(form = ~ long + lat, nugget=T), method = "ML",control=lmeControl(opt = "optim"), data=TT1z)

#without log.chl_ugl
mT1_1 <- lme(fixed = log.ABUNDANCE ~ sampledepth + k + log.chl_ugl + logit.Anabaena + logit.Oocystis, 
           random = list(~1|SPP,~ 0+logit.Anabaena|SPP),
           correlation = corExp(form = ~ long + lat, nugget=T), method = "ML",control=lmeControl(opt = "optim"), data=TT1z)

#without logit.Anabaena
mT1_2 <- lme(fixed = log.ABUNDANCE ~ sampledepth + k + log.chl_ugl + logit.Anabaena + logit.Oocystis, 
           random = list(~1|SPP,~ 0+log.chl_ugl|SPP),
           correlation = corExp(form = ~ long + lat, nugget=T), method = "ML",control=lmeControl(opt = "optim"), data=TT1z)

# if deviance in LMER is the same as (-2)*logLik, this is correct:
LLT1 <- c((-2)*logLik(mT1)[1], (-2)*logLik(mT1_1)[1], (-2)*logLik(mT1_2)[1])
mlm.pvalsT1 <- c(1-pchisq(LLT1[2] - LLT1[1],1), 1-pchisq(LLT1[3] - LLT1[1],1))
mlm.pvalsT1 #Table1


####################################################################################
## NUL MODEL
####################################################################################

nT1cn <- lme(fixed = log.ABUNDANCE ~ 1, random = ~1|SPP,
           correlation = corExp(form = ~ long + lat, nugget=T), method = "ML", data=TT1z) #AIC: 600 #<--- BEST MODEL BY FAR
nT1c <- lme(fixed = log.ABUNDANCE ~ 1, random = ~1|SPP,
            correlation = corExp(form = ~ long + lat, nugget=F), method = "ML", data=TT1z)#AIC: 612
nT1 <- lme(fixed = log.ABUNDANCE ~ 1, random = ~1|SPP, method = "ML", data=TT1z)          #AIC: 705

AIC(nT1cn,nT1c,nT1)
anova(nT1cn,nT1c)
anova(nT1cn,nT1)
summary(nT1cn)


####################################################################################
## CORRELATION STRUCTURE BY SPECIES
####################################################################################

PT1 <- cbind(PT1, dummy=rep(1, 31))
a <- lme(fixed = log(dalo_nL+0.5) ~ 1, random = ~1|dummy,correlation = corExp(form = ~ long + lat, nugget=T), method = "ML", data=PT1,control=lmeControl(opt = "optim"))
b <- lme(fixed = log(dalo_nL+0.5) ~ 1, random = ~1|dummy,correlation = corExp(form = ~ long + lat, nugget=F), method = "ML", data=PT1,control=lmeControl(opt = "optim"))
c <- lme(fixed = log(dalo_nL+0.5) ~ 1, random = ~1|dummy, method = "ML", data=PT1,control=lmeControl(opt = "optim"))          

anova(a,c)
anova(a,b)
summary(a)
summary(b)

# MORAN : parametric test 
library(ade4)
stationT1.dists <- dist(cbind(PT1$long, PT1$lat))

library(ape)
stationT1.dists.inv <- 1/as.matrix(stationT1.dists)
diag(stationT1.dists.inv) <- 0
Moran.I(log(PT1$dalo_nL+0.5), stationT1.dists.inv) 

#Mantel test: semi-parametric
a <- dist(PT1$dalo_nL)
mantel.rtest(stationT1.dists, a, nrepet = 999) 

####################################################################################
## CORRELATION 
####################################################################################

PT1$sampledepth2 <- PT1$sampledepth^2
cT1z <- subset(PT1, select=c("do","k","ph","sampledepth2","exp.wtemp","log.chl_ugl","log.PCYV","logit.Oocystis","logit.Anabaena","logit.Diatoms","logit.OtherChlorophyta","logit.Chrysophyceae","SpCond","dist","horn"))
c1=cor(cT1z, use="complete.obs")


####################################################################################
## POISSON REPLICATES 
####################################################################################

T2=read.csv('T2_replicates.csv', header=T)
T3=read.csv('T3_replicates.csv', header=T)

library(lme4)
library(reshape)

T2m <- melt(T2[,-1], id=(c("sampleid","replicate")))
T2mr <- data.frame(reshape(T2m, direction="wide", idvar=c("sampleid", "variable"), timevar="replicate"))
names(T2mr)<-c("sampleid","taxa","rep1","rep2")

T3m <- melt(T3[,-1], id=(c("sampleid","replicate")))
T3mr <- data.frame(reshape(T3m, direction="wide", idvar=c("sampleid", "variable"), timevar="replicate"))
names(T3mr)<-c("sampleid","taxa","rep1","rep2")

plot(rep1 ~ rep2 , data = T2mr)
summary(m2 <- glm(rep1 ~ rep2 + taxa, family = "poisson", data = T2mr))

plot(rep1 ~ rep2 , data = T3mr)
summary(m2 <- glm(rep1 ~ rep2 + taxa, family = "poisson", data = T3mr))

##################
TT <- T2
TT$total <- rowSums(TT[,4:14])
R2 <- NULL
for(i in 4:15){
  x <- TT[,i]
  V.total <- var(x)
  xx <- (x[TT$replicate==1] + x[TT$replicate==2])/2
  V.site <- (14/15)*(30/29)*var(xx)
  R2 <- rbind(R2,c(mean(x),1-V.site/V.total))
}

TT <- T3
TT$total <- rowSums(TT[,4:14])
R2 <- NULL
for(i in 4:15){
  x <- TT[,i]
  V.total <- var(x)
  xx <- (x[TT$replicate==1] + x[TT$replicate==2])/2
  V.site <- (14/15)*(30/29)*var(xx)
  R2 <- rbind(R2,c(mean(x),1-V.site/V.total))
}

var2 <- T2mr$rep1^2 + T2mr$rep2^2 - 2*((T2mr$rep1+T2mr$rep2)/2)^2
VM <- var2/((T2mr$rep1+T2mr$rep2)/2)
cbind(T2mr,VM)
hist(VM)