LARVAE
#use the raw data:
larv <- read.csv(file="C:/Users/rjolma/OneDrive - NIOZ/Databases and forms/Larvae07092023.csv",header=TRUE, sep=",")
#the time from start to the detection of the first infectious stage in each temperature (same steps also for other responses):
larv$date <- as.Date(larv$date, format = "%d.%m.%Y")
library(lubridate)
seasons = function(x){
  if(x %in% 2:5) return("Spring")
  if(x %in% 6:7) return("Summer")
  if(x %in% 8:10) return("Fall")
  if(x %in% c(11,12,1)) return("Winter")
  
}
larv$Season = sapply(month(larv$date), seasons)

larv$red10 <- as.Date(larv$red10, format = "%d.%m.%Y")
larv$red14 <- as.Date(larv$red14, format = "%d.%m.%Y")
larv$red18 <- as.Date(larv$red18, format = "%d.%m.%Y")
larv$red22 <- as.Date(larv$red22, format = "%d.%m.%Y") 
larv$red26 <- as.Date(larv$red26, format = "%d.%m.%Y")

larv$fh10 <- as.Date(larv$fh10, format = "%d.%m.%Y")
larv$fh14 <- as.Date(larv$fh14, format = "%d.%m.%Y")
larv$fh18 <- as.Date(larv$fh18, format = "%d.%m.%Y")
larv$fh22 <- as.Date(larv$fh22, format = "%d.%m.%Y")
larv$fh26 <- as.Date(larv$fh26, format = "%d.%m.%Y")

larv$fc10 <- as.Date(larv$fc10, format = "%d.%m.%Y")
larv$fc14 <- as.Date(larv$fc14, format = "%d.%m.%Y")
larv$fc18 <- as.Date(larv$fc18, format = "%d.%m.%Y")
larv$fc22 <- as.Date(larv$fc22, format = "%d.%m.%Y")
larv$fc26 <- as.Date(larv$fc26, format = "%d.%m.%Y")

larv$ad10 <- as.Date(larv$ad10, format = "%d.%m.%Y")
larv$ad14 <- as.Date(larv$ad14, format = "%d.%m.%Y")
larv$ad18 <- as.Date(larv$ad18, format = "%d.%m.%Y")
larv$ad22 <- as.Date(larv$ad22, format = "%d.%m.%Y")
larv$ad26 <- as.Date(larv$ad26, format = "%d.%m.%Y")

larv$fh10d<-as.numeric(difftime(larv$fh10,larv$date))
larv$fh14d<-as.numeric(difftime(larv$fh14,larv$date))
larv$fh18d<-as.numeric(difftime(larv$fh18,larv$date))
larv$fh22d<-as.numeric(difftime(larv$fh22,larv$date))
larv$fh26d<-as.numeric(difftime(larv$fh26,larv$date))

larv$fc10d<-as.numeric(difftime(larv$fc10,larv$date))
larv$fc14d<-as.numeric(difftime(larv$fc14,larv$date))
larv$fc18d<-as.numeric(difftime(larv$fc18,larv$date))
larv$fc22d<-as.numeric(difftime(larv$fc22,larv$date))
larv$fc26d<-as.numeric(difftime(larv$fc26,larv$date))

larv$ad10d<-as.numeric(difftime(larv$ad10,larv$date))
larv$ad14d<-as.numeric(difftime(larv$ad14,larv$date))
larv$ad18d<-as.numeric(difftime(larv$ad18,larv$date))
larv$ad22d<-as.numeric(difftime(larv$ad22,larv$date))
larv$ad26d<-as.numeric(difftime(larv$ad26,larv$date))

#one came out as seconds and was then specified
larv$fh18d<-as.numeric(difftime(larv$fh18,larv$date, units = "days"))

larv$X50d10 <- as.Date(larv$X50d10, format = "%d.%m.%Y")
larv$X50d14 <- as.Date(larv$X50d14, format = "%d.%m.%Y")
larv$X50d18 <- as.Date(larv$X50d18, format = "%d.%m.%Y")
larv$X50d22 <- as.Date(larv$X50d22, format = "%d.%m.%Y")
larv$X50d26 <- as.Date(larv$X50d26, format = "%d.%m.%Y")

larv$X50d10d<-as.numeric(difftime(larv$X50d10,larv$date))
larv$X50d14d<-as.numeric(difftime(larv$X50d14,larv$date))
larv$X50d18d<-as.numeric(difftime(larv$X50d18,larv$date))
larv$X50d22d<-as.numeric(difftime(larv$X50d22,larv$date))
larv$X50d26d<-as.numeric(difftime(larv$X50d26,larv$date))

#time from first copepodite to 50% dead:
larv$life10<-as.numeric(difftime(larv$X50d10,larv$fc10, units = "days"))
larv$life14<-as.numeric(difftime(larv$X50d14,larv$fc14, units = "days"))
larv$life18<-as.numeric(difftime(larv$X50d18,larv$fc18, units = "days"))
larv$life22<-as.numeric(difftime(larv$X50d22,larv$fc22, units = "days"))
larv$life26<-as.numeric(difftime(larv$X50d26,larv$fc26, units = "days"))


#Go from each egg sac (“number”) on a separate row to each observation on a separate row. Thus 5 rows per egg sac in the new format. We also want to keep information about the parasite species and season onboard:
library(reshape2)
dat.c<- melt(larv,id.vars=c("number", "ori", "Season"), measure.vars=c('fc10d','fc14d','fc18d','fc22d','fc26d'), variable.name = "temperature", 
             value.name = "time")
dat.h <- melt(larv,id.vars=c("number", "ori", "Season"), measure.vars=c('fh10d','fh14d','fh18d','fh22d','fh26d'), variable.name = "temperature", value.name = "time")


dat.d <- melt(larv,id.vars=c("number", "ori", "Season"), measure.vars=c('ad10d','ad14d','ad18d','ad22d','ad26d'), variable.name = "temperature", 
              value.name = "time")
dat.50d <- melt(larv,id.vars=c("number", "ori", "Season"), measure.vars=c('X50d10d','X50d14d','X50d18d','X50d22d','X50d26d'), variable.name = "temperature", 
                value.name = "time")

dat.a <- melt(larv,id.vars=c("number", "ori", "Season"), measure.vars=c('life10','life14','life18','life22','life26'), variable.name = "temperature", value.name = "time")


#make a success column in which all observations that have a number (a infectious stage appeared at some point) are 1 and all with NA are 0 (no infectious stage seen)
dat.c$time<-as.numeric(dat.c$value)
dat.c$success<-dat.c$time
library(dplyr)
dat.c <- dat.c %>% mutate(success = ifelse(is.na(success), 0, success))
dat.c$success[dat.c$success>0]<-1
#species separated:

#success models :
#first looking at the effect of season separately:
fisher.test(ori$success, ori$Season)
fisher.test(int$success, int$Season)

library(lme4)
sucmodel3<-glmer(success ~ ori * temperature+ (1|number), family=binomial,data=dat.c,nAGQ=2)
summary(sucmodel3)

sucmodel4<-glmer(success ~ ori + temperature+ (1|number), family=binomial,data=dat.c,nAGQ=2)
summary(sucmodel4)

anova(sucmodel3, sucmodel4)

ori<-subset(dat.c, ori==1)
int<-subset(dat.c, ori==0)

sucmodeli<-glmer(success ~  temperature+ (1|number),family=binomial,data=int,nAGQ=2)
summary(sucmodeli)

sucmodelo<-glmer(success ~  temperature+ (1|number), family=binomial,data=ori,nAGQ=2)
summary(sucmodelo)

#Poisson models:
#first hatch


dat.h$ori<-as.factor(dat.h$ori)
summary(h1 <- glmer(time ~ temperature + ori + (1 | number), family="poisson", data=dat.h))
summary(h2 <- glmer(time ~ temperature * ori + (1 | number), family="poisson", data=dat.h))
anova(h1, h2)


#First copepodite
dat.c$ori<-as.factor(dat.c$ori)
summary(c1 <- glmer(time ~ temperature + ori + (1 | number), family="poisson", data=dat.c))
summary(c2 <- glmer(time ~ temperature * ori + (1 | number), family="poisson", data=dat.c))
anova(c1, c2)

#50% dead:
dat.50d$ori<-as.factor(dat.50d$ori)
summary(d1 <- glmer(time ~ temperature + ori + (1 | number), family="poisson", data=dat.50d))
summary(d2 <- glmer(time ~ temperature * ori + (1 | number), family="poisson", data=dat.50d))
anova (d1, d2)

# copepodite lifespan
dat.a$ori<-as.factor(dat.a$ori)
summary(cl1 <- glmer(time ~ temperature + ori + (1 | number), family="poisson", data=dat.a))
summary(cl2 <- glmer(time ~ temperature * ori + (1 | number), family="poisson", data=dat.a))
anova(cl1, cl2)

dat.c$temperature<- as.character(dat.c$temperature)
dat.c["temperature"][dat.c["temperature"] == "f10d"] <- "10"
dat.c["temperature"][dat.c["temperature"] == "f14d"] <- "14"
dat.c["temperature"][dat.c["temperature"] == "f18d"] <- "18"
dat.c["temperature"][dat.c["temperature"] == "f22d"] <- "22"
dat.c["temperature"][dat.c["temperature"] == "f26d"] <- "26"

#separate poissons by species

dat.hi<- subset(dat.h, ori==0)
dat.ho<- subset(dat.h, ori==1)

dat.ci<- subset(dat.c, ori==0)
dat.co<- subset(dat.c, ori==1)

dat.50di<- subset(dat.50d, ori==0)
dat.50do<- subset(dat.50d, ori==1)

dat.ai<- subset(dat.a, ori==0)
dat.ao<- subset(dat.a, ori==1)

summary(hi1 <- glmer(time ~ temperature + (1 | number), family="poisson", data=dat.hi))
summary(ho1 <- glmer(time ~ temperature + (1 | number), family="poisson", data=dat.ho))

summary(ci1 <- glmer(time ~ temperature + (1 | number), family="poisson", data=dat.ci))
summary(co1 <- glmer(time ~ temperature + (1 | number), family="poisson", data=dat.co))

summary(d50i1 <- glmer(time ~ temperature + (1 | number), family="poisson", data=dat.50di))
summary(d50o1 <- glmer(time ~ temperature + (1 | number), family="poisson", data=dat.50do))

summary(ai1 <- glmer(time ~ temperature + (1 | number), family="poisson", data=dat.ai))
summary(ao1 <- glmer(time ~ temperature + (1 | number), family="poisson", data=dat.ao))

#Figure 2 (response times)
dat.h <- melt(larv,id.vars=c("number", "ori"), measure.vars=c('fh10d','fh14d','fh18d','fh22d','fh26d'), variable.name = "temperature", value.name = "time")
dat.d <- melt(larv,id.vars=c("number", "ori"), measure.vars=c('ad10d','ad14d','ad18d','ad22d','ad26d'), variable.name = "temperature", value.name = "time")
dat.50d <- melt(larv,id.vars=c("number", "ori"), measure.vars=c('X50d10d','X50d14d','X50d18d','X50d22d','X50d26d'), variable.name = "temperature",  value.name = "time")

names(dat.h)[names(dat.h) == "ori"] <- "Species"
dat.h$Species<-as.character(dat.h$Species)
dat.h["Species"][dat.h["Species"] == 0] <- "M. intestinalis"
dat.h["Species"][dat.h["Species"] == 1] <- "M. orientalis"

names(dat.50d)[names(dat.50d) == "ori"] <- "Species"
dat.50d$Species<-as.character(dat.50d$Species)
dat.50d["Species"][dat.50d["Species"] == 0] <- "M. intestinalis"
dat.50d["Species"][dat.50d["Species"] == 1] <- "M. orientalis"

names(dat.d)[names(dat.d) == "ori"] <- "Species"
dat.d$Species<-as.character(dat.d$Species)
dat.d["Species"][dat.d["Species"] == 0] <- "M. intestinalis"
dat.d["Species"][dat.d["Species"] == 1] <- "M. orientalis"

dat.50d$temperature<- as.character(dat.50d$temperature)
dat.50d["temperature"][dat.50d["temperature"] == "X50d10d"] <- "10"
dat.50d["temperature"][dat.50d["temperature"] == "X50d14d"] <- "14"
dat.50d["temperature"][dat.50d["temperature"] == "X50d18d"] <- "18"
dat.50d["temperature"][dat.50d["temperature"] == "X50d22d"] <- "22"
dat.50d["temperature"][dat.50d["temperature"] == "X50d26d"] <- "26"
dat.50d$temperature<- as.factor(dat.50d$temperature)
dat.50d$Species<-as.factor(dat.50d$Species)

dat.h$temperature<- as.character(dat.h$temperature)
dat.h["temperature"][dat.h["temperature"] == "fh10d"] <- "10"
dat.h["temperature"][dat.h["temperature"] == "fh14d"] <- "14"
dat.h["temperature"][dat.h["temperature"] == "fh18d"] <- "18"
dat.h["temperature"][dat.h["temperature"] == "fh22d"] <- "22"
dat.h["temperature"][dat.h["temperature"] == "fh26d"] <- "26"
dat.h$temperature<- as.factor(dat.h$temperature)
dat.h$Species<-as.factor(dat.h$Species)

dat.d$temperature<- as.character(dat.d$temperature)
dat.d["temperature"][dat.d["temperature"] == "ad10d"] <- "10"
dat.d["temperature"][dat.d["temperature"] == "ad14d"] <- "14"
dat.d["temperature"][dat.d["temperature"] == "ad18d"] <- "18"
dat.d["temperature"][dat.d["temperature"] == "ad22d"] <- "22"
dat.d["temperature"][dat.d["temperature"] == "ad26d"] <- "26"
dat.d$temperature<- as.factor(dat.d$temperature)
dat.d$Species<-as.factor(dat.d$Species)


ph<-ggplot(dat.h, aes(x = temperature, y = time, fill=Species)) +
  labs(x = "Temperature (°C)", y = "Days") +
  theme(text = element_text(size = 15))   +
  scale_color_manual(values=c("blue3","orange2")) +
  geom_boxplot()+ 
  ggtitle("First hatching") +
  ylim(0, 50)
ph
pc<-ggplot(dat.c, aes(x = tempf, y = value, fill=Species)) +
  labs(x = "Temperature (°C)", y = "Days") +
  theme(text = element_text(size = 15))   +
  scale_color_manual(values=c("blue3","orange2")) +
  geom_boxplot()+ 
  ggtitle("Infectious stage") +
  ylim(0, 50) 
pc
pd<-ggplot(dat.d, aes(x = temperature, y = time, fill=Species)) +
  labs(x = "Temperature (°C)", y = "Days") +
  theme(text = element_text(size = 15))   +
  scale_color_manual(values=c("blue3","orange2")) +
  geom_boxplot()+ 
  ggtitle("All dead") +
  ylim(0, 50)
pd
p50d<-ggplot(dat.50d, aes(x = temperature, y = time, fill=Species)) +
  labs(x = "Temperature (°C)", y = "Days") +
  theme(text = element_text(size = 15))   +
  scale_color_manual(values=c("blue3","orange2")) +
  geom_boxplot()+ 
  ggtitle("50% dead") +
  ylim(0, 50)
p50d
ggarrange(ph, pc, p50d, pd,
          ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom")

#lifespan:
#time from first copepodite to 50% dead:
larv$life10<-as.numeric(difftime(larv$X50d10,larv$fc10, units = "days"))
larv$life14<-as.numeric(difftime(larv$X50d14,larv$fc14, units = "days"))
larv$life18<-as.numeric(difftime(larv$X50d18,larv$fc18, units = "days"))
larv$life22<-as.numeric(difftime(larv$X50d22,larv$fc22, units = "days"))
larv$life26<-as.numeric(difftime(larv$X50d26,larv$fc26, units = "days"))

library(reshape2)
dat.a <- melt(larv,id.vars=c("number", "ori", "Season"), measure.vars=c('life10','life14','life18','life22','life26'), variable.name = "temperature", value.name = "time")

dat.a <- melt(larv,id.vars=c("number", "ori", "Season"), measure.vars=c('life10','life14','life18','life22','life26'), variable.name = "temperature", value.name = "time")
dat.a$temperature <- sub("life", "", dat.a$temperature)
str(dat.a)
dat.a$temperature<- as.factor(dat.a$temperature)
dat.a$Season<- as.factor(dat.a$Season)
dat.a$ori<- as.factor(dat.a$ori)
pa<-ggplot(dat.a, aes(x = temperature, y = time, fill=Species)) +
  labs(x = "Temperature (°C)", y = "Days") +
  theme(text = element_text(size = 15))   +
  scale_color_manual(values=c("blue3","orange2")) +
  geom_boxplot()+ 
  ggtitle("Copepodite lifespan") +
  ylim(0, 50)
pa





















#STAGES INSIDE HOSTS:  

par <- read.csv(file="C:/Users/rjolma/OneDrive - NIOZ/Documenten/Project planning/R digging/adults.csv",header=TRUE, sep=",")
par$date <- as.Date(par$date, format = "%d/%m/%Y")
par$end_date <- as.Date(par$end_date, format = "%d/%m/%Y")
par$days<-difftime(par$end_date,par$date)
par$days<-as.numeric(par$days)
par$tempi <- par$temp
par$tempi <- as.factor(par$tempi)
par$plength <- as.numeric(par$plength)
parin<- subset(par, in24h==1)
parinalive<- subset(parin, died==0)
pars<-subset(parinalive, inend==1)
pars$g<-as.factor(pars$g)

aso<-subset(pars, ori==1)
asi<-subset(pars, int==1)

aso10<- subset(aso, temp==10)
aso14<- subset(aso, temp==14)
aso18<- subset(aso, temp==18)
aso22<- subset(aso, temp==22)
aso26<- subset(aso, temp==26)

asi10<- subset(asi, temp==10)
asi14<- subset(asi, temp==14)
asi18<- subset(asi, temp==18)
asi22<- subset(asi, temp==22)
asi26<- subset(asi, temp==26)
#host entry success

#with linear temperature:
sucf242<-glmer(in24h ~ ori * temp + (1|number), family=binomial, data=par)
sucf243<-glmer(in24h ~ ori + temp + (1|number), family=binomial, data=par)
anova(sucf242, sucf243)
summary(sucf242)
summary(sucf243)

#with categorical temperature:
sucfl<-glmer(in24h ~ ori * tempi + (1|number), family=binomial, data=par)
sucfl2<-glmer(in24h ~ ori + tempi + (1|number), family=binomial, data=par)
anova(sucfl, sucfl2)
summary(sucfl)
summary(sucfl2)

#by species:
orip<-subset(par, ori==1)
intp<-subset(par, ori==0)

suci<-glmer(in24h ~  tempi + (1|number), family=binomial, data=intp)
suco<-glmer(in24h ~  tempi + (1|number), family=binomial, data=orip)
summary(suci)
summary(suco)

#confidence intervals:
cc24 <-confint(sucfl2,parm="beta_",method="Wald")
ctab24 <- cbind(est=fixef(sucfl2),cc24)
rtab24 <- exp(ctab24)
print(rtab24, digits=3)

cc24i <-confint(suci,parm="beta_",method="Wald")
ctab24i <- cbind(est=fixef(suci),cc24i)
rtab24i <- exp(ctab24i)
print(rtab24i, digits=3)

cc24o <-confint(suco,parm="beta_",method="Wald")
ctab24o <- cbind(est=fixef(suco),cc24o)
rtab24o <- exp(ctab24o)
print(rtab24o, digits=3)

#infection success:
#removed parasites that were clearly of second generation:

parinalive1 <- parinalive[(parinalive$X1stgen == 1 | is.na(parinalive$X1stgen) | parinalive$X1stgen != 0), ]
str(parinalive1)

#with linear temperature:
sucin1<-glmer(inend ~ ori * temp + (1|number), family=binomial, data=parinalive1)
sucin2<-glmer(inend ~ ori + temp + (1|number), family=binomial, data=parinalive1)
anova(sucin1, sucin2)
summary(sucin1)
summary(sucin2)

#with categorical temperature:
sucin3<-glmer(inend ~ ori * tempi + (1|number), family=binomial, data=parinalive1)
sucin4<-glmer(inend ~ ori + tempi + (1|number), family=binomial, data=parinalive1)
anova(sucin3, sucin4)
summary(sucin3)
summary(sucin4)

#by species:
oripin<-subset(parinalive1, ori==1)
intpin<-subset(parinalive1, ori==0)

suciin<-glmer(inend ~  tempi + (1|number), family=binomial, data=intpin)
sucoin<-glmer(inend ~  tempi + (1|number), family=binomial, data=oripin)
summary(suciin)
summary(sucoin)

#confidence intervals:
ccin <-confint(sucin4,parm="beta_",method="Wald")
ctabin<- cbind(est=fixef(sucin4),ccin)
rtabin <- exp(ctabin)
print(rtabin, digits=3)

cci <-confint(suciin,parm="beta_",method="Wald")
ctabi <- cbind(est=fixef(suciin),cci)
rtabi <- exp(ctabi)
print(rtabi, digits=3)

cco <-confint(sucoin,parm="beta_",method="Wald")
ctabo <- cbind(est=fixef(sucoin),cco)
rtabo <- exp(ctabo)
print(rtabo, digits=3)


#Figure 3:

p1<- ggplot(aso10, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle(expression(italic("Mytilicola orientalis") ~ "10°C")) +
  scale_color_manual(values=c("#00BFC4", "purple")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(50, 150) +
  geom_hline(yintercept = 3.55, linetype = "dashed")
p1
p2<- ggplot(aso14, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("14°C") +
  scale_color_manual(values=c("#00BFC4", "purple")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(50, 150) +
  geom_hline(yintercept = 3.55, linetype = "dashed")

p3<- ggplot(aso18, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("18°C") +
  scale_color_manual(values=c("#00BFC4", "purple")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(50, 150) +
  geom_hline(yintercept = 3.55, linetype = "dashed")

p4<- ggplot(aso22, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("22°C") +
  scale_color_manual(values=c("#00BFC4", "purple")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(50, 150) +
  geom_hline(yintercept = 3.55, linetype = "dashed")

p5<- ggplot(aso26, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("26°C") +
  scale_color_manual(values=c("#00BFC4", "purple")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(50, 150) +
  geom_hline(yintercept = 3.55, linetype = "dashed")



p6 <- ggplot(asi10, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle(expression(italic("Mytilicola intestinalis") ~ "10°C")) +
  scale_color_manual(values=c("#F8766D","purple")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8) +
  xlim(50, 150) +
  geom_hline(yintercept = 4.23, linetype = "dashed")

p6

p6
p7<- ggplot(asi14, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("14°C") +
  scale_color_manual(values=c("#F8766D","purple")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(50, 150) +
  geom_hline(yintercept = 4.23, linetype = "dashed")
p7
p8<- ggplot(asi18, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("18°C") +
  scale_color_manual(values=c("#F8766D","purple")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(50, 150) +
  geom_hline(yintercept = 4.23, linetype = "dashed")
p8
p9<- ggplot(asi22, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("22°C") +
  scale_color_manual(values=c("#F8766D","purple")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(50, 150) +
  geom_hline(yintercept = 4.23, linetype = "dashed")

p10<- ggplot(asi26, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("26°C") +
  scale_color_manual(values=c("#F8766D", "purple")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(50, 150) +
  geom_hline(yintercept = 4.23, linetype = "dashed")

ggarrange(p6, p1, p7, p2, p8, p3, p9, p4, p10, p5,
          ncol = 2, nrow = 5, common.legend = TRUE, legend = "bottom")
# exported as a tiff with dimensions of 700x1100

#LENGTH REGRESSION:

#subsets for until reproductive size
#theoretical as does not get there: aso10r = aso10
aso14r<- subset(aso14, days<100)
aso18r<- subset(aso18, days<57)
aso22r<- subset(aso22, days<57)
aso26r<- subset(aso26, days<57)

asi10r<- subset(asi10, days<122)
asi14r<- subset(asi14, days<122)
asi18r<- subset(asi18, days<78)
asi22r<- subset(asi22, days<57)
asi26r<- subset(asi26, days<57)

#subsets for until maximum size 
aso10m<-subset(aso10, days<122)
aso14m<- subset(aso14, days<140)
aso18m<- subset(aso18, days<122)
aso22m<- subset(aso22, days<122)
aso26m<- subset(aso26, days<57)

asi10m<- subset(asi10, days<140)
asi14m<- subset(asi14, days<122) #same as repro
asi18m<- subset(asi18, days<122)
asi22m<- subset(asi22, days<140)
asi26m<- subset(asi26, days<122)

#with also the non-relevant:
#until the reproductive size:
p1<- ggplot(aso10, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 10C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 3.55)

p2<- ggplot(aso14r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 14C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 3.55)

p3<- ggplot(aso18r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 18C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 3.55)

p4<- ggplot(aso22r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 22C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 3.55)

p5<- ggplot(aso26r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 26C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 3.55)

p6<- ggplot(asi10r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 10C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 4.23)

p7<- ggplot(asi14r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 14C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 4.23)

p8<- ggplot(asi18r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 18C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 4.23)

p9<- ggplot(asi22r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 22C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 4.23)

p10<- ggplot(asi26r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 26C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 4.23)
p1
p2
p3
p4
p5
p6
p7
p8
p9
p10

ggarrange(p1, p6, p2, p7, p3, p8, p4, p9, p5, p10,
          ncol = 2, nrow = 5, common.legend = TRUE, legend = "bottom")

#until maximal size:
p1<- ggplot(aso10m, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 10C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 3.55)

p2<- ggplot(aso14m, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 14C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 3.55)

p3<- ggplot(aso18m, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 18C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 3.55)

p4<- ggplot(aso22m, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 22C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 3.55)

p5<- ggplot(aso26m, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 26C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 3.55)

p6<- ggplot(asi10m, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 10C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 4.23)

p7<- ggplot(asi14m, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 14C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 4.23)

p8<- ggplot(asi18m, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 18C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 4.23)

p9<- ggplot(asi22m, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 22C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 4.23)

p10<- ggplot(asi26m, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 26C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_hline(yintercept = 4.23)
p1
p2
p3
p4
p5
p6
p7
p8
p9
p10

ggarrange(p1, p6, p2, p7, p3, p8, p4, p9, p5, p10,
          ncol = 2, nrow = 5, common.legend = TRUE, legend = "bottom")

#linear regression with and without log transformation:

intercept <- 0.27
lmo10 <- lm(I(plength - intercept) ~ 0 + days, aso10)
summary(lmo10)
lmo14r <- lm(I(plength - intercept) ~ 0 + days, aso14r)
summary(lmo14r)
lmo18r <- lm(I(plength - intercept) ~ 0 + days, aso18r)
summary(lmo18r)
lmo22r <- lm(I(plength - intercept) ~ 0 + days, aso22r)
summary(lmo22r)
lmo26r <- lm(I(plength - intercept) ~ 0 + days, aso26r)
summary(lmo26r)




p1<- ggplot(aso10, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 10C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_abline(intercept = intercept, slope = coef(lmo10), color = "red") +  
  geom_hline(yintercept = 3.55)
p1

p2<- ggplot(aso14r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 14C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_abline(intercept = intercept, slope = coef(lmo14r), color = "red") +  
  geom_hline(yintercept = 3.55)
p2


p3<- ggplot(aso18r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 18C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_abline(intercept = intercept, slope = coef(lmo18r), color = "red") +  
  geom_hline(yintercept = 3.55)
p3


p4<- ggplot(aso22r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 22C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_abline(intercept = intercept, slope = coef(lmo22r), color = "red") +  
  geom_hline(yintercept = 3.55)
p4

p5<- ggplot(aso26r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. orientalis 26C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_abline(intercept = intercept, slope = coef(lmo26r), color = "red") +  
  geom_hline(yintercept = 3.55)
p5



intercept <- 0.47
lmi10r <- lm(I(plength - intercept) ~ 0 + days, asi10r)
summary(lmi10r)
lmi14r <- lm(I(plength - intercept) ~ 0 + days, asi14r)
summary(lmi14r)
lmi18r <- lm(I(plength - intercept) ~ 0 + days, asi18r)
summary(lmi18r)
lmi22r <- lm(I(plength - intercept) ~ 0 + days, asi22r)
summary(lmi22r)
lmi26r <- lm(I(plength - intercept) ~ 0 + days, asi26r)
summary(lmi26r)


p6<- ggplot(asi10r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 10C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_abline(intercept = intercept, slope = coef(lmi10r), color = "red") +  
  geom_hline(yintercept = 4.23)
p6

p7<- ggplot(asi14r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 14C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_abline(intercept = intercept, slope = coef(lmi14r), color = "red") +  
  geom_hline(yintercept = 4.23)
p7


p8<- ggplot(asi18r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 18C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_abline(intercept = intercept, slope = coef(lmi18r), color = "red") +  
  geom_hline(yintercept = 4.23)
p8


p9<- ggplot(asi22r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 22C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_abline(intercept = intercept, slope = coef(lmi22r), color = "red") +  
  geom_hline(yintercept = 4.23)
p9

p10<- ggplot(asi26r, aes(x = days, y = plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("M. intestinalis 26C") +
  scale_color_manual(values=c("steelblue", "orange2")) +
  labs(x = "Days", y = "Length (mm)") +
  ylim(0, 8)+
  xlim(0, 150) +
  geom_abline(intercept = intercept, slope = coef(lmi26r), color = "red") +  
  geom_hline(yintercept = 4.23)
p10

ggarrange(p1, p6, p2, p7, p3, p8, p4, p9, p5, p10,
          ncol = 2, nrow = 5, common.legend = TRUE, legend = "bottom")

#The same with log transformation:
  
  # Specify the desired intercept
  intercept <- log(0.27)  # Log-transformed intercept value
  intercept

# Log-transform the length measurements
aso10$log_plength <- log(aso10$plength)
aso14r$log_plength <- log(aso14r$plength)
aso18r$log_plength <- log(aso18r$plength)
aso22r$log_plength <- log(aso22r$plength)
aso26r$log_plength <- log(aso26r$plength)

lmo10 <- lm(I(log_plength - intercept) ~ 0 + days, aso10)
summary(lmo10)
lmo14r <- lm(I(log_plength - intercept) ~ 0 + days, aso14r)
summary(lmo14r)
lmo18r <- lm(I(log_plength - intercept) ~ 0 + days, aso18r)
summary(lmo18r)
lmo22r <- lm(I(log_plength - intercept) ~ 0 + days, aso22r)
summary(lmo22r)
lmo26r <- lm(I(log_plength - intercept) ~ 0 + days, aso26r)
summary(lmo26r)

p1<- ggplot(aso10, aes(x = days, y = log_plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
    ggtitle(expression(italic("Mytilicola orientalis") ~ "10°C")) +
  scale_color_manual(values=c("#00BFC4", "purple")) +
  labs(x = "Days", y = "Length (log (mm))") +
  ylim(-2, 4)+
  xlim(0, 200) +
  geom_abline(intercept = intercept, slope = coef(lmo10), color = "red") +  
  geom_hline(yintercept = 1.266948)
p1

p2<- ggplot(aso14r, aes(x = days, y = log_plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("14°C ") +
  scale_color_manual(values=c("#00BFC4", "purple ")) +
  labs(x = "Days", y = "Length (log (mm))") +
  ylim(-2, 4)+
  xlim(0, 200) +
  geom_abline(intercept = intercept, slope = coef(lmo14r), color = "red") +  
  geom_hline(yintercept = 1.266948)
p2


p3<- ggplot(aso18r, aes(x = days, y = log_plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("18°C ") +
  scale_color_manual(values=c("#00BFC4", "purple")) +
  labs(x = "Days", y = "Length (log (mm))") +
  ylim(-2, 8)+
  xlim(0, 200) +
  geom_abline(intercept = intercept, slope = coef(lmo18r), color = "red")+
  geom_hline(yintercept = 1.266948)
p3


p4<- ggplot(aso22r, aes(x = days, y = log_plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("22°C ") +
  scale_color_manual(values=c("#00BFC4", "purple")) +
  labs(x = "Days", y = "Length (log (mm))") +
  ylim(-2, 8)+
  xlim(0, 200) +
  geom_abline(intercept = intercept, slope = coef(lmo22r), color = "red") +  
  geom_hline(yintercept = 1.266948)
p4

p5<- ggplot(aso26r, aes(x = days, y = log_plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("26°C") +
  scale_color_manual(values=c("#00BFC4", "purple ")) +
  labs(x = "Days", y = "Length (log (mm))") +
  ylim(-2, 8)+
  xlim(0, 200) +
  geom_abline(intercept = intercept, slope = coef(lmo26r), color = "red") +  
  geom_hline(yintercept = 1.266948)
p5

# Specify the desired intercept
intercept <- log(0.47)  # Log-transformed intercept value
intercept

# Log-transform the length measurements
asi10r$log_plength <- log(asi10r$plength)
asi14r$log_plength <- log(asi14r$plength)
asi18r$log_plength <- log(asi18r$plength)
asi22r$log_plength <- log(asi22r$plength)
asi26r$log_plength <- log(asi26r$plength)
# Fit a linear regression model with the log-transformed intercept
fit <- lm(log_plength ~ 0 + days, data = aso10)


lmi10r <- lm(I(log_plength - intercept) ~ 0 + days, asi10r)
summary(lmi10r)
lmi14r <- lm(I(log_plength - intercept) ~ 0 + days, asi14r)
summary(lmi14r)
lmi18r <- lm(I(log_plength - intercept) ~ 0 + days, asi18r)
summary(lmi18r)
lmi22r <- lm(I(log_plength - intercept) ~ 0 + days, asi22r)
summary(lmi22r)
lmi26r <- lm(I(log_plength - intercept) ~ 0 + days, asi26r)
summary(lmi26r)


log(4.23)

p6<- ggplot(asi10r, aes(x = days, y = log_plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
ggtitle(expression(italic("Mytilicola intestinalis") ~ "10°C")) +  scale_color_manual(values=c("#F8766D", "orange2")) +
  labs(x = "Days", y = "Length (log (mm))") +
  ylim(-2, 8)+
  xlim(0, 200) +
  geom_abline(intercept = intercept, slope = coef(lmi10r), color = "red") +  
  geom_hline(yintercept = 1.442202)
p6

p7<- ggplot(asi14r, aes(x = days, y = log_plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("14°C ") +
  scale_color_manual(values=c("#F8766D", "orange2")) +
  labs(x = "Days", y = "Length (log (mm))") +
  ylim(-2, 8)+
  xlim(0, 200) +
  geom_abline(intercept = intercept, slope = coef(lmi14r), color = "red") +  
  geom_hline(yintercept = 1.442202)
p7


p8<- ggplot(asi18r, aes(x = days, y = log_plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("18°C ") +
  scale_color_manual(values=c("#F8766D", "orange2")) +
  labs(x = "Days", y = "Length (log (mm))") +
  ylim(-2, 8)+
  xlim(0, 200) +
  geom_abline(intercept = intercept, slope = coef(lmi18r), color = "red") +  
  geom_hline(yintercept = 1.442202)
p8


p9<- ggplot(asi22r, aes(x = days, y = log_plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("22°C ") +
  scale_color_manual(values=c("#F8766D", "orange2")) +
  labs(x = "Days",y = "Length (log (mm))") +
  ylim(-2, 8)+
  xlim(0, 200) +
  geom_abline(intercept = intercept, slope = coef(lmi22r), color = "red") +  
  geom_hline(yintercept = 1.442202)
p9

p10<- ggplot(asi26r, aes(x = days, y = log_plength, group=g)) +
  geom_point(aes(shape=g, color=g))+
  ggtitle("26°C ") +
  scale_color_manual(values=c("#F8766D", "orange2")) +
  labs(x = "Days", y = "Length (log (mm))") +
  ylim(-2, 8)+
  xlim(0, 200) +
  geom_abline(intercept = intercept, slope = coef(lmi26r), color = "red") +  
  geom_hline(yintercept = 1.442202)
p10


ggarrange(p6, p1, p7, p2, p8, p3, p9, p4, p10, p5,
          ncol = 2, nrow = 5, common.legend = TRUE, legend = "bottom")
# exported as a tiff with dimensions of 700x1100


#analysing the impact of temperature using linear temperature on the data in each temperature until the first detection of a reproductive size parasite:
# Using rbind to combine data frames
combr <- rbind(aso10, aso14r, aso18r, aso22r, aso26r, asi10r, asi14r, asi18r, asi22r, asi26r)

#including only rows that contained data for the parasite length:
combr2 <- combr[complete.cases(combr$log_plength), ]

model1 <- lmer(plength ~  days *temp *ori + (1 | tank/number), data = combr)
summary(model1)
step1<-step(model1)
step1
model2 <- lmer(plength ~  days *temp *ori + (1 | number), data = combr)
summary(model2)
step2<-step(model2)
step2

model3<-lmer(plength ~ days + temp + ori + (1 | number) + days:ori, data =combr)
summary(model3)