#load packages

library(emmeans)
library(ggplot2)
library(nlme)
library(car)
library(lme4)
library(multcomp)
library(MASS)

#import data
mr <- read.csv("/Users/laurennadler/Desktop/MyRData/Euha-Fupa Metabolic Rate_w AS.csv", header=T)
indiv <- read.csv("/Users/laurennadler/Desktop/MyRData/Euha-Fupa Metabolic Rate_Indiv 2.csv", header=T)
magnitude <- read.csv("/Users/laurennadler/Desktop/MyRData/Euha-Fupa Metabolic Rate_Magnitude.csv", header=T)
enzyme <- read.csv("/Users/laurennadler/Desktop/MyRData/Euha-Fupa Enzyme Activity.csv", header=T)
nka <- read.csv("/Users/Lauren/Desktop/MyRData/Killifish-NKA-WB.csv", header=T)

#Analysis of aerobic metabolic rate

logsmr <- glm(log(SMR)~BM*Infection+Tank, data=mr, na.action = "na.omit")
summary(logsmr)
anova(logsmr, test="F")
plot(logsmr)
qqnorm(logsmr$residuals)
emmeans(logsmr, , specs = pairwise~Infection, adjust = c("none"), type = "response")

logrmr <- glm(log(RMR)~BM*Infection+Metacercaria+Tank, data=mr, na.action = "na.omit")
summary(logsmr)
anova(logrmr, test="F")
plot(logrmr)
qqnorm(logrmr$residuals)
emmeans(logrmr, specs = pairwise~Infection, adjust = c("none"), type = "response")

mmr <- glm(MMR~BM*Infection+Tank, data=mr, na.action = "na.omit")
summary(mmr)
anova(mmr, test="F")
plot(mmr)
qqnorm(mmr$residuals)
emmeans(mmr, specs = pairwise~Infection, adjust = c("none"), type = "response")

as <- glm(AS~BM*Infection+Tank, data=mr, na.action = "na.omit")
summary(as)
anova(as, test="F")
plot(as)
qqnorm(as$residuals)
emmeans(as, specs = pairwise~Infection, adjust = c("none"), type = "response")

#Analysis of the number of new parasite infections

pc <- glm.nb(Postcercaria~BM*Infection*Activity+Tank, data=mr, na.action = "na.omit")
summary(pc)
anova(pc)
r.squaredGLMM(pc)

#Analysis of activity and MRacute

activity.glmer = glmer(Activity ~ Infection*Treatment + (1|Tank:Individual), data=indiv, family=poisson)
summary(activity.glmer)
Anova(activity.glmer)
emmeans(activity.glmer, specs = pairwise~Infection+Treatment, adjust = c("none"), type = "response")
activity1.glmer = glmer(Activity ~ Interaction + (1|Tank:Individual), data=indiv, family=poisson)
summary(glht(activity1.glmer, linfct=mcp(Interaction="Tukey")))
emmeans(activityint.glmer, pairwise~Infection+Treatment, adjust = c("none"))

mracute.lme = lme(MRacute~BM*Infection*Treatment, random =~1|Tank/Individual, data=indiv, na.action = "na.omit")
summary(mracute.lme)
anova(mracute.lme)
plot(mracute.lme)
qqnorm(mracute.lme$residuals)
summary(glht(mracute.lme, linfct=mcp(Treatment="Tukey")))
emmeans(mracute.lme, pairwise~Treatment+Infection, adjust = c("none"))

#Analysis of delta MRacute

absmag.lme = lme(Abs.Mag~BM*Infection*Treatment, random =~1|Tank/Individual, data=magnitude, na.action = "na.omit")
summary(absmag.lme)
anova(absmag.lme)
plot(absmag.lme)
qqnorm(absmag.lme$residuals)

#Metabolic enzyme analysis 

enzyme.lme = lme(Ratio ~ Treatment*Tissue*Assay, random =~1|Tank/Fish, data=enzyme, na.action = "na.omit")
summary(enzyme.lme) 
anova(enzyme.lme)
plot(enzyme.lme)
qqnorm(enzyme.lme$residuals)
enzymeI.lme = lme(Activity ~ Interaction, random =~1|Tank/Fish, data=enzyme, na.action = "na.omit")
summary(glht(enzymeI.lme, linfct=mcp(Interaction="Tukey")))

cs.lme = lme(CS ~ Treatment*Tissue, random =~1|Tank/Fish, data=enzyme, na.action = "na.omit")
summary(cs.lme)
anova(cs.lme)
plot(cs.lme)
qqnorm(cs.lme$residuals)
csI.lme = lme(CS ~ Interaction, random =~1|Tank/Fish, data=enzyme, na.action = "na.omit")
anova(csI.lme)
summary(glht(csI.lme, linfct=mcp(Interaction="Tukey")))

ldh.lme = lme(LDH ~ Treatment*Tissue, random =~1|Tank/Fish, data=enzyme, na.action = "na.omit")
summary(ldh.lme)
anova(ldh.lme)
plot(ldh.lme)
qqnorm(ldh.lme$residuals)
ldhI.lme = lme(LDH ~ Interaction, random =~1|Tank/Fish, data=enzyme, na.action = "na.omit")
anova(ldhI.lme)
library(multcomp)
summary(glht(ldhI.lme, linfct=mcp(Interaction="Tukey")))

#Analysis of relative NKA abundance

nka.glm <- glm(sqrt(NKA)~Infection+Tank, data=nka, na.action = "na.omit")
summary(nka.glm)
anova(nka.glm, test="F")
plot(nka.glm)
qqnorm(nka.glm$residuals)