setwd("D:/Work/Travail_Postdoc/Manuscript d'article/Gojelly/Catostylus tagi/Life cycle draft/Figures and data raw/Submission")


#### Packages ####
# Packages stat
library(car)
library(lme4)
library(Matrix)
library(nlme)
library(MASS)
library(pscl)
library(lmtest)
library(AER)       # to calcultate overdispertion in GLM

# graphes
library(ggplot2)
library(ggpubr)


###################################
###  1. Planula plankton time  ####
###################################

# Data #
{
Planula = read.csv("Planula plankton time.csv", header = TRUE)

Planula$Temperature <- as.factor(Planula$Temperature)
Planula$Salinity <- as.factor(Planula$Salinity)
}

### 1.1. two-ways ANOVA non transformed data ####  
# -> not validated because of the non homogeneity of the variance  
{
# ANOVA
aov1 <- aov(Planktonic ~ Temperature * Salinity, data = Planula)
  
# Normal distribution test
aov_residuals1 <- residuals(object = aov1)
hist(aov_residuals1)
shapiro.test(x = aov_residuals1)

# Variance Homogeneity test
leveneTest(Planktonic ~ Temperature * Salinity, data = Planula)
}

### 1.2. two-ways ANOVA on sqrt transformed data ####  
# ->  Validated

{
# ANOVA
aov2 <- aov(PlankSQRT ~ Temperature * Salinity, data = Planula)
  

# Normal distribution test
aov_residuals2 <- residuals(object = aov2)
hist(aov_residuals2)
shapiro.test(x = aov_residuals2)

# Variance Homogeneity test
leveneTest(PlankSQRT ~ Temperature * Salinity, data = Planula)

# ANOVA output
summary(aov2)

# TukeyHSD test
TukeyHSD(aov2)
tukey.plot.aov<-aov(PlankSQRT ~ Temperature : Salinity, data=Planula)
tukey.plot.test<-TukeyHSD(tukey.plot.aov)
plot(tukey.plot.test, las = 1)
}

##############################
###  2. Pre-strobilation  ####
##############################

# Data #
{
Strobilation = read.csv("Strobilation_Prestr.csv", header = TRUE)

Strobilation$Temperature <- as.factor(Strobilation$Temperature)
Strobilation$Feed <- as.factor(Strobilation$Feed)

#plot the data
ggplot(Strobilation, aes(x=Prestr)) +  geom_histogram() # -> poisson distribution
}

### 2.1. two-ways ANOVA on non transformed data ####
{
# ANOVA
aov3 <- aov(Prestr ~ Temperature * Feed, data = Strobilation)

# Normal distribution test
aov_residuals3 <- residuals(object = aov3)
hist(aov_residuals3)
shapiro.test(x = aov_residuals3)

# Variance Homogeneity test
leveneTest(Prestr ~ Temperature * Feed, data = Strobilation)
}

### 2.2. GLM Poisson distribution ####
# -> the model ended with an overdispertion. So we moved to a GLM negative binomial distribution
{
Strobilation$Temperature <- as.numeric(Strobilation$Temperature)
Strobilation$Feed <- as.numeric(Strobilation$Feed)

# GLM model
fit.p <- glm(Prestr ~ Temperature + Feed, data=Strobilation, family=poisson())
summary(fit.p)

# Check for over/underdispersion in the model
dispersiontest(fit.p,trafo=1)
}

### 2.3. GLM Negative binomial distribution ####
# -> the model ended with an underdispersion. 
{
# GLM model
fit.nb <- glm.nb(Prestr ~ Temperature + Feed, data=Strobilation)
summary(fit.nb)

# Check for over/underdispersion in the model
E1 <- resid(fit.nb, type = "pearson")
N1  <- nrow(Strobilation)
p1  <- length(coef(fit.nb)) + 1  
sum(E1^2) / (N1 - p1)  
}

# --> Both GLMS showed a high overdispertion with migth be explained by the high proportion of zero: 72.2% of the data are zeros
100*sum(Strobilation$Prestr == 0)/nrow(Strobilation) 
# Therefore, we moved to Zero-inlfated models


### 2.3. Zero-inflated Poisson distribution ####
  
{
# Model
fit.zp <- zeroinfl(Prestr ~ Temperature + Feed |Temperature + Feed, dist = 'poisson', data=Strobilation)
summary (fit.zp)

# Dispersion statistic

E2 <- resid(fit.zp, type = "pearson")
N2  <- nrow(Strobilation)
p2  <- length(coef(fit.zp))  
sum(E2^2) / (N2 - p2)
}

### 2.4. Zero-inflated Negative Binomial distribution ####
  
{
# Model
fit.znb <- zeroinfl(Prestr ~ Temperature + Feed |Temperature + Feed, dist = 'negbin', data=Strobilation)
summary (fit.znb)

# Dispersion statistic

E3 <- resid(fit.znb, type = "pearson")
N3  <- nrow(Strobilation)
p3  <- length(coef(fit.znb)) + 1  
sum(E2^2) / (N3 - p3)
}

# compare between ZINP,ZINB,GML Poisson & GLM Negative Binomial
{
lrtest(fit.p, fit.nb, fit.zp, fit.znb)
vuong(fit.zp,fit.znb)

}

################################################################
###  3. Podocysts: Temperature & feed #####
################################################################

# Data 
{
Podocysts = read.csv("Podocyst.csv", header = TRUE)
  
Podocysts$Temperature <- as.factor(Podocysts$Temperature)
Podocysts$Feed <- as.factor(Podocysts$Feed)
}

### 3.1. two-ways ANOVA non transformed data ####  
# ->   
{
# ANOVA
aov4 = aov(Podocyst ~ Temperature * Feed, data = Podocysts)
  
# Normal distribution test
aov_residuals4 <- residuals(object = aov4)
hist(aov_residuals4)
shapiro.test(x = aov_residuals4)
  
# Variance Homogeneity test
leveneTest(Podocyst ~ Temperature * Feed, data = Podocysts)
  
# ANOVA output
summary(aov4)
}

###############################
###  4. Strobilation_bet  #####
###############################
#-> no significant difference
{
# Data 
Strob.bet = read.csv("Strobilation_bet&str.csv", header = TRUE)
  
## Shapiro-Wilk normality test 
with(Strob.bet, tapply(betstr, Feed, shapiro.test))

  
# F-test to test for homogeneity in variances
bet.ftest <- var.test(betstr ~ Feed, data = Strob.bet)
bet.ftest
  
# wilcox test
bet.test <- wilcox.test(betstr ~ Feed, data = Strob.bet, exact =FALSE)
bet.test
}

###############################
###  5. Strobilation_str  #####
###############################
#-> significant difference
{
## Shapiro-Wilk normality test 
with(Strob.bet, tapply(str, Feed, shapiro.test))

# F-test to test for homogeneity in variances
str.ftest <- var.test(str ~ Feed, data = Strob.bet)
str.ftest
  
# t-test
str.test <- t.test(str ~ Feed, data = Strob.bet, alternative = "two.sided", var.equal = TRUE)
str.test
}

###############################
###  6. Ephyrae released  #####
###############################
#-> no significant difference
{
# Data 
Ephyrae = read.csv("Ephyrae.csv", header = TRUE)

## Shapiro-Wilk normality test 
with(Ephyrae, tapply(Ephyra, Feed, shapiro.test))

# F-test to test for homogeneity in variances
eph.ftest <- var.test(Ephyra ~ Feed, data = Ephyrae)
eph.ftest

# t-test
eph.test <- t.test(Ephyra ~ Feed, data = Ephyrae, alternative = "two.sided", var.equal = TRUE)
eph.test
}