---
title: 'D. galeata & Fish Kairomones -part I: GxE-'
author: "Verena Tams"
date: "August 2018"
output: html_document
---
In this first part of the analysis I test for the 'Genotypex Environment interaction (response ~ treatment * clone +(1|round)). 

For hypothesis testing (significant 'Genotype x Environment' interaction) I use the function Anova(model, type=2). This test performs a Wald Chi-Square-Test with the underlying null hypothesis that there is no significant interaction between the two explanatory variables.

```{r SetWd, include=FALSE}
#setwd()
#citation()
#sessionInfo()
```

```{r LoadData, include=FALSE}
FKn15<- read.table("./Input/FKmaster.txt", header = T)
rates<- read.table("./Input/surv_repro_relfit.txt", header = T)
FKn15.no0<- FKn15[ which(FKn15$AFR > 5), ] 
FKn15.no_na<- subset(FKn15, !is.na(length))
FKn15.no0$Indiv <- factor(1:nrow(FKn15.no0)) 
```

```{r libraries, include=FALSE}
library(lme4)
library(ordinal)
library(MASS)
library(car)
library(RVAideMemoire)
```


###Age at First Reproduction (AFR)
####1. Data exploration
```{r PlotAFR1, echo=FALSE, cache=FALSE, fig.show='hide', warning=FALSE, fig.width=6, fig.height=4}
#boxplot:treatment wise
plot(AFR ~ treatment, data = FKn15.no0)
#data distribution
#plot(table(FKn15.no0$AFR)) 
```

```{r PlotAFR2, echo=FALSE, cache=FALSE, fig.show='hide', fig.cap="Plot of Age at First Reproduction", warning=FALSE, fig.width=6, fig.height=4}
#distribution diagnostics
#qqp(FKn15.no0$AFR, "norm") 
#qqp(FKn15.no0$AFR, "lnorm") 
#binom<-fitdistr(FKn15.no0$AFR, "Binomial") 
#qqp(FKn15.no0$AFR, "binom", size = binom$estimate[[1]], mu = binom$estimate[[2]])
#nbinom<-fitdistr(FKn15.no0$AFR, "Negative Binomial") 
#qqp(FKn15.no0$AFR, "nbinom", size = nbinom$estimate[[1]], mu = nbinom$estimate[[2]])
#poisson<-fitdistr(FKn15.no0$AFR, "Poisson") 
#qqp(FKn15.no0$AFR, "pois", poisson$estimate)
#gamma <- fitdistr(FKn15.no0$AFR, "gamma")  
#qqp(FKn15.no0$AFR, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])
```

```{r PlotAFR3, echo=FALSE, cache=FALSE, fig.cap="Plot of Age at First Reproduction, with fitted curves for Poisson and Gamma distributions", warning=FALSE, fig.width=6, fig.height=4, fig.show='hide'}
#data distribution
At <- sort(unique(FKn15.no0$AFR))
Mean.AFR <- mean(FKn15.no0$AFR)
Var.AFR <- var(FKn15.no0$AFR)
shape.AFR <- Mean.AFR^2/Var.AFR
scale.AFR <- Mean.AFR/Var.AFR
Fitted.Pois <- dpois(At, lambda=Mean.AFR)*nrow(FKn15.no0)
Fitted.Gamma <- dgamma(At, shape=shape.AFR, scale=1/scale.AFR)*nrow(FKn15.no0)
Tab.AFR <- table(FKn15.no0$AFR)
par(mfrow=c(1,1), mar=c(4,4,1,1))
plot(Tab.AFR, yaxs="i", ylim=c(0,1.04*max(Tab.AFR)), lwd=20, lend=4, xlab="AFR", ylab="Frequency")
lines(At, Fitted.Pois, col=3, lwd=2)
lines(At, Fitted.Gamma, col=4, lwd=2)
legend(10.5, 220, c("Poisson", "Gamma"), lty=1, lwd=2, col=c(3,4))
```

####2. Test of Normality/Test for gaussian distribution
Shapiro-Wilk-Test [shapiro()]: The null hypothesis of normally distributed data is rejected, since the p value is < 2.2e-16. 
The qqplot confirms a non-normal distribution.
```{r TestAFR, echo=FALSE, cache=FALSE, results='hide', plot='hide', warning=FALSE}
##############2. Test of Normality ###############################
##formal hypothesis testing in R: Shapiro-Wilk Test
#H(0)= data is normally distributed => null hypothesis
#H(A)= data is not normally distributed => alternative hypothesis

#reject H(0), if p<0.05 (less) => data is NOT normally distributed
#fail to reject H(0), if p>0.05(greater) => can assume normality
#qqplot() normal probability plot
#qqline() adds a line which passes through the first and third quartiles

AFR<- FKn15.no0[c(3:5, 10)] # select 3rd through 5th and 10th variables (column)
AFR_m<- data.matrix(AFR)
#qqnorm(AFR_m); qqline(AFR_m, col = 2, lwd=2, lty=2)
summary(AFR_m)
shapiro.test(AFR_m) #W = 0.82707, p-value < 2.2e-16 => reject H0
```
=> AFR: gamma distribution 

####3. Model application
```{r ModelAFR}
#random intercept model, single fixed factor
AFR.model <- glm(AFR ~ treatment*clone, data=FKn15.no0)
Anova(AFR.model, type=2) #***

#random intercept model plus random factor 'round'
AFR.model.gamma <- glmer(AFR ~ treatment *clone +(1|round), family=Gamma(link="log"), data=FKn15.no0)
AFR.null.gamma <- glmer(AFR ~ treatment+clone +(1|round), family=Gamma(link="log"), data=FKn15.no0)
anova(AFR.null.gamma,AFR.model.gamma) 
anova(AFR.model,AFR.model.gamma) #round? #?
Anova(AFR.model.gamma, type=2) #GxE?***, G***, E***
```


###Total Number of Broods (broods)
####1. Data exploration
```{r PlotBroods1, echo=FALSE, cache=FALSE, warning=FALSE, fig.width=6, fig.height=4, fig.show='hide'}
#boxplot:treatment wise
#plot(broods ~ treatment, data = FKn15.no0)
#data distribution
plot(table(FKn15.no0$broods)) 
```

```{r PlotBroods2, echo=FALSE, cache=FALSE, fig.cap="Plot of Total Number of Broods", warning=FALSE, fig.width=6, fig.height=4, fig.show='hide'}
#distribution diagnostics
#qqp(FKn15.no0$broods, "norm")
#qqp(FKn15.no0$broods, "lnorm")
#qqp(FKn15.no0$broods, "ordinal")
#binom<-fitdistr(FKn15.no0$broods, "Binomial")
#qqp(FKn15.no0$broods, "binom", size = binom$estimate[[1]], mu = binom$estimate[[2]])
#nbinom<-fitdistr(FKn15.no0$broods, "Negative Binomial")
#qqp(FKn15.no0$broods, "nbinom", size = nbinom$estimate[[1]], mu = nbinom$estimate[[2]])
#poisson<-fitdistr(FKn15.no0$broods, "Poisson")
#qqp(FKn15.no0$broods, "pois", poisson$estimate)
#gamma <- fitdistr(FKn15.no0$broods, "gamma")
#qqp(FKn15.no0$broods, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])
#ordinal<-fitdistr(FKn15.no0$broods, "ordinal")#error: nicht unterstütze Verteilung
#qqp(FKn15.no0$broods, "ordinal", ordinal$esimate)
```

####2. Test of Normality/Test for gaussian distribution
Shapiro-Wilk-Test [shapiro()]: The null hypothesis of normally distributed data is rejected, since the p value is < 2.2e-16. 
The qqplot confirms a non-normal distribution.

```{r TestBroods, echo=FALSE, cache=FALSE, results='hide', fig.show='hide', warning=FALSE}
broods<- FKn15.no0[c(3:5, 11)]
broods_m<- data.matrix(broods)
#hist(broods_m) # skewed to the left 
#qqnorm(broods_m); qqline(broods_m, col = 2, lwd=2, lty=2)
shapiro.test(broods_m) #W = 0.63, p-value < 2.2e-16 => reject H0 
```
=> broods: ordinal distribution

####3. Model application
```{r ModelBroods}
# fit the model
FKn15.no0$broodsF <- factor(FKn15.no0$broods, levels = sort(unique(FKn15.no0$broods)), ordered=TRUE)

#random intercept model plus random factor 'round'
broods.model.clmm <- clmm(broodsF ~ treatment *clone +(1|round), data=FKn15.no0, link="probit")
broods.null.clmm <- clmm(broodsF ~ treatment+clone + (1|round), data=FKn15.no0, link="probit")
anova(broods.null.clmm, broods.model.clmm) # ***
Anova.clmm(broods.model.clmm) #GxE?
```

###Total Number of Offspring (offspring)
####1. Data exploration
```{r PlotOff1, echo=FALSE, cache=FALSE, warning=FALSE, fig.show='hide'}
#boxplot:treatment wise
#plot(offspring ~ treatment, data = FKn15.no0)
#data distribution
plot(table(FKn15.no0$offspring )) 
```

```{r PlotOff2, echo=FALSE, cache=FALSE, fig.cap="Plot of Total Number of Offspring", warning=FALSE, fig.show='hide'}
#distribution diagnostics
#qqp(FKn15.no0$offspring, "norm")
#qqp(FKn15.no0$offspring, "lnorm")
#binom<-fitdistr(FKn15.no0$offspring, "Binomial")
#qqp(FKn15.no0$broods, "binom", size = binom$estimate[[1]], mu = binom$estimate[[2]])
#nbinom<-fitdistr(FKn15.no0$offspring, "Negative Binomial")
#qqp(FKn15.no0$offspring, "nbinom", size = nbinom$estimate[[1]], mu = nbinom$estimate[[2]])
poisson<-fitdistr(FKn15.no0$offspring, "Poisson") 
qqp(FKn15.no0$offspring, "pois", poisson$estimate)
#gamma <- fitdistr(FKn15.no0$offspring, "gamma")
#qqp(FKn15.no0$offspring, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])
```

####2. Test of Normality/Test for gaussian distribution
Shapiro-Wilk-Test [shapiro()]: The null hypothesis of normally distributed data is rejected, since the p value is < 2.2e-16. 
The qqplot confirms a non-normal distribution.
```{r TestOff, echo=FALSE, cache=FALSE, results='hide', fig.show='hide', warning=FALSE}
offspring<- FKn15.no0[c(3:5, 12)]
offspring_m<- data.matrix(offspring)
#hist(offspring_m) # skewed to the left zero inflation=> subset
#qqnorm(offspring_m);qqline(offspring_m, col = 2, lwd=2, lty=2)
shapiro.test(offspring_m) #W = 0.81643, p-value < 2.2e-16 => reject H0 
```
=> offspring: poisson distribution plus overdsispersion term

####3. Model application
```{r ModelOff}
#random intercept model plus random factor 'round'
offspring.model.glmm <- glmer(offspring ~ treatment *clone +(1|round), data=FKn15.no0, family = poisson())
offspring.null.glmm <- glmer(offspring ~ treatment + clone +(1|round), data=FKn15.no0, family = poisson())
anova(offspring.null.glmm,offspring.model.glmm) 
Anova(offspring.model.glmm, type=2) #GxE?***, G***, E***
```


###Total Number of Offspring of First Brood (brood1)
####1. Data exploration
```{r PlotBrood1.1, echo=FALSE, cache=FALSE, warning=FALSE, fig.show='hide'}
#boxplot:treatment wise
#plot(brood1 ~ treatment, data = FKn15.no0)
#data distribution
plot(table(FKn15.no0$brood1)) 
```

```{r PlotBrood1.2, echo=FALSE, cache=FALSE, fig.cap="Plot of Total Number of Offspring First Brood", warning=FALSE, fig.show='hide'}
#distribution diagnostics
#qqp(FKn15.no0$brood1, "norm")
#qqp(FKn15.no0$brood1, "lnorm")
#nbinom<-fitdistr(FKn15.no0$brood1, "Negative Binomial")
#qqp(FKn15.no0$brood1, "nbinom", size = nbinom$estimate[[1]], mu = nbinom$estimate[[2]])
#poisson<-fitdistr(FKn15.no0$brood1, "Poisson")
#qqp(FKn15.no0$brood1, "pois", poisson$estimate)
#gamma <- fitdistr(FKn15.no0$brood1, "gamma")
#qqp(FKn15.no0$brood1, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])
```

####2. Test of Normality/Test for gaussian distribution
Shapiro-Wilk-Test [shapiro()]: The null hypothesis of normally distributed data is rejected, since the p value is < 2.2e-16. 
The qqplot confirms a non-normal distribution.
```{r TestBrood1, echo=FALSE, cache=FALSE, results='hide', fig.show='hide', warning=FALSE}
brood1<- FKn15.no0[c(3:5, 6)]
brood1_m<- data.matrix(brood1)
#hist(brood1_m) # skewed to the left
#qqnorm(brood1_m); qqline(brood1_m, col = 2, lwd=2, lty=2)
shapiro.test(brood1_m) #W = 0.69999, p-value < 2.2e-16 => reject H0
```
=> brood1: poisson distribution plus overdispersion term

####3. Model application
```{r ModelBrood1}
#random intercept model plus random factor 'round'
brood1.model.glmm <- glmer(brood1 ~ treatment *clone + (1|round), data=FKn15.no0, family = poisson())
brood1.null.glmm <- glmer(brood1 ~ treatment +clone+ (1|round), data=FKn15.no0, family = poisson())
anova(brood1.null.glmm, brood1.model.glmm)
Anova(brood1.model.glmm, type=2) #G***
```


###Survival (surv)
####1. Data exploration
```{r PlotSurvived1, echo=FALSE, cache=FALSE, warning=FALSE, fig.width=6, fig.height=4, fig.show='hide'}
#boxplot:treatment wise
#plot(survived_t14 ~ treatment, data = FKn15)
#data distribution
#plot(table(FKn15$survived_t14))
```

```{r PlotSurvived2, echo=FALSE, cache=FALSE, fig.cap="Plot of Survival", fig.width=6, fig.height=4, warning=FALSE, fig.show='hide'}
#distribution diagnostics
#qqp(FKn15$survived_t14, "norm")
#qqp(FKn15$survived_t14, "lnorm")
#nbinom<-fitdistr(FKn150$survived_t14, "Negative Binomial")
#qqp(FKn15$survived_t14, "nbinom", size = nbinom$estimate[[1]], mu = nbinom$estimate[[2]])
#poisson<-fitdistr(FKn15$survived_t14, "Poisson")
#qqp(FKn15$survived_t14, "pois", poisson$estimate)
#gamma <- fitdistr(FKn15$survived_t14, "gamma")
#qqp(FKn15$survived_t14, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])
```

####2. Test of Normality/Test for gaussian distribution
Shapiro-Wilk-Test [shapiro()]: The null hypothesis of normally distributed data is rejected, since the p value is < 2.2e-16. 
The qqplot confirms a non-normal distribution.
```{r TestSurvived, echo=FALSE, cache=FALSE, results='hide', fig.show='hide', warning=FALSE}
surv<- FKn15[c(3:5, 13)]
surv_m<- data.matrix(surv)
#hist(surv_m) # skewed to the left 
#qqnorm(surv_m); qqline(surv_m, col = 2, lwd=2, lty=2)
shapiro.test(surv_m) #W = 0.61741, p-value < 2.2e-16 => reject H0 
```
=> survived_t14: binomial distribution

####3. Model application
```{r ModelSurv}
#random intercept model plus random factor 'round'
surv.model <- glm(survived_t14 ~ treatment *clone + (1|round), data=FKn15, family = binomial(link = "logit"))
surv.null <- glm(survived_t14 ~ treatment + clone + (1|round), data=FKn15, family = binomial(link = "logit"))
anova(surv.null,surv.model) 
Anova(surv.model, type=2) #n.s., n.s., n.s.
```


###relative fitness within each population (relnest)
####1. Data exploration
```{r PlotRel.nest1, echo=FALSE, cache=FALSE,fig.show='hide', warning=FALSE, fig.show='hide'}
#boxplot:treatment wise
plot(relfit.nested ~ treatment, data = rates)
#data distribution
plot(table(rates$relfit.nested)) 
```

```{r PlotRel.nest2, echo=FALSE, cache=FALSE, fig.cap="Relative Fitness within each population", warning=FALSE, fig.show='hide'}
#distribution diagnostics
qqp(rates$relfit.nested, "norm")
#qqp(rates$relfit.nested, "lnorm")
#nbinom<-fitdistr(rates$relfit.nested, "Negative Binomial")
#qqp(rates$relfit.nested, "nbinom", size = nbinom$estimate[[1]], mu = nbinom$estimate[[2]])
#poisson<-fitdistr(rates$relfit.nested, "Poisson")
#qqp(rates$relfit.nested, "pois", poisson$estimate)
gamma <- fitdistr(rates$relfit.nested, "gamma")
qqp(rates$relfit.nested, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])
```

####2. Test of Normality/Test for gaussian distribution
Shapiro-Wilk-Test [shapiro()]: The null hypothesis of normally distributed data is rejected, since the p value is < 2.2e-16. 
The qqplot confirms a non-normal distribution.
```{r TestRel.nest, echo=FALSE, cache=FALSE, results='hide', fig.show='hide', warning=FALSE}
rel.nest<- rates[c(2:4, 7)]
rel.nest_m<- data.matrix(rel.nest)
#hist(rel.nest_m) # skewed to the left
#qqnorm(rel.nest_m); qqline(rel.nest_m, col = 2, lwd=2, lty=2)
shapiro.test(rel.nest_m) #W = 0.64043, p-value < 2.2e-16 => reject H0 
```
=> relfit.nested: gamma distribution

####3. Model application
```{r ModelRelnest}
#random intercept model plus random factor 'round'
#relnest.model.glm <- glmer(relfit.nested ~ treatment * clone + (1|round), data=rates, family = Gamma(link="log")) #error
#relnest.model.glm <- glm(relfit.nested ~ treatment * clone + (1|round), data=rates, family = Gamma(link="log")) #NaNs
#relnest.model.glm <- glm(relfit.nested ~ treatment * clone, data=rates, family = Gamma(link="log")) #NaNs

relnest.null.glm <- glmer(relfit.nested ~ treatment + clone + (1|round), data=rates, family = Gamma(link="log"))
#anova(relnest.null.glm, relnest.model.glm) 
Anova(relnest.null.glm, type=2) 
```
  

###relative fitness among all populations (relclone)
####1. Data exploration
```{r PlotRel.clone1, echo=FALSE, cache=FALSE, fig.show='hide', warning=FALSE}
#boxplot:treatment wise
plot(relfit.clone ~ treatment, data = rates)
#data distribution
plot(table(rates$relfit.clone)) 
```

```{r PlotRel.clone2, echo=FALSE, cache=FALSE, fig.cap="Relative Fitness among populations", fig.show='hide', warning=FALSE}
#distribution diagnostics
qqp(rates$relfit.clone, "norm")
#qqp(rates$relfit.clone, "lnorm")
#nbinom<-fitdistr(rates$relfit.clone, "Negative Binomial")
#qqp(rates$relfit.clone, "nbinom", size = nbinom$estimate[[1]], mu = nbinom$estimate[[2]])
#poisson<-fitdistr(rates$relfit.clone, "Poisson")
#qqp(rates$relfit.clone, "pois", poisson$estimate)
gamma <- fitdistr(rates$relfit.clone, "gamma")
qqp(rates$relfit.clone, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])
```

####2. Test of Normality/Test for gaussian distribution
Shapiro-Wilk-Test [shapiro()]: The null hypothesis of normally distributed data is rejected, since the p value is < 2.2e-16. 
The qqplot confirms a non-normal distribution.
```{r TestRel.clone, echo=FALSE, cache=FALSE, results='hide', fig.show='hide', warning=FALSE}
rel.clone<- rates[c(2:4, 8)]
rel.clone_m<- data.matrix(rel.clone)
#hist(rel.clone_m) # skewed to the left zero inflation?
#qqnorm(rel.clone_m); qqline(rel.clone_m, col = 2, lwd=2, lty=2)
shapiro.test(rel.clone_m) #W = 0.6477, p-value < 2.2e-16 => reject H0 
```
=> relfit.nested: gamma distribution

####3. Model application
```{r ModelRelclone}
#random intercept model plus random factor 'round'
#relclone.model <- glmer(relfit.clone ~ treatment * clone + (1|round), data=rates, family = Gamma(link="log")) #error
relclone.null <- glmer(relfit.clone ~ treatment + clone + (1|round), data=rates, family = Gamma(link="log"))
#anova(relclone.null, relclone.model)
Anova(relclone.null, type=2) 
```
    

###Somatic Growth rate (SGR)
####1. Data exploration
```{r PlotSGR1, echo=FALSE, cache=FALSE, fig.show='hide', warning=FALSE}
#boxplot:treatment wise
plot(SGR ~ treatment, data = FKn15)
#data distribution 
plot(table(FKn15$SGR)) 
```

```{r PlotSGR2, echo=FALSE, cache=FALSE, fig.cap="Somatic Growth Rate (µm/day)", fig.show='hide', warning=FALSE}
#distribution diagnostics
qqp(FKn15$SGR,"norm")
#qqp(FKn15$SGR, "lnorm")
#nbinom<-fitdistr(FKn15$SGR, "Negative Binomial")
#qqp(FKn15$SGR, "nbinom", size = nbinom$estimate[[1]], mu = nbinom$estimate[[2]])
#poisson<-fitdistr(FKn15$SGR, "Poisson") 
#qqp(FKn15$SGR, "pois", poisson$estimate)
#gamma <- fitdistr(FKn15$SGR, "gamma")
#qqp(FKn15$SGR, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])
```

```{r PlotSGR3, echo=FALSE, cache=FALSE, fig.cap="Somatic Growth Rate (µm/day), with fitted curves for Poisson and Gamma distributions", fig.width=6, fig.height=4, warning=FALSE, fig.show='hide'}
#data distribution
At <- sort(unique(FKn15$SGR))
Mean.SGR <- mean(FKn15$SGR)
Var.SGR <- var(FKn15$SGR)
shape.SGR <- Mean.SGR^2/Var.SGR
scale.SGR<- Mean.SGR/Var.SGR
Fitted.Pois <- dpois(At, lambda=Mean.SGR)*nrow(FKn15)
Fitted.Gamma <- dgamma(At, shape=shape.SGR, scale=1/scale.SGR)*nrow(FKn15)
Tab.SGR<- table(FKn15$SGR)
par(mfrow=c(1,1), mar=c(4,4,1,1))
plot(Tab.SGR, yaxs="i", ylim=c(0,1.04*max(Tab.SGR)), lwd=20, lend=4, xlab="SGR", ylab="Frequency")
lines(At, Fitted.Pois, col=3, lwd=2)
lines(At, Fitted.Gamma, col=4, lwd=2)
legend(10.5, 220, c("Poisson", "Gamma"), lty=1, lwd=2, col=c(3,4))
```
####2. Test of Normality/Test for gaussian distribution
Shapiro-Wilk-Test [shapiro()]: The null hypothesis of normally distributed data is rejected, since the p value is < 2.2e-16. 
The qqplot confirms a non-normal distribution.
```{r TestSGR, echo=FALSE, cache=FALSE, results='hide', fig.show='hide', warning=FALSE}
SGR<- FKn15[c(3:5, 14)]
SGR_m<- data.matrix(SGR)
#hist(SGR_m) # skewed to the left
#qqnorm(SGR_m); qqline(SGR_m, col = 2, lwd=2, lty=2)
shapiro.test(SGR_m) #W = 0.69504, p-value < 2.2e-16 => reject H0 
```
=> SGR: gaussian distribution

####3. Model application
```{r ModelSGR}
#random intercept model plus random factor 'round'
SGR.model.glmm <- lmer(SGR ~ treatment *clone +(1|round), data=FKn15, REML=FALSE)
SGR.null.glmm <- lmer(SGR ~ treatment + clone +(1|round), data=FKn15, REML=FALSE)
anova(SGR.null.glmm, SGR.model.glmm) 
Anova(SGR.model.glmm, type=2) #GxE***, G***, E***
```

###Body Length (size)
####1. Data exploration
```{r PlotSize1, echo=FALSE, cache=FALSE, fig.show='hide', warning=FALSE}
#boxplot:treatment wise
plot(length ~ treatment, data = FKn15.no_na)
#data distribution
plot(table(FKn15.no_na$length)) 
```

```{r PlotSize2, echo=FALSE, cache=FALSE, fig.cap="Somatic Growth Rate (µm/day)", fig.show='hide', warning=FALSE}
#distribution diagnostics
#qqp(FKn15.no_na$length,"norm") #=)
#qqp(FKn15.no_na$length, "lnorm")
#nbinom<-fitdistr(FKn15.no_na$length, "Negative Binomial")
#qqp(FKn15.no_na$length, "nbinom", size = nbinom$estimate[[1]], mu = nbinom$estimate[[2]])
#poisson<-fitdistr(FKn15.no_na$length, "Poisson") #=)
#qqp(FKn15.no_na$length, "pois", poisson$estimate)
#gamma <- fitdistr(FKn15.no_na$length, "gamma")
#qqp(FKn15.no_na$length, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])
```

```{r PlotSize3, echo=FALSE, cache=FALSE, fig.cap="Somatic Growth Rate (µm/day), with fitted curves for Poisson and Gamma distributions", fig.width=6, fig.height=4, warning=FALSE, fig.show='hide'}
#data distribution
At <- sort(unique(FKn15.no_na$length))
Mean.size <- mean(FKn15.no_na$length)
Var.size <- var(FKn15.no_na$length)
shape.size <- Mean.size^2/Var.size
scale.size<- Mean.size/Var.size
Fitted.Pois <- dpois(At, lambda=Mean.size)*nrow(FKn15.no_na)
Fitted.Gamma <- dgamma(At, shape=shape.size, scale=1/scale.size)*nrow(FKn15.no_na)
Tab.size<- table(FKn15.no_na$length)
par(mfrow=c(1,1), mar=c(4,4,1,1))
plot(Tab.size, yaxs="i", ylim=c(0,1.04*max(Tab.size)), lwd=20, lend=4, xlab="Size", ylab="Frequency")
lines(At, Fitted.Pois, col=3, lwd=2)
lines(At, Fitted.Gamma, col=4, lwd=2)
legend(10.5, 220, c("Poisson", "Gamma"), lty=1, lwd=2, col=c(3,4))
```

####2. Test of Normality/Test for gaussian distribution

Shapiro-Wilk-Test [shapiro()]: The null hypothesis of normally distributed data is rejected, since the p value is < 2.2e-16. 
The qqplot confirms a non-normal distribution.
```{r TestSize, echo=FALSE, cache=FALSE, results='hide', fig.show='hide', warning=FALSE}
Size<- FKn15.no_na[c(3:5, 15)]
Size_m<- data.matrix(Size)
#hist(Size_m) # skewed to the left
#qqnorm(Size_m); qqline(Size_m, col = 2, lwd=2, lty=2)
shapiro.test(Size_m) #W = 0.57806, p-value < 2.2e-16 => reject H0 
```
=> Size: gamma distribution
####3. Model application
```{r ModelSize}
#random intercept model plus random factor 'round'
length.model.glmm <- glmer(length ~ treatment *clone +(treatment|round), data=FKn15.no_na, family = Gamma(link="log"))
length.null.glmm <- glmer(length ~ treatment + clone +(treatment|round), data=FKn15.no_na, family = Gamma(link="log"))
anova(length.null.glmm, length.model.glmm) 
Anova(length.model.glmm, type=2) #GxE***, G***, E***
```