#### Nicolakis et al. 2025 - BMC Biology - "Pre- and post-copulatory traits are affected by experimental inbreeding, but they are not correlated" ####
#### -------------------------------------------------------------------------------------------------------------------------------------------- ####


########################################################################################################################################
#---------------------------------------------------------------------------------------------------------------------------------------
#### R-code to calculate ANOSIM and crate nMDS graphs to compare Syllable type usage between breeding treatments (O vs I vs II) or inbred groups (Inbred vs Outbred),
#### separately for Introduction phase (using 15 and 10 syllable types) and Interaction phase (using 10 syllable types)
#---------------------------------------------------------------------------------------------------------------------------------------
########################################################################################################################################


######### load data #####################

# load data for introduction using 15 syllable types
library(readxl)
wall_USV_syllabletypes <- read_excel("C:/Documents/stats/wall_USV syllabletypes.xlsx")
View(wall_USV_syllabletypes)

# load data for introduction using 10 syllable classes (pooled simple types to make it comparable with interaction)
library(readxl)
wall_USV_syllabletypes_pooled <- read_excel("C:/Documents/stats/wall_USV syllabletypes_pooled.xlsx")
View(wall_USV_syllabletypes_pooled)

# load data for interaction using 10 syllable types
library(readxl)
contact_USV_syllabletypes <- read_excel("C:/Documents/stats/contact_USV syllabletypes.xlsx")
View(contact_USV_syllabletypes)

# load data with grouping variables
library(readxl)
groups <- read_excel("C:/Documents/stats/groups.xlsx")
View(groups)


library(vegan)

#rename datasets
wall <- wall_USV_syllabletypes
wall.pooled <- wall_USV_syllabletypes_pooled
contact <- contact_USV_syllabletypes
groups <- groups

# load vegan package
library("vegan", lib.loc="~/R/win-library/3.5")
#citation("vegan")

#################################################
### INTRODUCTION Phase - Using 15 syllable types #################
##################################################


############## STATS - Anosim #################

##### ----- comparing 3 breeding treatments (OO vs OI vs II)----- ######

#anosim for wall for breeding treatment
wall.dist <- vegdist(wall)
wall.ano <- anosim(wall.dist, groups$Genetic.Status_OO.1_OI.2_II.3)
summary(wall.ano)
plot(wall.ano)


#to test which of the syllable types contributes the most to the observed differences --> which species are responsible for the difference between groups
sim<-simper(wall, group=groups$Genetic.Status_OO.1_OI.2_II.3)
sim
summary(sim)


####################  GRAPHS: nMDS-graphs  #########################

##### ----- comparing 3 breeding treatments (OO-OI-II)----- ######


#nMDS for wall on USVs vs breeding treatment
wall.nMDS<-metaMDS(wall, distance = "bray")


# Create a variable to assign symbols and colors in the graph
symbols <- groups$Genetic.Status_OO.1_OI.2_II.3
symbols
symbols[symbols==1]<-16 # filled circle
symbols[symbols==2]<-17 # filled triangle
symbols[symbols==3]<-17 # filled triangle
symbols


color <- groups$Genetic.Status_OO.1_OI.2_II.3
color
color[color==1]<-"deepskyblue3" # t?rkis, blue
color[color==2]<-"firebrick1" # middle red
color[color==3]<-"firebrick4" # dark red
color

# for the legend
inbreed.status <- c("Outbred","Inbred 1.gen", "Inbred 2.gen")
legendcolor <- c("deepskyblue3", "firebrick1", "firebrick4")

fig<-plot(wall.nMDS, type = "n")
points(fig, "sites", cex=c(1.5),pch = symbols,col = color)
text(fig, "species", font=c(2), cex=c(1), col = "black")
legend("bottomright", legend = inbreed.status, cex = c(1.3), col = legendcolor, pch = c(16,17))
#ordispider(fig, groups=groups$Genetic.Status_OO.1_OI.2_II.3, col = legendcolor, lty="longdash") # lines connect to the centroid
ordihull(fig, groups=groups$Genetic.Status_OO.1_OI.2_II.3, col = legendcolor, lty="dotted") # circle around each group
title("Introduction")


############################################################################

############## STATS - Anosim #################

##### -----Pooling Inbred groups (Outbred vs inbred)----- ######

# anosim for wall for breeding treatment
wall.dist <- vegdist(wall)
wall.ano <- anosim(wall.dist, groups$PooledGenetic_OO.1_OI.II.2)
summary(wall.ano)
plot(wall.ano)


#to test which of the syllable types contributes the most to the observed differences --> which species are responsible for the difference between groups
sim<-simper(wall, group=groups$PooledGenetic_OO.1_OI.II.2)
sim
summary(sim)


####################  GRAPHS: nMDS-graphs  ##########################

##### ----- comparing 2 groups (OB vs IB)----- ######


# nMDS for wall on USVs vs breeding groups
wall.nMDS<-metaMDS(wall, distance = "bray")


# Create a variable to assign symbols and colors in the graph
symbols <- groups$PooledGenetic_OO.1_OI.II.2
symbols
symbols[symbols==1]<-16 # filled circle
symbols[symbols==2]<-17 # filled triangle
symbols


color <- groups$PooledGenetic_OO.1_OI.II.2
color
color[color==1]<-"deepskyblue3" # t?rkis, blue
color[color==2]<-"firebrick1" # middle red
color

#for the legend
inbreed.group <- c("Outbred","Inbred")
legendcolor <- c("deepskyblue3", "firebrick1")

fig<-plot(wall.nMDS, type = "n")
points(fig, "sites", cex=c(1.5),pch = symbols,col = color)
text(fig, "species", font=c(2), cex=c(1), col = "black")
legend("bottomright", legend = inbreed.group, cex = c(1.3), col = legendcolor, pch = c(16,17))
#ordispider(fig, groups=groups$PooledGenetic_OO.1_OI.II.2, col = legendcolor, lty="longdash") # lines connect to the centroid
ordihull(fig, groups=groups$PooledGenetic_OO.1_OI.II.2, col = legendcolor, lty="dotted") # circle around each group
title("Introduction")





#################################################
### INTRODUCTION Phase - Using 10 syllable types (pooled simple types to compare with contact 10 types) #################
##################################################



############## STATS - Anosim  #################

##### ----- comparing 3 breeding treatments (OO vs OI vs II)----- ######

# anosim for wall for breeding treatments
wallpooled.dist <- vegdist(wall.pooled)
wallpooled.ano <- anosim(wallpooled.dist, groups$Genetic.Status_OO.1_OI.2_II.3)
summary(wallpooled.ano)
plot(wallpooled.ano)


# test which of the syllable types contributes the most to the observed differences --> which species are responsible for the difference between groups
sim<-simper(wall.pooled, group=groups$Genetic.Status_OO.1_OI.2_II.3)
sim
summary(sim)


####################  GRAPHS: nMDS-graphs  ##########################

##### ----- comparing 3 breeding treatments (OO-OI-II)----- ######


# nMDS for wall on USVs vs breeding treatments
wallpooled.nMDS<-metaMDS(wall.pooled, distance = "bray")


# Create a variable to assign symbols and colors in the graph
symbols <- groups$Genetic.Status_OO.1_OI.2_II.3
symbols
symbols[symbols==1]<-16 # filled circle
symbols[symbols==2]<-17 # filled triangle
symbols[symbols==3]<-17 # filled triangle
symbols


color <- groups$Genetic.Status_OO.1_OI.2_II.3
color
color[color==1]<-"deepskyblue3" # t?rkis, blue
color[color==2]<-"firebrick1" # middle red
color[color==3]<-"firebrick4" # dark red
color

# for the legend
inbreed.status <- c("Outbred","Inbred 1.gen", "Inbred 2.gen")
legendcolor <- c("deepskyblue3", "firebrick1", "firebrick4")

fig<-plot(wallpooled.nMDS, type = "n")
text(fig, "species", font=c(2), cex=c(1), col = "black")
legend("bottomright", legend = inbreed.status, cex = c(1.3), col = legendcolor, pch = c(16,17,17))
ordispider(fig, groups=groups$Genetic.Status_OO.1_OI.2_II.3, col = legendcolor, lty="longdash") # lines connect to the centroid
#ordihull(fig, groups=groups$Genetic.Status_OO.1_OI.2_II.3, col = legendcolor, lty="dotted") # circle around each group
title("Introduction")


############################################################################

##### ----- comparing Inbred groups (Outbred vs inbred)----- ######

# anosim for wall for Inbred groups
wallpooled.dist <- vegdist(wall.pooled)
wallpooled.ano <- anosim(wallpooled.dist, groups$PooledGenetic_OO.1_OI.II.2)
summary(wallpooled.ano)
plot(wallpooled.ano)

#to test which of the syllable types contributes the most to the observed differences --> which species are responsible for the difference between groups
sim<-simper(wall.pooled, group=groups$PooledGenetic_OO.1_OI.II.2)
sim
summary(sim)


####################  GRAPHS: nMDS-graphs  ##########################

##### ----- comparing 2 Inbred groups (OB vs IB)----- ######


# nMDS for wall on USVs vs inbreeding groups
wallpooled.nMDS<-metaMDS(wall.pooled, distance = "bray")

# Create a variable to assign symbols and colors in the graph
symbols <- groups$PooledGenetic_OO.1_OI.II.2
symbols
symbols[symbols==1]<-16 # filled circle
symbols[symbols==2]<-17 # filled triangle
symbols


color <- groups$PooledGenetic_OO.1_OI.II.2
color
color[color==1]<-"deepskyblue3" # t?rkis, blue
color[color==2]<-"firebrick1" # middle red
color

#  for the legend
inbreed.group <- c("Outbred","Inbred")
legendcolor <- c("deepskyblue3", "firebrick1")

fig<-plot(wallpooled.nMDS, type = "n")
points(fig, "sites", cex=c(1.5),pch = symbols,col = color)
text(fig, "species", font=c(2), cex=c(1), col = "black")
legend("bottomright", legend = inbreed.group, cex = c(1.3), col = legendcolor, pch = c(16,17))
#ordispider(fig, groups=groups$PooledGenetic_OO.1_OI.II.2, col = legendcolor, lty="longdash") # lines connect to the centroid
ordihull(fig, groups=groups$PooledGenetic_OO.1_OI.II.2, col = legendcolor, lty="dotted") # circle around each group
title("Introduction")




##########-----------------------------------#############################--------------------------------################################################


#################################################
### INTERACTION Phase - Using 10 syllable types #################
##################################################



############## STATS - Anosim & Adonis #################

##### ----- comparing 3 breeding treatments (OO vs OI vs II)----- ######

# anosim for contact for breeding treatments
contact.dist <- vegdist(contact)
contact.ano <- anosim(contact.dist, groups$Genetic.Status_OO.1_OI.2_II.3)
summary(contact.ano)
plot(contact.ano)


#to test which of the syllable types contributes the most to the observed differences --> which species are responsible for the difference between groups
sim<-simper(contact, group=groups$Genetic.Status_OO.1_OI.2_II.3)
sim
summary(sim)


####################  GRAPHS: nMDS-graphs  ##########################

##### -----for comparing 3 breeding treatments (OO-OI-II)----- ######


#nMDS for contact on USVs vs inbreeding status
contact.nMDS<-metaMDS(contact, distance = "bray")


# Create a variable to assign symbols and colors in the graph
symbols <- groups$Genetic.Status_OO.1_OI.2_II.3
symbols
symbols[symbols==1]<-16 # filled circle
symbols[symbols==2]<-17 # filled triangle
symbols[symbols==3]<-17 # filled triangle
symbols


color <- groups$Genetic.Status_OO.1_OI.2_II.3
color
color[color==1]<-"deepskyblue3" # t?rkis, blue
color[color==2]<-"firebrick1" # middle red
color[color==3]<-"firebrick4" # dark red
color

# for the legend
inbreed.status <- c("Outbred","Inbred 1.gen", "Inbred 2.gen")
legendcolor <- c("deepskyblue3", "firebrick1", "firebrick4")

fig<-plot(contact.nMDS, type = "n")
points(fig, "sites", cex=c(1.5),pch = symbols,col = color)
text(fig, "species", font=c(2), cex=c(1), col = "black")
legend("bottomright", legend = inbreed.status, cex = c(1.3), col = legendcolor, pch = c(16,17))
#ordispider(fig, groups=groups$Genetic.Status_OO.1_OI.2_II.3, col = legendcolor, lty="longdash") # lines connect to the centroid
ordihull(fig, groups=groups$Genetic.Status_OO.1_OI.2_II.3, col = legendcolor, lty="dotted") # circle around each group
title("Interaction")


############################################################################

##### -----for comparing Outbred vs Inbred groups (Outbred vs inbred)----- ######

# anosim for contact for OB vs IB
contact.dist <- vegdist(contact)
contact.ano <- anosim(contact.dist, groups$PooledGenetic_OO.1_OI.II.2)
summary(contact.ano)
plot(contact.ano)


#to test which of the syllable types contributes the most to the observed differences --> which species are responsible for the difference between groups
sim<-simper(contact, group=groups$PooledGenetic_OO.1_OI.II.2)
sim
summary(sim)


####################  GRAPHS: nMDS-graphs  ##########################

##### -----for comparing Outbred vs Inbred groups (OB vs IB)----- ######


# nMDS for contact on USVs for OB vs IB
contact.nMDS<-metaMDS(contact, distance = "bray")


# Create a variable to assign symbols and colors in the graph
symbols <- groups$PooledGenetic_OO.1_OI.II.2
symbols
symbols[symbols==1]<-16 # filled circle
symbols[symbols==2]<-17 # filled triangle
symbols


color <- groups$PooledGenetic_OO.1_OI.II.2
color
color[color==1]<-"deepskyblue3" # t?rkis, blue
color[color==2]<-"firebrick1" # middle red
color

#for the legend
inbreed.group <- c("Outbred","Inbred")
legendcolor <- c("deepskyblue3", "firebrick1")

fig<-plot(contact.nMDS, type = "n")
points(fig, "sites", cex=c(1.5),pch = symbols,col = color)
text(fig, "species", font=c(2), cex=c(1), col = "black")
legend("bottomright", legend = inbreed.group, cex = c(1.3), col = legendcolor, pch = c(16,17))
#ordispider(fig, groups=groups$PooledGenetic_OO.1_OI.II.2, col = legendcolor, lty="longdash") # lines connect to the centroid
ordihull(fig, groups=groups$PooledGenetic_OO.1_OI.II.2, col = legendcolor, lty="dotted") # circle around each group
title("Interaction")







##########-----------------------------------#############################--------------------------------################################################
##########-----------------------------------#############################--------------------------------################################################
##########-----------------------------------#############################--------------------------------################################################
##########-----------------------------------#############################--------------------------------################################################
##########-----------------------------------#############################--------------------------------################################################
##########-----------------------------------#############################--------------------------------################################################







########################################################################################################################################
#---------------------------------------------------------------------------------------------------------------------------------------
#### R-code to calculate Models to compare spectro-temporal features for single USV datapoints (Long table) between breeding treatments (O vs I vs II) or inbred groups (Inbred vs Outbred) for the Interaction phase
#---------------------------------------------------------------------------------------------------------------------------------------
########################################################################################################################################


library(MuMIn)
library(MASS)

library(multcomp)
library(lme4)
library(car)
library(GeneralizedHyperbolic)
library(actuar)

library('fitdistrplus')

setwd("C:/Documents/stats")
Longdata=read.csv("LongTable.csv",na.strings="", dec = ".", sep=";")
View(Longdata)

Longdata$genetic1=as.factor(Longdata$Genetic)
levels(Longdata$genetic1)=c("OO","OI","II")

Longdata$genPooled=as.factor(Longdata$PooledGenetic)
levels(Longdata$genPooled)=c("Outbred","Inbred")



library(fiwidat)


################################################################################
#  Length - USV length [ms]


#plot the data
plot(Longdata$Length_ms~Longdata$genetic1,log="y",ylab="Length (ms)",xlab="Genetic status",cex.lab=1.5) # y axis is log-based (not transformed the values, just show raw values on log axis)
plot(Longdata$Length_ms~Longdata$genetic1,ylab="Length (ms)",xlab="Genetic status",cex.lab=1.5) # normal y axis

tapply(Longdata$Length_ms,Longdata$genetic1,mean) # calculate the mean for each group
#      OO       OI       II 
# 33.86564 31.63200 39.26674 

tapply(Longdata$Length_ms,Longdata$genetic1,median)  # calculate the median for each group
#    OO     OI     II 
# 31.060 27.530 33.565 


sem <- function(x) sqrt(var(x)/length(x))

tapply(Longdata$Length_ms,Longdata$genetic1,sem) # calculate the standard error of the mean for each group
#      OO        OI        II 
# 0.2580223 0.2803565 0.3053640 


hist(Longdata$Length_ms)


###calculate model with 3 genetic groups
model=glmer(Length_ms~genetic1+(1|Replicate),data=Longdata,family="inverse.gaussian") # calculate generalized mixed model with genetic1 as fixed effect and Replicate as random effect
summary(model)

# Post Hoc Test for multiple comparison:
summary(glht(model, mcp(genetic1="Tukey"))) 

resi=residuals(model)
hist(resi) 
plot(fitdist(resi+3, "invgauss", start = list(mean = 3, shape = 1))) 

# calculate Null-Model
model0=glmer(Length_ms~1+(1|Replicate), data=Longdata,family="inverse.gaussian")
summary(model0)


### calculate model with 2 Groups (IB vs OB)
model2=glmer(Length_ms~genPooled+(1|Replicate),data=Longdata,family="inverse.gaussian")
summary(model2)

# Post Hoc Test for multiple comparison for model 2:
summary(glht(model2, mcp(genPooled="Tukey"))) 


tapply(Longdata$Length_ms,Longdata$genPooled,mean)
tapply(Longdata$Length_ms,Longdata$genPooled,sem)
# Outbred   Inbred 
# 33.86564 35.96169 
# 0.2580223 0.2141663 


AICc(model0) 
AICc(model)  
AICc(model2)

# calculate glm model without replicate as random effect
modelg=glm (Length_ms~genetic1,data=Longdata,family='inverse.gaussian')
summary(modelg)

AICc(modelg)
AICc(modelg)-AICc(model)
anova(model,modelg)

par(mar=c(5,5,4,4))
plot(Length_ms~genetic1,data=Longdata,log="y",xlab="Genetic status",cex.lab=1.5,ylab="USV length [ms]", ylim=c(1,700))

par(mar=c(5,5,4,4))
plot(Length_ms~genetic1,data=Longdata,log="y",xlab="Genetic status",cex.lab=1.5,ylab="USV length [ms]", ylim=c(1,700))




################################################################################
# fmean - Mean Frequency [kHz]

hist(Longdata$fMean)
qqnorm(Longdata$fMean)

#plot the data
plot(Longdata$fMean~Longdata$genetic1,log="y",ylab="fMean",xlab="Genetic status",cex.lab=1.5) # y axis is log-based (not transformed the values, just show raw velues on log axis)
plot(Longdata$fMean~Longdata$genetic1,ylab="fMean",xlab="Genetic status",cex.lab=1.5) # normal y axis

tapply(Longdata$fMean,Longdata$genetic1,mean,na.rm=T)
#      OO       OI       II 
# 65.41552 65.39148 63.53678 

tapply(Longdata$fMean,Longdata$genetic1,sd,na.rm=T)
#     OO       OI       II 
# 8.232778 9.685964 8.401025 

tapply(Longdata$fMean,Longdata$genetic1,median,na.rm=T)
#   OO    OI    II 
# 66.15 64.58 63.20 


### calculate model with 3 genetic groups
model=glmer ((Longdata$fMean/100)~genetic1+(1|Replicate),data=Longdata,family='inverse.gaussian',
             control=glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

summary(model)

resi=residuals(model)
hist(resi) 

plot(fitdist(resi+60,"invgauss",start=list(mean=2,shape=1))) 

# Post Hoc Test for multiple comparison:
summary(glht(model, mcp(genetic1="Tukey")))

# calculate Null model 
model0=glmer ((Longdata$fMean/100)~1+(1|Replicate),data=Longdata,family='inverse.gaussian',
              control=glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))
summary(model0)



### calculate model with 2 Groups (IB vs OB)
model2=glmer ((Longdata$fMean/100)~genPooled+(1|Replicate),data=Longdata,family='inverse.gaussian',
              control=glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))
summary(model2)

# Post Hoc Test for multiple comparison for model2:
summary(glht(model2, mcp(genPooled="Tukey"))) 


AICc(model)  
AICc(model0) 
AICc(model2) 

# calculate glm model without replicate as random effect
modelg=glm((Longdata$fMean/100)~genetic1,data=Longdata,family='inverse.gaussian')
summary (modelg)
AICc(modelg)
anova(model,modelg)
AICc(modelg)-AICc(model)


plot(Longdata$fMean~Longdata$genetic1,xlab="Genetic status",ylab="Mean Frequency (kHZ)")

tapply(Longdata$fMean,Longdata$genPooled,mean,na.rm=T)
# Outbred   Inbred 
# 65.41552 64.33165 

tapply(Longdata$fMean,Longdata$genPooled,sd,na.rm=T)
# Outbred   Inbred 
# 8.232778 9.020703 

leveneTest(Longdata$fMean~Longdata$genetic1)
leveneTest(Longdata$fMean~Longdata$genPooled)

fMeandev=abs(Longdata$fMean-median(Longdata$fMean,na.rm=T))
devres=fMeandev+1

hist(devres)

modelv=glmer(devres~genetic1+(1|Replicate), family="Gamma",data=Longdata)
summary(modelv)
modelv0=glmer(devres~1+(1|Replicate), family="Gamma",data=Longdata)
modelv2=glmer(devres~genPooled+(1|Replicate), family="Gamma",data=Longdata)

AICc(modelv0) 
AICc(modelv) 
AICc(modelv2)

par(mar=c(5,5,4,4))
plot(fMeandev~genetic1,data=Longdata,xlab="Genetic status",cex.lab=1.5,ylab="Variation in mean frequency",xaxt= "n", ylim=c(-5,70))




################################################################################
#  ISI - Inter-syllable interval (s)


#plot the data
plot(Longdata$ISI~Longdata$genetic1,log="y",ylab="ISI (s)",xlab="Genetic status",cex.lab=1.5) # y axis is log-based (not transformed the values, just show raw velues on log axis)
plot(Longdata$ISI~Longdata$genetic1,ylab="ISI (s)",xlab="Genetic status",cex.lab=1.5) # normal y axis

tapply(Longdata$ISI,Longdata$genetic1,mean,na.rm=TRUE) # calculate the mean for each group, after removing NA
#      OO       OI       II 
# 0.1116960 0.1177070 0.1130647  

tapply(Longdata$ISI,Longdata$genetic1,median,na.rm=TRUE)  # calculate the median for each group, after removing NA
#       OO       OI       II 
#   0.079620 0.086995 0.080265   

hist(Longdata$ISI)


### calculate model with 3 genetic groups
model=glmer(ISI~genetic1+(1|Replicate),data=Longdata,family="inverse.gaussian") # calculate generalized mixed model with genetic1 as fixed effect and Replicate as random effect
summary(model)

# Post Hoc Test for multiple comparison:
summary(glht(model, mcp(genetic1="Tukey")))

resi=residuals(model)
hist(resi) 
plot(fitdist(resi+3, "invgauss", start = list(mean = 3, shape = 1))) 

# calculate Null-Model
model0=glmer(ISI~1+(1|Replicate), data=Longdata,family="inverse.gaussian")
summary(model0)


### calculate model with 2 Groups (IB vs OB)
model2=glmer(ISI~genPooled+(1|Replicate),data=Longdata,family="inverse.gaussian")
summary(model2)


# Post Hoc Test for multiple comparison for model 2:
summary(glht(model2, mcp(genPooled="Tukey")))


AICc(model0)
AICc(model)  
AICc(model2)

### calculate glm model without replicate as random effect
modelg=glm (ISI~genetic1,data=Longdata,family='inverse.gaussian')
summary(modelg)

AICc(modelg)
AICc(modelg)-AICc(model)
anova(model,modelg)


anova(model,model0)

par(mar=c(5,5,4,4))
plot(ISI~genetic1,data=Longdata,log="y",xlab="Genetic status",cex.lab=1.5,ylab="Inter syllable interval (ISI) [s]")






################################################################################
#  IBI - Inter-bout interval (s)


#plot the data
plot(Longdata$IBI~Longdata$genetic1,log="y",ylab="IBI (s)",xlab="Genetic status",cex.lab=1.5) # y axis is log-based (not transformed the values, just show raw velues on log axis)
plot(Longdata$IBI~Longdata$genetic1,ylab="IBI (s)",xlab="Genetic status",cex.lab=1.5) # normal y axis

tapply(Longdata$IBI,Longdata$genetic1,mean,na.rm=TRUE) # calculate the mean for each group, after removing NA
#      OO       OI       II 
# 4.671909 5.127112 3.969187 

tapply(Longdata$IBI,Longdata$genetic1,median,na.rm=TRUE)  # calculate the median for each group, after removing NA
#       OO       OI       II 
#   0.885500 0.831495 0.724710  


hist(Longdata$IBI) 


### calculate model with 3 genetic groups
model=glmer(IBI~genetic1+(1|Replicate),data=Longdata,family="inverse.gaussian") # calculate generalized mixed model with genetic1 as fixed effect and Replicate as random effect
summary(model)

# Post Hoc Test for multiple comparison:
summary(glht(model, mcp(genetic1="Tukey")))

resi=residuals(model)
hist(resi) # check residuals - look quasi/nearly normal distributed
plot(fitdist(resi+3, "invgauss", start = list(mean = 3, shape = 1))) # all plots similar to priming data

# calculate Null-Model
model0=glmer(IBI~1+(1|Replicate), data=Longdata,family="inverse.gaussian")
summary(model0)



### calculate model with 2 Groups (IB vs OB)
model2=glmer(IBI~genPooled+(1|Replicate),data=Longdata,family="inverse.gaussian")
summary(model2)

# Post Hoc Test for multiple comparison for model 2:
summary(glht(model2, mcp(genPooled="Tukey"))) 



AICc(model0) 
AICc(model) 
AICc(model2) 

# calculate glm model without replicate as random effect
modelg=glm (IBI~genetic1,data=Longdata,family='inverse.gaussian')
summary(modelg)

AICc(modelg)
AICc(modelg)-AICc(model)
anova(model,modelg)


anova(model,model0)

par(mar=c(5,5,4,4))
plot(IBI~genetic1,data=Longdata,log="y",xlab="Genetic status",cex.lab=1.5,ylab="Inter bout interval (IBI) [s]", ylim=c(1,700))