Life history alterations upon oral and haemocoelic bacterial exposure in the butterfly Melitaea cinxia
Luisa Woestmann, Dimitri Stucki, and Marjo Saastamoinen
1) Introduction
1.1) Assumption and outlier policy
The focus for this project is on the direct outcome that the treatments have on the dependent variables (i.e. the effect). Statistical analyses are performed as consultative measures, and their outcome are not considered the result. As a consequence the assumption and outlier policies are slightly more lenient in the analyses, than for classi hypothesis testing.
Outliers are only removed if the notes indicate a known problem during the experiment stage. Otherwise, far out observations are taken as legit observations and are included in the analysis.
Violations of model assumptions are only corrected in very severe cases. In particular non-normality has been shown to be of lower importance to a model than usually considered (Gelman and Hill, 2006, Data analysis using regression and multilevel/hierarchical models. Cambridge University Press). This relaxation of the normality assumption was done to increase comparability (in particular among genes, where transformation of some genes would bias comparability to the other genes).
2) Preamble
# Load packages
library(psych) # (describeBy)
library(Hmisc) # (errbar)
library(multcomp) # (glht)
library(lme4)
library(lmerTest)
library(dplyr) #(%>%)
library(tidyr)
library(lubridate)
library(knitr)
library(stringr)
library(survival)
library(coxme)
library(influence.ME)
library(VGAM)
library(emmeans)
library(boot)### Standard error
se<-function(x) {sd(x,na.rm=TRUE)/sqrt(length(x[which(x!="NA")]))}
### Confidence interval (eqial variation)
CI<-function(x) {qt(0.975,df=length(x[which(x!="NA")])-1)*(sd(x,na.rm=TRUE)/sqrt(length(x[which(x!="NA")])))}
### Check for overdispersion (@Ben_Bolker's function)
overdisp_fun <- function(model) {
rdf <- df.residual(model)
rp <- residuals(model,type="pearson")
Pearson.chisq <- sum(rp^2)
prat <- Pearson.chisq/rdf
pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}3) Data and preliminary manipulations
All .csv files have to be located in the same directory as this script.
No other external manipulations are required.
4) Read & clean data
4.1) Read
Read all the data.
OE.GE <- read.csv("1a_OralExposure_GeneExpression.csv")
OE.EN <- read.csv("1b_OralExposure_Encapsulation.csv")
OE.RP <- read.csv("1c_OralExposure_Reproduction.csv")
OE.LS <- read.csv("1d_OralExposure_Lifespan.csv")
HE.GE <- read.csv("2a_HaemExposure_GeneExpression.csv")
HE.EN <- read.csv("2b_HaemExposure_Encapsulation.csv")
HE.RP <- read.csv("2c_HaemExposure_Reproduction.csv")
HE.LS <- read.csv("2d_HaemExposure_Lifespan.csv")4.2) Check the structure
Results are suppressed due to the size of the output. Run the script for details.
str(OE.GE)
str(OE.EN)
str(OE.RP)
str(OE.LS)
str(HE.GE)
str(HE.EN)
str(HE.RP)
str(HE.LS)All of these look mostly okay. Some unification and adjustments may be required.
4.3) Clean
Unify the name structure
Unify the column names to feature similar structure.
Replace dots with underscores and collapse double underscores.
colnames(OE.GE)<- tolower(colnames(OE.GE)) %>%
gsub("^(\\w)(\\w+)", "\\U\\1\\L\\2", ., perl = TRUE)
colnames(OE.EN) <- tolower(colnames(OE.EN)) %>%
gsub("^(\\w)(\\w+)", "\\U\\1\\L\\2", ., perl = TRUE)
colnames(OE.RP) <- tolower(colnames(OE.RP)) %>%
gsub("^(\\w)(\\w+)", "\\U\\1\\L\\2", ., perl = TRUE)
colnames(OE.LS) <- tolower(colnames(OE.LS)) %>%
gsub("^(\\w)(\\w+)", "\\U\\1\\L\\2", ., perl = TRUE)
colnames(HE.GE) <- tolower(colnames(HE.GE)) %>%
gsub("^(\\w)(\\w+)", "\\U\\1\\L\\2", ., perl = TRUE)
colnames(HE.EN) <- tolower(colnames(HE.EN)) %>%
gsub("^(\\w)(\\w+)", "\\U\\1\\L\\2", ., perl = TRUE)
colnames(HE.RP) <- tolower(colnames(HE.RP)) %>%
gsub("^(\\w)(\\w+)", "\\U\\1\\L\\2", ., perl = TRUE)
colnames(HE.LS) <- tolower(colnames(HE.LS)) %>%
gsub("^(\\w)(\\w+)", "\\U\\1\\L\\2", ., perl = TRUE)
colnames(HE.RP) <- gsub("\\.","_",colnames(HE.RP))
colnames(HE.RP) <- gsub("__","_",colnames(HE.RP))Convert all dates to date format
## OE.GE: no dates
OE.EN$Birth_date <- as.Date(OE.EN$Birth_date,format="%m/%d/%Y")
OE.RP$Mating_date <- as.Date(OE.RP$Mating_date,format="%m/%d/%Y")
OE.RP$Laying_date <- as.Date(OE.RP$Laying_date,format="%m/%d/%Y")
OE.RP$Hatching_date <- as.Date(OE.RP$Hatching_date,format="%m/%d/%Y")
OE.LS$Birth_date <- as.Date(OE.LS$Birth_date,format="%m/%d/%Y")
OE.LS$Death_date <- as.Date(OE.LS$Death_date,format="%m/%d/%Y")
## HE.GE: no dates
HE.EN$Birth_date <- as.Date(HE.EN$Birth_date,format="%m/%d/%Y")
HE.RP$Mating_date <- as.Date(HE.RP$Mating_date,format="%m/%d/%Y")
HE.RP$Laying_date <- as.Date(HE.RP$Laying_date,format="%m/%d/%Y")
HE.RP$Hatching_date <- as.Date(HE.RP$Hatching_date,format="%m/%d/%Y")
HE.LS$Birth_date <- as.Date(HE.LS$Birth_date,format="%m/%d/%Y")
HE.LS$Death_date <- as.Date(HE.LS$Death_date,format="%m/%d/%Y")Subset to females only
Remove the males, as they are not included in this study. The male data is only available for a subset of the traits studied here. The manuscript focuses on life history traits which are missing for the males.
OE.GE.Fem <- OE.GE[which(OE.GE$Sex=="Female"),]
OE.GE.Fem$Sex <- droplevels(OE.GE.Fem$Sex)
OE.EN.Fem <- OE.EN[which(OE.EN$Sex=="Female"),]
OE.EN.Fem$Sex <- droplevels(OE.EN.Fem$Sex)
OE.LS.Fem <- OE.LS[which(OE.LS$Sex=="Female"),]
OE.LS.Fem$Sex <- droplevels(OE.LS.Fem$Sex)
HE.GE.Fem <- HE.GE[which(HE.GE$Sex=="Female"),]
HE.GE.Fem$Sex <- droplevels(HE.GE.Fem$Sex)
HE.EN.Fem <- HE.EN[which(HE.EN$Sex=="Female"),]
HE.EN.Fem$Sex <- droplevels(HE.EN.Fem$Sex)
HE.LS.Fem <- HE.LS[which(HE.LS$Sex=="female"),]
HE.LS.Fem$Sex <- droplevels(HE.LS.Fem$Sex)4.3.2) Normalize the gene expression data
Normalize Ct-values (NORMA-Gene)
The following code is based on the methods developed in the publication by Heckman et al.
NORMA-Gene: A simple and robust method for qPCR normalization based on target gene data (2011) Lars-Henrik Heckmann, Peter B Sørensen, Paul Henning Krogh and Jesper G Sørensen
OE.GE.Fem.Norm <- data.frame("ID"=NULL,
"FullTreatment"=NULL,"Gene"=NULL,
"NG.dCt"=NULL)
TheTable <- OE.GE.Fem
TheTable$IDs <- with(TheTable, interaction(Larval_id,Family))
TheTable$FullFac <- with(TheTable, interaction(Timepoint,Adult_treatment))
AllTreats <- levels(TheTable$FullFac)
TargetGenes <- levels(TheTable$Target)
for(a in 1:length(AllTreats)){#
SampleTable <- subset(TheTable,FullFac==AllTreats[a])
Larvae <- levels(droplevels(SampleTable$IDs))
N <- dim(SampleTable)[1]/length(TargetGenes)
for(j in 1:N){# loop samples within treatment
z<-0
for(i in 1:length(TargetGenes)){# loop target genes
AveLogi<-mean(log(SampleTable[which(SampleTable$Target==TargetGenes[i]),
"Mean_ct_corr_single_set"]),na.rm=TRUE)
Logij<-log(SampleTable[which(SampleTable$IDs==Larvae[j] &
SampleTable$Target==TargetGenes[i]),
"Mean_ct_corr_single_set"])
z<-z+(AveLogi-ifelse(is.na(Logij),0,Logij))
}
Logaj<-z/sum(!is.na(SampleTable[which(SampleTable$IDs==Larvae[j]),
"Mean_ct_corr_single_set"]))
CtTemp<-cbind(ID=rep(Larvae[j],length(TargetGenes)),
FullTreatment=rep(AllTreats[a],length(TargetGenes)),
Gene=TargetGenes,
NG.dCt=exp(Logaj)*SampleTable[
which(SampleTable$IDs==Larvae[j]),
"Mean_ct_corr_single_set"])
OE.GE.Fem.Norm<-rbind(OE.GE.Fem.Norm,CtTemp,deparse.level=0)
}
}
str(OE.GE.Fem.Norm)## 'data.frame': 600 obs. of 4 variables:
## $ ID : Factor w/ 60 levels "L0316-153.10",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ FullTreatment: Factor w/ 6 levels "24.10","72.10",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Gene : Factor w/ 10 levels "Attacin","beta-1,3-glucan",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ NG.dCt : Factor w/ 600 levels "21.2798879366077",..: 5 10 4 6 2 8 7 9 1 3 ...head(OE.GE.Fem.Norm)## ID FullTreatment Gene NG.dCt
## 1 L0316-153.10 24.10 Attacin 25.114858212461
## 2 L0316-153.10 24.10 beta-1,3-glucan 31.7479503421108
## 3 L0316-153.10 24.10 Histone 23.820807851243
## 4 L0316-153.10 24.10 Lysozyme 25.7865462024402
## 5 L0316-153.10 24.10 mrpL37 22.3922922698447
## 6 L0316-153.10 24.10 Pelle 29.3081800153221# Some housekeeping
OE.GE.Fem.Norm$NG.dCt<-(as.numeric(as.character(OE.GE.Fem.Norm$NG.dCt))*(-1))
OE.GE.Fem.Mat <- spread(OE.GE.Fem.Norm,Gene,NG.dCt)
OE.GE.Fem.Mat$FullTreatment <- gsub("24","t24",OE.GE.Fem.Mat$FullTreatment)
OE.GE.Fem.Mat$FullTreatment <- gsub("72","t72",OE.GE.Fem.Mat$FullTreatment)
OE.GE.Fem.Mat$FullTreatment <- gsub("5","c5",OE.GE.Fem.Mat$FullTreatment)
OE.GE.Fem.Mat$FullTreatment <- gsub("10","c10",OE.GE.Fem.Mat$FullTreatment)
OE.GE.Fem.Mat[,c("Larval_id","Family")] <- str_split_fixed(OE.GE.Fem.Mat$ID, "\\.", 2)
OE.GE.Fem.Mat[,c("Timepoint","Treatment")] <- str_split_fixed(OE.GE.Fem.Mat$FullTreatment, "\\.", 2)
OE.GE.Fem.Mat$Treatment <- factor(OE.GE.Fem.Mat$Treatment,levels=c("control","c5","c10"))
colnames(OE.GE.Fem.Mat) <- gsub("beta-1,3-glucan","B1.3g",colnames(OE.GE.Fem.Mat))
colnames(OE.GE.Fem.Mat) <- gsub("PGRP-LC","PGRP",colnames(OE.GE.Fem.Mat))HE.GE.Fem.Norm <- data.frame("ID"=NULL,
"FullTreatment"=NULL,"Gene"=NULL,
"NG.dCt"=NULL)
TheTable <- HE.GE.Fem
TheTable$IDs <- with(TheTable, interaction(Larval_id,Family))
TheTable$FullFac <- with(TheTable, interaction(Timepoint,Injection_type))
AllTreats <- levels(TheTable$FullFac)
TargetGenes <- levels(TheTable$Target)
for(a in 1:length(AllTreats)){#
SampleTable <- subset(TheTable,FullFac==AllTreats[a])
Larvae <- levels(droplevels(SampleTable$IDs))
N <- dim(SampleTable)[1]/length(TargetGenes)
for(j in 1:N){# loop samples within treatment
z<-0
for(i in 1:length(TargetGenes)){# loop target genes
AveLogi<-mean(log(SampleTable[which(SampleTable$Target==TargetGenes[i]),
"Mean_ct_corr_single_set"]),na.rm=TRUE)
Logij<-log(SampleTable[which(SampleTable$IDs==Larvae[j] &
SampleTable$Target==TargetGenes[i]),
"Mean_ct_corr_single_set"])
z<-z+(AveLogi-ifelse(is.na(Logij),0,Logij))
}
Logaj<-z/sum(!is.na(SampleTable[which(SampleTable$IDs==Larvae[j]),
"Mean_ct_corr_single_set"]))
CtTemp<-cbind(ID=rep(Larvae[j],length(TargetGenes)),
FullTreatment=rep(AllTreats[a],length(TargetGenes)),
Gene=TargetGenes,
NG.dCt=exp(Logaj)*SampleTable[
which(SampleTable$IDs==Larvae[j]),
"Mean_ct_corr_single_set"])
HE.GE.Fem.Norm<-rbind(HE.GE.Fem.Norm,CtTemp,deparse.level=0)
}
}
str(HE.GE.Fem.Norm)## 'data.frame': 590 obs. of 4 variables:
## $ ID : Factor w/ 59 levels "L0116-202.1",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ FullTreatment: Factor w/ 6 levels "24.bacteria",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Gene : Factor w/ 10 levels "Attacin","beta-1,3-glucan",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ NG.dCt : Factor w/ 590 levels "17.8155350226666",..: 1 8 5 6 3 7 9 10 4 2 ...head(HE.GE.Fem.Norm)## ID FullTreatment Gene NG.dCt
## 1 L0116-202.1 24.bacteria Attacin 17.8155350226666
## 2 L0116-202.1 24.bacteria beta-1,3-glucan 26.0467523418135
## 3 L0116-202.1 24.bacteria Histone 23.9514140468485
## 4 L0116-202.1 24.bacteria Lysozyme 24.5717109624664
## 5 L0116-202.1 24.bacteria mrpL37 23.4913507741716
## 6 L0116-202.1 24.bacteria Pelle 25.4249473246108# Some housekeeping
HE.GE.Fem.Norm$NG.dCt<-(as.numeric(as.character(HE.GE.Fem.Norm$NG.dCt))*(-1))
HE.GE.Fem.Mat <- spread(HE.GE.Fem.Norm,Gene,NG.dCt)
HE.GE.Fem.Mat$FullTreatment <- gsub("24","t24",HE.GE.Fem.Mat$FullTreatment)
HE.GE.Fem.Mat$FullTreatment <- gsub("72","t72",HE.GE.Fem.Mat$FullTreatment)
HE.GE.Fem.Mat[,c("Larval_id","Family")] <- str_split_fixed(HE.GE.Fem.Mat$ID, "\\.", 2)
HE.GE.Fem.Mat[,c("Timepoint","Treatment")] <- str_split_fixed(HE.GE.Fem.Mat$FullTreatment, "\\.", 2)
HE.GE.Fem.Mat$Treatment <- factor(HE.GE.Fem.Mat$Treatment,levels=c("naive","PBS","bacteria"))
colnames(HE.GE.Fem.Mat) <- gsub("beta-1,3-glucan","B1.3g",colnames(HE.GE.Fem.Mat))
colnames(HE.GE.Fem.Mat) <- gsub("PGRP-LC","PGRP",colnames(HE.GE.Fem.Mat))Adjust treatment levels
The following blocks specify the names of treatments and the sequence of the factor levels. Furthermore the weight of the eggs is adjusted to ng, and the time until laying is calculated.
In addition, a new data frame is created containing the amount of clutches each female produced in its lifetime, and lifetime egg production is calculated.
Variables of oral and haemocoelic exposure are adjusted to feature the same naming.
OE.EN.Fem$Treatment <- gsub("5mg/mL","c5",OE.EN.Fem$Treatment)
OE.EN.Fem$Treatment <- gsub("10mg/mL","c10",OE.EN.Fem$Treatment)
OE.EN.Fem$Treatment <- factor(OE.EN.Fem$Treatment,levels=c("Control","c5","c10"))
OE.RP$Treatment <- gsub("bacteria","Bacteria",OE.RP$Treatment)
OE.RP$Treatment <- gsub("control","Control",OE.RP$Treatment)
OE.RP$Treatment <- factor(OE.RP$Treatment,levels=c("Control","Bacteria"))
OE.RP$Weight_per_egg <- OE.RP$Weight_per_egg*1000
OE.RP$Time_until_laying <- as.numeric(OE.RP$Laying_date-OE.RP$Mating_date)
## Clutches per female
OE.CPF <- aggregate(OE.RP$Clutch_rank,by=list(OE.RP$Female_id,OE.RP$Treatment),max)
## Lifetime egg production
OE.LEP <- aggregate(OE.RP$Number_eggs,by=list(OE.RP$Female_id,OE.RP$Treatment),sum)
OE.LS.Fem$Feeding_treatment <- gsub("5","c5",OE.LS.Fem$Feeding_treatment)
OE.LS.Fem$Feeding_treatment <- gsub("10","c10",OE.LS.Fem$Feeding_treatment)
OE.LS.Fem$Feeding_treatment <- factor(OE.LS.Fem$Feeding_treatment,levels=c("Control","c5","c10"))
HE.EN.Fem$Treatment <- HE.EN.Fem$Treatment %>%
trimws() %>%
{factor(.,levels=c("Control","PBS","Bacterial"))}
colnames(HE.RP) <- gsub("Number_of_eggs","Number_eggs",colnames(HE.RP))
HE.RP[order(HE.RP$Female_id,HE.RP$Clutch),"Clutch_rank"] <-
unlist(tapply(HE.RP$Clutch,HE.RP$Female_id,rank))
HE.RP$Treatment <- HE.RP$Treatment %>%
{gsub("control","Control",.)} %>%
{factor(.,levels=c("Control","PBS","Bacteria"))}
HE.RP$Time_until_laying <- as.numeric(HE.RP$Laying_date-HE.RP$Mating_date)
HE.RP$Weight_per_egg <- as.numeric(as.character(
HE.RP$Egg_weight_mg_weight_before_counting))/HE.RP$Number_eggs
HE.RP$Hatching_percentage <- (HE.RP$Number_hatched/HE.RP$Number_eggs)*100
## Clutches per female
HE.CPF <- aggregate(HE.RP$Clutch_rank,by=list(HE.RP$Female_id,HE.RP$Treatment),max)
## Lifetime egg production
HE.LEP <- aggregate(HE.RP$Number_eggs,by=list(HE.RP$Female_id,HE.RP$Treatment),sum)5) Analysis: Oral exposure
5.1) Gene expression
5.1.1) Assess PCA
## Target gene columns: 3,4,6,8,9,10,12
cortest.bartlett(OE.GE.Fem.Mat[,c(3,4,6,8,9,10,12)])## $chisq
## [1] 140.6298
##
## $p.value
## [1] 1.040403e-19
##
## $df
## [1] 21det(cor(as.matrix(OE.GE.Fem.Mat[,c(3,4,6,8,9,10,12)])))## [1] 0.08056082KMO(OE.GE.Fem.Mat[,c(3,4,6,8,9,10,12)])## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = OE.GE.Fem.Mat[, c(3, 4, 6, 8, 9, 10, 12)])
## Overall MSA = 0.46
## MSA for each item =
## Attacin B1.3g Lysozyme Pelle PGRP proPO Serpin
## 0.44 0.38 0.34 0.43 0.60 0.49 0.65Colinearity is okay, but the sampling adequacy is not. Removal of genes may help, but with so few genes the remaining PCA would not tell much more than a gene-by-gene analysis.
I run the gene-by-gene analysis directly.
5.1.2) Gene by gene; mixed linear models
Only the model specification stepp is shown here. These models are primarily used for the interaction effects, as the pairwise comparisons are of primary interest.
OE.model.adhoc.Attacin.Full <- lmer(Attacin~(Treatment+Timepoint)^2+(1|Family),
data=OE.GE.Fem.Mat)
OE.model.adhoc.B1.3g.Full <- lmer(B1.3g~(Treatment+Timepoint)^2+(1|Family),
data=OE.GE.Fem.Mat)
OE.model.adhoc.Lysozyme.Full <- lmer(Lysozyme~(Treatment+Timepoint)^2+(1|Family),
data=OE.GE.Fem.Mat)
OE.model.adhoc.Pelle.Full <- lmer(Pelle~(Treatment+Timepoint)^2+(1|Family),
data=OE.GE.Fem.Mat)
OE.model.adhoc.PGRP.Full <- lmer(PGRP~(Treatment+Timepoint)^2+(1|Family),
data=OE.GE.Fem.Mat)
OE.model.adhoc.proPO.Full <- lmer(proPO~(Treatment+Timepoint)^2+(1|Family),
data=OE.GE.Fem.Mat)
OE.model.adhoc.Serpin.Full <- lmer(Serpin~(Treatment+Timepoint)^2+(1|Family),
data=OE.GE.Fem.Mat)Adjust p-values
Adjust all p-values for the comparison of seven genes using false discovery rates.
pvals.adho <- rbind(Attacin=coef(summary(OE.model.adhoc.Attacin.Full))[,5],
B1.3g=coef(summary(OE.model.adhoc.B1.3g.Full))[,5],
Lysozyme=coef(summary(OE.model.adhoc.Lysozyme.Full))[,5],
Pelle=coef(summary(OE.model.adhoc.Pelle.Full))[,5],
PGRP=coef(summary(OE.model.adhoc.PGRP.Full))[,5],
proPO=coef(summary(OE.model.adhoc.proPO.Full))[,5],
Serpin=coef(summary(OE.model.adhoc.Serpin.Full))[,5])
qvals.adho<-apply(pvals.adho,2,p.adjust,method="fdr")Validate the models
1) Attacin:
## Attacin
## Linearity
plot(resid(OE.model.adhoc.Attacin.Full),fitted(OE.model.adhoc.Attacin.Full))
## Homogeneity
plot(OE.model.adhoc.Attacin.Full)
## Normality
qqnorm(residuals(OE.model.adhoc.Attacin.Full))
## Outlier
plot(hatvalues(OE.model.adhoc.Attacin.Full)~residuals(OE.model.adhoc.Attacin.Full))
- Linearity is fine.
- Homogeneity is fine.
- The qq-line looks a bit drunk, but it looks okay for analysis.
- No notable outliers.
2) B1.3g
## B1.3g
## Linearity
plot(resid(OE.model.adhoc.B1.3g.Full),fitted(OE.model.adhoc.B1.3g.Full))
## Homogeneity
plot(OE.model.adhoc.B1.3g.Full)
## Normality
qqnorm(residuals(OE.model.adhoc.B1.3g.Full))
## Outlier
plot(hatvalues(OE.model.adhoc.B1.3g.Full)~residuals(OE.model.adhoc.B1.3g.Full))
- Linearity is fine.
- Homogeneity is fine.
- The qq-line has a slight drop at its beginning. But I think this is okay, as linearity is not the most important assumption.
- No notable outliers.
3) Lysozyme
## Lysozyme
## Linearity
plot(resid(OE.model.adhoc.Lysozyme.Full),fitted(OE.model.adhoc.Lysozyme.Full))
## Homogeneity
plot(OE.model.adhoc.Lysozyme.Full)
## Normality
qqnorm(residuals(OE.model.adhoc.Lysozyme.Full))
## Outlier
plot(hatvalues(OE.model.adhoc.Lysozyme.Full)~residuals(OE.model.adhoc.Lysozyme.Full))
- Linearity is fine.
- Homogeneity could be better.
- The qq-line has a slight bending. But all other assumptions are okay and normality is the least influential assumption in linear models. I’ll keep it.
- No notable outliers.
4) Pelle
## Pelle
## Linearity
plot(resid(OE.model.adhoc.Pelle.Full),fitted(OE.model.adhoc.Pelle.Full))
## Homogeneity
plot(OE.model.adhoc.Pelle.Full)
## Normality
qqnorm(residuals(OE.model.adhoc.Pelle.Full))
## Outlier
plot(hatvalues(OE.model.adhoc.Pelle.Full)~residuals(OE.model.adhoc.Pelle.Full))
- Linearity is fine.
- Homogeneity is fine.
- The qq-line has strong bends at the beginning and end. The rest of the line looks normal, and normality is not a very important assumption.
- No notable outliers.
5) PGRP
## PGRP
## Linearity
plot(resid(OE.model.adhoc.PGRP.Full),fitted(OE.model.adhoc.PGRP.Full))
## Homogeneity
plot(OE.model.adhoc.PGRP.Full)
## Normality
qqnorm(residuals(OE.model.adhoc.PGRP.Full))
## Outlier
plot(hatvalues(OE.model.adhoc.PGRP.Full)~residuals(OE.model.adhoc.PGRP.Full))
- Linearity looks just okay. Outlier?
- Homogeneity looks okay. Outlier?
- The qq-line looks a bit drunk, but okay. Outlier?
- Two possible outliers.
which(hatvalues(OE.model.adhoc.PGRP.Full)>0.5)## 26 27
## 26 27OE.GE.Fem.Mat[26:27,]## ID FullTreatment Attacin B1.3g Histone Lysozyme
## 26 L0316-197.59 t24.c5 -22.21241 -28.95806 -22.84675 -25.68251
## 27 L0316-352.68 t24.c5 -18.88389 -28.00184 -23.84010 -23.98594
## mrpL37 Pelle PGRP proPO S28rpS24 Serpin Larval_id
## 26 -23.80555 -27.20844 -28.32669 -29.11996 -22.90269 -23.53343 L0316-197
## 27 -22.25695 -31.06983 -32.90216 -29.78290 -21.83648 -24.29103 L0316-352
## Family Timepoint Treatment
## 26 59 t24 c5
## 27 68 t24 c5OE.model.adhoc.PGRP.Full.NoOut <- lmer(PGRP~(Treatment+Timepoint)^2+(1|Family),
data=OE.GE.Fem.Mat[-c(26,27),])
## Linearity
plot(resid(OE.model.adhoc.PGRP.Full.NoOut),fitted(OE.model.adhoc.PGRP.Full.NoOut))
## Homogeneity
plot(OE.model.adhoc.PGRP.Full.NoOut)
## Normality
qqnorm(residuals(OE.model.adhoc.PGRP.Full.NoOut))
## Outlier
plot(hatvalues(OE.model.adhoc.PGRP.Full.NoOut)~residuals(OE.model.adhoc.PGRP.Full.NoOut))
Removal of the possible outliers fixes the graphs, but there is no information in the notes on possible reasons for these outliers. So technically they are valid observations. I checked the influence on the final model and there is a slight increase in significance in the effect of 5mg/ml bacteria.
Given that there is no information on why to remove these, I’ll keep the two extreme values in.
6) ProPO
## Linearity
plot(resid(OE.model.adhoc.proPO.Full),fitted(OE.model.adhoc.proPO.Full))
## Homogeneity
plot(OE.model.adhoc.proPO.Full)
## Normality
qqnorm(residuals(OE.model.adhoc.proPO.Full))
## Outlier
plot(hatvalues(OE.model.adhoc.proPO.Full)~residuals(OE.model.adhoc.proPO.Full))
- Linearity is fine.
- Homogeneity is fine.
- The qq-line has a slight drop at its beginning. But I think this is just okay.
- No notable outliers.
7) Serpin
## Serpin
## Linearity
plot(resid(OE.model.adhoc.Serpin.Full),fitted(OE.model.adhoc.Serpin.Full))
## Homogeneity
plot(OE.model.adhoc.Serpin.Full)
## Normality
qqnorm(residuals(OE.model.adhoc.Serpin.Full))
## Outlier
plot(hatvalues(OE.model.adhoc.Serpin.Full)~residuals(OE.model.adhoc.Serpin.Full))
- Linearity is fine.
- Homogeneity is fine.
- The qq-line looks messed. But I think this is just okay.
- No notable outliers.
5.1.3) Gene by gene; pairwise comparisons
Create comparisons for a subset of pre-specified comparisons.
Adjust the p-values for an analysis of seven genes, using the fdr method.
comps<-c(
"t24.control - t24.c5 = 0","t24.control - t24.c10 = 0",
"t72.control - t72.c5 = 0","t72.control - t72.c10 = 0",
"t24.control - t72.control = 0","t24.c5 - t72.c5 = 0",
"t24.c10 - t72.c10 = 0"
)
OE.model.posthoc.Attacin <- lmer(Attacin~FullTreatment+(1|Family),data=OE.GE.Fem.Mat)
summ.Att <- summary(glht(OE.model.posthoc.Attacin,linfct=mcp(FullTreatment=comps)),test=adjusted(type="none"))
OE.model.posthoc.B1.3g <- lmer(B1.3g~FullTreatment+(1|Family),data=OE.GE.Fem.Mat)
summ.B13 <- summary(glht(OE.model.posthoc.B1.3g,linfct=mcp(FullTreatment=comps)),test=adjusted(type="none"))
OE.model.posthoc.Lysozyme <- lmer(Lysozyme~FullTreatment+(1|Family),data=OE.GE.Fem.Mat)
summ.Lys <- summary(glht(OE.model.posthoc.Lysozyme,linfct=mcp(FullTreatment=comps)),test=adjusted(type="none"))
OE.model.posthoc.Pelle <- lmer(Pelle~FullTreatment+(1|Family),data=OE.GE.Fem.Mat)
summ.Pel <- summary(glht(OE.model.posthoc.Pelle,linfct=mcp(FullTreatment=comps)),test=adjusted(type="none"))
OE.model.posthoc.PGRP <- lmer(PGRP~FullTreatment+(1|Family),data=OE.GE.Fem.Mat)
summ.PGR <- summary(glht(OE.model.posthoc.PGRP,linfct=mcp(FullTreatment=comps)),test=adjusted(type="none"))
OE.model.posthoc.proPO <- lmer(proPO~FullTreatment+(1|Family),data=OE.GE.Fem.Mat)
summ.pro <- summary(glht(OE.model.posthoc.proPO,linfct=mcp(FullTreatment=comps)),test=adjusted(type="none"))
OE.model.posthoc.Serpin <- lmer(Serpin~FullTreatment+(1|Family),data=OE.GE.Fem.Mat)
summ.Ser <- summary(glht(OE.model.posthoc.Serpin,linfct=mcp(FullTreatment=comps)),test=adjusted(type="none"))
pvals.poho <- rbind(Attacin=summ.Att$test$pvalues,
B1.3g=summ.B13$test$pvalues,
Lysozyme=summ.Lys$test$pvalues,
Pelle=summ.Pel$test$pvalues,
PGRP=summ.PGR$test$pvalues,
proPO=summ.pro$test$pvalues,
Serpin=summ.Ser$test$pvalues)
qvals.poho<-t(apply(pvals.poho,2,p.adjust,method="fdr"))Create table
Create a list with the results and write it to a table.
OE.GE.Fem.GxG.PH.Results <- list(Attacin=
summ.Att$test[c("coefficients","sigma","tstat","pvalues")] %>%
{cbind(estimate=
{paste(format(round(.$coefficients,digits=2),nsmall=2),
format(round(.$sigma,digits=2),nsmall=2),sep="±")},
tvalue=format(round(.$tstat,digits=4),nsmall=4),
pvalue=format(round(.$pvalues,digits=4),nsmall=4),
qvalue=format(round(qvals.poho[,"Attacin"],digits=4),nsmall=4))},
B1.3g=
summ.B13$test[c("coefficients","sigma","tstat","pvalues")] %>%
{cbind(estimate=
{paste(format(round(.$coefficients,digits=2),nsmall=2),
format(round(.$sigma,digits=2),nsmall=2),sep="±")},
tvalue=format(round(.$tstat,digits=4),nsmall=4),
pvalue=format(round(.$pvalues,digits=4),nsmall=4),
qvalue=format(round(qvals.poho[,"B1.3g"],digits=4),nsmall=4))},
Lysozyme=
summ.Lys$test[c("coefficients","sigma","tstat","pvalues")] %>%
{cbind(estimate=
{paste(format(round(.$coefficients,digits=2),nsmall=2),
format(round(.$sigma,digits=2),nsmall=2),sep="±")},
tvalue=format(round(.$tstat,digits=4),nsmall=4),
pvalue=format(round(.$pvalues,digits=4),nsmall=4),
qvalue=format(round(qvals.poho[,"Lysozyme"],digits=4),nsmall=4))},
Pelle=
summ.Pel$test[c("coefficients","sigma","tstat","pvalues")] %>%
{cbind(estimate=
{paste(format(round(.$coefficients,digits=2),nsmall=2),
format(round(.$sigma,digits=2),nsmall=2),sep="±")},
tvalue=format(round(.$tstat,digits=4),nsmall=4),
pvalue=format(round(.$pvalues,digits=4),nsmall=4),
qvalue=format(round(qvals.poho[,"Pelle"],digits=4),nsmall=4))},
PGRP=
summ.PGR$test[c("coefficients","sigma","tstat","pvalues")] %>%
{cbind(estimate=
{paste(format(round(.$coefficients,digits=2),nsmall=2),
format(round(.$sigma,digits=2),nsmall=2),sep="±")},
tvalue=format(round(.$tstat,digits=4),nsmall=4),
pvalue=format(round(.$pvalues,digits=4),nsmall=4),
qvalue=format(round(qvals.poho[,"PGRP"],digits=4),nsmall=4))},
proPO=
summ.pro$test[c("coefficients","sigma","tstat","pvalues")] %>%
{cbind(estimate=
{paste(format(round(.$coefficients,digits=2),nsmall=2),
format(round(.$sigma,digits=2),nsmall=2),sep="±")},
tvalue=format(round(.$tstat,digits=4),nsmall=4),
pvalue=format(round(.$pvalues,digits=4),nsmall=4),
qvalue=format(round(qvals.poho[,"proPO"],digits=4),nsmall=4))},
Serpin=
summ.Ser$test[c("coefficients","sigma","tstat","pvalues")] %>%
{cbind(estimate=
{paste(format(round(.$coefficients,digits=2),nsmall=2),
format(round(.$sigma,digits=2),nsmall=2),sep="±")},
tvalue=format(round(.$tstat,digits=4),nsmall=4),
pvalue=format(round(.$pvalues,digits=4),nsmall=4),
qvalue=format(round(qvals.poho[,"Serpin"],digits=4),nsmall=4))}
)
OE.GE.Fem.GxG.PH.Results## $Attacin
## estimate tvalue pvalue qvalue
## t24.control - t24.c5 "-2.06±0.94" "-2.1928" "0.0283" "0.0991"
## t24.control - t24.c10 "-0.32±0.89" "-0.3528" "0.7242" "0.9382"
## t72.control - t72.c5 "-0.32±0.89" "-0.3640" "0.7159" "0.9584"
## t72.control - t72.c10 " 0.43±0.99" " 0.4314" "0.6662" "0.7772"
## t24.control - t72.control "-0.64±0.93" "-0.6891" "0.4908" "0.6387"
## t24.c5 - t72.c5 " 1.10±0.91" " 1.2131" "0.2251" "0.3919"
## t24.c10 - t72.c10 " 0.11±0.98" " 0.1083" "0.9138" "0.9932"
##
## $B1.3g
## estimate tvalue pvalue qvalue
## t24.control - t24.c5 "-1.13±0.47" "-2.3845" "0.0171" "0.0991"
## t24.control - t24.c10 " 0.59±0.46" " 1.2840" "0.1991" "0.3485"
## t72.control - t72.c5 " 0.32±0.46" " 0.7012" "0.4832" "0.9584"
## t72.control - t72.c10 " 1.27±0.50" " 2.5280" "0.0115" "0.0698"
## t24.control - t72.control "-0.48±0.47" "-1.0193" "0.3081" "0.6387"
## t24.c5 - t72.c5 " 0.97±0.46" " 2.0985" "0.0359" "0.2510"
## t24.c10 - t72.c10 " 0.19±0.49" " 0.3922" "0.6949" "0.9729"
##
## $Lysozyme
## estimate tvalue pvalue qvalue
## t24.control - t24.c5 "-0.51±0.39" "-1.3029" "0.1926" "0.3371"
## t24.control - t24.c10 " 0.03±0.36" " 0.0876" "0.9302" "0.9382"
## t72.control - t72.c5 "-0.15±0.36" "-0.4119" "0.6804" "0.9584"
## t72.control - t72.c10 "-0.07±0.41" "-0.1618" "0.8714" "0.8714"
## t24.control - t72.control " 0.10±0.38" " 0.2681" "0.7886" "0.7886"
## t24.c5 - t72.c5 " 0.46±0.37" " 1.2364" "0.2163" "0.3919"
## t24.c10 - t72.c10 " 0.00±0.40" " 0.0085" "0.9932" "0.9932"
##
## $Pelle
## estimate tvalue pvalue qvalue
## t24.control - t24.c5 " 0.28±0.48" " 0.5816" "0.5608" "0.6543"
## t24.control - t24.c10 " 1.08±0.45" " 2.4167" "0.0157" "0.0548"
## t72.control - t72.c5 " 0.48±0.44" " 1.0757" "0.2821" "0.9584"
## t72.control - t72.c10 " 0.49±0.51" " 0.9542" "0.3400" "0.5950"
## t24.control - t72.control " 0.28±0.47" " 0.6016" "0.5474" "0.6387"
## t24.c5 - t72.c5 " 0.48±0.46" " 1.0354" "0.3005" "0.3919"
## t24.c10 - t72.c10 "-0.32±0.50" "-0.6263" "0.5311" "0.9729"
##
## $PGRP
## estimate tvalue pvalue qvalue
## t24.control - t24.c5 "0.30±0.39" "0.7714" "0.4404" "0.6166"
## t24.control - t24.c10 "0.95±0.33" "2.8728" "0.0041" "0.0285"
## t72.control - t72.c5 "0.40±0.32" "1.2275" "0.2196" "0.9584"
## t72.control - t72.c10 "0.92±0.39" "2.3274" "0.0199" "0.0698"
## t24.control - t72.control "0.25±0.35" "0.7076" "0.4792" "0.6387"
## t24.c5 - t72.c5 "0.35±0.36" "0.9622" "0.3359" "0.3919"
## t24.c10 - t72.c10 "0.21±0.39" "0.5436" "0.5867" "0.9729"
##
## $proPO
## estimate tvalue pvalue qvalue
## t24.control - t24.c5 "-0.80±0.60" "-1.3278" "0.1842" "0.3371"
## t24.control - t24.c10 " 0.05±0.59" " 0.0776" "0.9382" "0.9382"
## t72.control - t72.c5 " 0.03±0.59" " 0.0522" "0.9584" "0.9584"
## t72.control - t72.c10 " 0.95±0.64" " 1.4935" "0.1353" "0.3157"
## t24.control - t72.control "-0.49±0.60" "-0.8204" "0.4120" "0.6387"
## t24.c5 - t72.c5 " 0.33±0.59" " 0.5715" "0.5677" "0.5677"
## t24.c10 - t72.c10 " 0.41±0.62" " 0.6621" "0.5079" "0.9729"
##
## $Serpin
## estimate tvalue pvalue qvalue
## t24.control - t24.c5 " 0.06±0.31" " 0.2083" "0.8350" "0.8350"
## t24.control - t24.c10 " 0.42±0.30" " 1.3807" "0.1674" "0.3485"
## t72.control - t72.c5 "-0.03±0.30" "-0.0852" "0.9321" "0.9584"
## t72.control - t72.c10 " 0.16±0.33" " 0.4916" "0.6230" "0.7772"
## t24.control - t72.control "-0.25±0.31" "-0.8209" "0.4117" "0.6387"
## t24.c5 - t72.c5 "-0.35±0.30" "-1.1386" "0.2549" "0.3919"
## t24.c10 - t72.c10 "-0.51±0.32" "-1.5866" "0.1126" "0.7883"#lapply(OE.GE.Fem.GxG.PH.Results, function(x) write.table( data.frame(x),
# "ResultsV3/OE.GE.Fem.GxG.PH.Results.csv", append=T, sep=',' ))5.1.4) Gene by gene; Plots
head(OE.GE.Fem.Mat)## ID FullTreatment Attacin B1.3g Histone Lysozyme
## 1 L0316-153.10 t24.c10 -25.11486 -31.74795 -23.82081 -25.78655
## 2 L0316-222.12 t24.c10 -25.95286 -30.74498 -23.35107 -25.66155
## 3 L0316-309.26 t24.c10 -23.25323 -31.29952 -23.18202 -26.54880
## 4 L0316-201.49 t24.c10 -25.30387 -30.80892 -23.64144 -25.89706
## 5 L0316-149.74 t24.c10 -28.64305 -30.54497 -22.99717 -25.77958
## 6 L0316-132.80 t24.c10 -25.27256 -31.97243 -23.63651 -25.32455
## mrpL37 Pelle PGRP proPO S28rpS24 Serpin Larval_id
## 1 -22.39229 -29.30818 -29.17201 -31.43168 -21.27989 -23.17691 L0316-153
## 2 -23.51957 -28.86817 -28.50096 -29.06331 -23.21788 -23.28185 L0316-222
## 3 -23.30021 -27.80881 -29.66233 -28.74173 -22.82930 -25.73194 L0316-309
## 4 -23.82217 -27.69176 -28.84018 -28.74969 -23.68400 -23.50914 L0316-201
## 5 -22.44880 -28.58916 -28.54557 -29.96492 -22.14189 -23.05537 L0316-149
## 6 -22.36366 -28.28929 -29.24741 -32.06730 -22.02840 -22.99295 L0316-132
## Family Timepoint Treatment
## 1 10 t24 c10
## 2 12 t24 c10
## 3 26 t24 c10
## 4 49 t24 c10
## 5 74 t24 c10
## 6 80 t24 c10OE.GE.Fem.Res.means <- aggregate(OE.GE.Fem.Mat[,c(3,4,6,8,9,10,12)],by=list(OE.GE.Fem.Mat$Treatment,
OE.GE.Fem.Mat$Timepoint),FUN=mean,na.rm=TRUE)
OE.GE.Fem.Res.CIs <- aggregate(OE.GE.Fem.Mat[,c(3,4,6,8,9,10,12)],by=list(OE.GE.Fem.Mat$Treatment,
OE.GE.Fem.Mat$Timepoint),FUN=CI)
xpos <- seq(0,2,1)+rep(c(0,4),each=3)+0.5
#svg("ResultsV3/OE.GE.Fem.Plot.svg",width=20,height=10)
par(mfrow=c(3,3))
par(mar=c(2,6,1,1))
for(i in 3:9){
ylimi <- round(range(c(floor(OE.GE.Fem.Res.means[,i]+OE.GE.Fem.Res.CIs[,i]*-1),ceiling(OE.GE.Fem.Res.means[,i]+OE.GE.Fem.Res.CIs[,i]))))
plot(NA,xlim=c(0,7),ylim=ylimi,axes=FALSE,xlab="",cex.lab=5,ylab="")
axis(side=1,at=c(-2,8),labels=c("",""),font=2)
axis(side=2,at=seq(ylimi[1]-2,ylimi[2]+1,1),labels=seq(ylimi[1]-2,ylimi[2]+1,1),
las=2,font=2,cex.axis=3)
axis(side=3,at=c(-2,8))
axis(side=4,at=c(ylimi[1]-2,ylimi[2]+1))
errbar(x=xpos,y=OE.GE.Fem.Res.means[,i],yplus=OE.GE.Fem.Res.means[,i]+OE.GE.Fem.Res.CIs[,i],yminus=OE.GE.Fem.Res.means[,i]-OE.GE.Fem.Res.CIs[,i],lwd=4,pch="",add=TRUE)
points(x=xpos,y=OE.GE.Fem.Res.means[,i],pch=c(21,22,23),col=c("black"),cex=4,lwd=4,bg=rep(c("black","gray"),each=3))
text(0,ylimi[1],labels=colnames(OE.GE.Fem.Res.means)[i],cex=2,adj=0)
}
#dev.off()
Version two; add all points
head(OE.GE.Fem.Mat)## ID FullTreatment Attacin B1.3g Histone Lysozyme
## 1 L0316-153.10 t24.c10 -25.11486 -31.74795 -23.82081 -25.78655
## 2 L0316-222.12 t24.c10 -25.95286 -30.74498 -23.35107 -25.66155
## 3 L0316-309.26 t24.c10 -23.25323 -31.29952 -23.18202 -26.54880
## 4 L0316-201.49 t24.c10 -25.30387 -30.80892 -23.64144 -25.89706
## 5 L0316-149.74 t24.c10 -28.64305 -30.54497 -22.99717 -25.77958
## 6 L0316-132.80 t24.c10 -25.27256 -31.97243 -23.63651 -25.32455
## mrpL37 Pelle PGRP proPO S28rpS24 Serpin Larval_id
## 1 -22.39229 -29.30818 -29.17201 -31.43168 -21.27989 -23.17691 L0316-153
## 2 -23.51957 -28.86817 -28.50096 -29.06331 -23.21788 -23.28185 L0316-222
## 3 -23.30021 -27.80881 -29.66233 -28.74173 -22.82930 -25.73194 L0316-309
## 4 -23.82217 -27.69176 -28.84018 -28.74969 -23.68400 -23.50914 L0316-201
## 5 -22.44880 -28.58916 -28.54557 -29.96492 -22.14189 -23.05537 L0316-149
## 6 -22.36366 -28.28929 -29.24741 -32.06730 -22.02840 -22.99295 L0316-132
## Family Timepoint Treatment
## 1 10 t24 c10
## 2 12 t24 c10
## 3 26 t24 c10
## 4 49 t24 c10
## 5 74 t24 c10
## 6 80 t24 c10GeneCols <-c(3,4,6,8,9,10,12)
OE.GE.Fem.Res.means <- aggregate(OE.GE.Fem.Mat[,c(3,4,6,8,9,10,12)],by=list(OE.GE.Fem.Mat$Treatment,
OE.GE.Fem.Mat$Timepoint),FUN=mean,na.rm=TRUE)
OE.GE.Fem.Res.CIs <- aggregate(OE.GE.Fem.Mat[,c(3,4,6,8,9,10,12)],by=list(OE.GE.Fem.Mat$Treatment,
OE.GE.Fem.Mat$Timepoint),FUN=CI)
xpos <- seq(0,2,1)+rep(c(0,4),each=3)+0.5
#svg("ResultsV3/OE.GE.Fem.Plot_Points.svg",width=20,height=10)
par(mfrow=c(3,3))
par(mar=c(2,6,1,1))
for(i in 3:9){
ylimi <- round(range(c(floor(OE.GE.Fem.Mat[,GeneCols[i-2]]),ceiling(OE.GE.Fem.Mat[,GeneCols[i-2]]))))
plot(NA,xlim=c(0,7),ylim=ylimi,axes=FALSE,xlab="",cex.lab=5,ylab="")
axis(side=1,at=c(-2,8),labels=c("",""),font=2)
axis(side=2,at=seq(ylimi[1]-2,ylimi[2]+1,1),labels=seq(ylimi[1]-2,ylimi[2]+1,1),
las=2,font=2,cex.axis=3)
axis(side=3,at=c(-2,8))
axis(side=4,at=c(ylimi[1]-2,ylimi[2]+1))
errbar(x=xpos,y=OE.GE.Fem.Res.means[,i],yplus=OE.GE.Fem.Res.means[,i]+OE.GE.Fem.Res.CIs[,i],yminus=OE.GE.Fem.Res.means[,i]-OE.GE.Fem.Res.CIs[,i],lwd=4,pch="",add=TRUE)
points(x=xpos,y=OE.GE.Fem.Res.means[,i],pch=c(21,22,23),col=c("black"),cex=4,lwd=4,bg=rep(c("black","gray"),each=3))
text(0,ylimi[1],labels=colnames(OE.GE.Fem.Res.means)[i],cex=2,adj=0)
points(x=rep(xpos,tapply(OE.GE.Fem.Mat[,GeneCols[i-2]],list(OE.GE.Fem.Mat$Treatment,OE.GE.Fem.Mat$Timepoint),length)),OE.GE.Fem.Mat[order(OE.GE.Fem.Mat$Timepoint,OE.GE.Fem.Mat$Treatment),GeneCols[i-2]],pch=4,col="gray32",cex=1.8,lwd=1.3)
}
#dev.off()
5.2) Encapsulation
5.2.1) Analysis
OE.model.EN.Full <- lmer(Corr_encaps~Treatment+(1|Family),data=OE.EN.Fem)
## Linearity
plot(resid(OE.model.EN.Full),fitted(OE.model.EN.Full))
## Homogeneity
plot(OE.model.EN.Full)
## Normality
qqnorm(residuals(OE.model.EN.Full))
## Influence
infl <- influence(OE.model.EN.Full, obs = TRUE)
plot(infl, which = "cook")
- Linearity is fine.
- Homogeneity is fine.
- The qq-line has a slight bend. But I think this is okay.
- No notable outliers.
Get the results
Display the model summary, and calculate the means of the gray values for each treatment.
As we are not interested in the comparison between the concentrations, the model summary provides sufficient information.
summary(OE.model.EN.Full)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Corr_encaps ~ Treatment + (1 | Family)
## Data: OE.EN.Fem
##
## REML criterion at convergence: 788.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8333 -0.5524 0.1255 0.6665 1.7232
##
## Random effects:
## Groups Name Variance Std.Dev.
## Family (Intercept) 0.0 0.00
## Residual 366.2 19.14
## Number of obs: 92, groups: Family, 19
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.03689 3.38276 89.00000 0.011 0.991
## Treatmentc5 2.76909 4.82237 89.00000 0.574 0.567
## Treatmentc10 -2.46860 4.90611 89.00000 -0.503 0.616
##
## Correlation of Fixed Effects:
## (Intr) Trtmn5
## Treatmentc5 -0.701
## Treatmntc10 -0.689 0.484
## convergence code: 0
## boundary (singular) fit: see ?isSingulartapply(OE.EN.Fem$Corr_encaps,OE.EN.Fem$Treatment,mean)## Control c5 c10
## 0.03689062 2.80598065 -2.43171034tapply(OE.EN.Fem$Corr_encaps,OE.EN.Fem$Treatment,se)## Control c5 c10
## 3.871824 3.165520 3.2162795.3) Life histories
5.3.1) Analysis
Prepare
Set the first clutch to be represented as zero. This way the intercept can be interpreted as the average difference in each life history trait at the first clutch.
OE.RP$Clutch_rank0 <- OE.RP$Clutch_rank-15.3.2) Laying time
How does tratement and clutch number affect the timing of the clutches?
OE.model.RP.LayTime.Full <- lmer(yday(Laying_date)~Treatment*Clutch_rank0+
(1|Female_id),data=OE.RP)
## Linearity
plot(resid(OE.model.RP.LayTime.Full),fitted(OE.model.RP.LayTime.Full))
## Homogeneity
plot(OE.model.RP.LayTime.Full)
## Normality
qqnorm(residuals(OE.model.RP.LayTime.Full))
## Influence
infl <- influence(OE.model.RP.LayTime.Full, obs = TRUE)
plot(infl, which = "cook")
- Linearity is fine.
- Homogeneity is fine.
- The qq-line is fine.
- One possible outlier, but the influence is acceptable (Cook’s D < 1).
Get the results
Display the model summary.
Get the relationship between clutch rank and egg laying time for each treatment separately.
summary(OE.model.RP.LayTime.Full)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: yday(Laying_date) ~ Treatment * Clutch_rank0 + (1 | Female_id)
## Data: OE.RP
##
## REML criterion at convergence: 321.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.68741 -0.59968 0.02047 0.72058 2.26887
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female_id (Intercept) 11.475 3.388
## Residual 1.362 1.167
## Number of obs: 82, groups: Female_id, 19
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 123.5645 1.1756 18.8283 105.107
## TreatmentBacteria -4.0863 1.6141 18.5390 -2.532
## Clutch_rank0 3.0621 0.1828 62.4144 16.749
## TreatmentBacteria:Clutch_rank0 -0.4344 0.2113 62.2335 -2.056
## Pr(>|t|)
## (Intercept) <2e-16 ***
## TreatmentBacteria 0.0206 *
## Clutch_rank0 <2e-16 ***
## TreatmentBacteria:Clutch_rank0 0.0439 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmnB Cltc_0
## TretmntBctr -0.728
## Clutch_rnk0 -0.195 0.142
## TrtmntB:C_0 0.169 -0.189 -0.865emtrends(OE.model.RP.LayTime.Full,pairwise~Treatment,var="Clutch_rank0")## $emtrends
## Treatment Clutch_rank0.trend SE df lower.CL upper.CL
## Control 3.06 0.183 62.4 2.70 3.43
## Bacteria 2.63 0.106 61.7 2.42 2.84
##
## Trends are based on the yday (transformed) scale
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Control - Bacteria 0.434 0.211 62.2 2.056 0.0439aggregate(yday(OE.RP$Laying_date), by=list(OE.RP$Treatment,OE.RP$Clutch_rank),FUN=mean)## Group.1 Group.2 x
## 1 Control 1 123.6667
## 2 Bacteria 1 119.2000
## 3 Control 2 126.8750
## 4 Bacteria 2 122.3000
## 5 Control 3 129.5714
## 6 Bacteria 3 124.4444
## 7 Control 4 131.0000
## 8 Bacteria 4 127.1111
## 9 Control 5 133.5000
## 10 Bacteria 5 130.2857
## 11 Bacteria 6 131.0000
## 12 Bacteria 7 136.0000aggregate(yday(OE.RP$Laying_date), by=list(OE.RP$Treatment,OE.RP$Clutch_rank),FUN=se)## Group.1 Group.2 x
## 1 Control 1 1.3437096
## 2 Bacteria 1 0.6110101
## 3 Control 2 1.7872315
## 4 Bacteria 2 0.6155395
## 5 Control 3 1.8107765
## 6 Bacteria 3 0.8012336
## 7 Control 4 2.4899799
## 8 Bacteria 4 0.9345891
## 9 Control 5 4.5000000
## 10 Bacteria 5 1.4261885
## 11 Bacteria 6 1.1547005
## 12 Bacteria 7 NAPlot
Plot the data and add a prediction line including 95% confidence interval.
newdat <- data.frame(Clutch_rank0=rep(0:6,2),Treatment=rep(c("Control","Bacteria"),each=7))
pfun <- function(fit) {
predict(fit,newdata=newdat,re.form=NA)
}
prediction <- cbind(newdat,pred=predict(OE.model.RP.LayTime.Full,newdata=newdat,re.form=NA))
prediction <- cbind(prediction,confint(OE.model.RP.LayTime.Full,method="boot",FUN=pfun))
colnames(prediction)[c(4,5)] <- c("lwr","upr")#svg("ResultsV3/OE.RP.LayDay.Plot.svg",width=10,height=10)
ylimi <- range(yday(OE.RP$Laying_date),na.rm=TRUE)
plot(yday(OE.RP$Laying_date)~jitter(OE.RP$Clutch_rank0,0.6),subset=OE.RP$Treatment=="Bacteria",ylim=ylimi,pch=16,col=2,xlab="Clutch number",ylab="Laying day",las=1)
points(yday(OE.RP$Laying_date)~jitter(OE.RP$Clutch_rank0,0.6),subset=OE.RP$Treatment=="Control",pch=16,col=1)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Bacteria",col=2)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Control",col=1)
with(prediction[prediction$Treatment=="Control",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("black", alpha.f = 0.10), border = NA))
with(prediction[prediction$Treatment=="Bacteria",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("red", alpha.f = 0.10), border = NA))
legend("bottomright",pch=16,col=c(1,2),legend=c("Control","5mg/ml"))
#dev.off()5.3.3) Number of clutches
How many clutches do the females of the different treatments lay in total?
OE.model.RP.NrClutch.Full <- glm(x~Group.2,data=OE.CPF,family=poisson())
## Overdispersion
overdisp_fun(OE.model.RP.NrClutch.Full)## chisq ratio rdf p
## 7.7290323 0.4546490 17.0000000 0.9720025## Linearity
plot(resid(OE.model.RP.NrClutch.Full),fitted(OE.model.RP.NrClutch.Full))
## Homogeneity
plot(OE.model.RP.NrClutch.Full)
## Normality
qqnorm(residuals(OE.model.RP.NrClutch.Full))
## Influence
#infl <- influence(OE.model.RP.NrClutch.Full, obs = TRUE)
#plot(infl, which = "cook")The data is quite severly under disperesed. This has to be fixed with a generalized poisson distribution.
OE.model.RP.NrClutch.Full <- vglm(x~Group.2,data=OE.CPF,family=genpoisson)
## Overdispersion
overdisp_fun(OE.model.RP.NrClutch.Full)## chisq ratio rdf p
## 34.9823713 0.9994963 35.0000000 0.4690395## Linearity
plot(resid(OE.model.RP.NrClutch.Full)[,2],fitted(OE.model.RP.NrClutch.Full))
## Homogeneity
plot(OE.model.RP.NrClutch.Full)
## Influence
hatvalues(OE.model.RP.NrClutch.Full)## rhobitlink(lambda) loglink(theta)
## 1 0.10338894 0.06167172
## 2 0.10338894 0.06167172
## 3 0.10338894 0.06167172
## 4 0.10338894 0.06167172
## 5 0.10338894 0.06167172
## 6 0.10338894 0.06167172
## 7 0.10338894 0.06167172
## 8 0.10338894 0.06167172
## 9 0.10338894 0.06167172
## 10 0.09483447 0.05661094
## 11 0.09483447 0.05661094
## 12 0.09483447 0.05661094
## 13 0.09483447 0.05661094
## 14 0.09483447 0.05661094
## 15 0.09483447 0.05661094
## 16 0.09483447 0.05661094
## 17 0.09483447 0.05661094
## 18 0.09483447 0.05661094
## 19 0.09483447 0.05661094
## attr(,"predictors.names")
## [1] "rhobitlink(lambda)" "loglink(theta)"
## attr(,"ncol.X.vlm")
## [1] 3- Overdispersion: fixed.
- Homogeneity is fine.
- No obvious outliers in the hatvalues.
Get the results
Display the model summary.
To be on the safe side, also perform a wilcoxon rank sum test. This should confirm the outcome.
Get the average number of clutches produced in each treatment.
summary(OE.model.RP.NrClutch.Full)##
## Call:
## vglm(formula = x ~ Group.2, family = genpoisson, data = OE.CPF)
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## rhobitlink(lambda) -0.7088 -0.6500 -0.30999 0.1715 2.357
## loglink(theta) -2.0007 -0.6157 -0.05135 0.7192 1.607
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -1.6636 0.9565 NA NA
## (Intercept):2 1.7645 0.2035 8.672 < 2e-16 ***
## Group.2Bacteria 0.3595 0.1294 2.779 0.00545 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: rhobitlink(lambda), loglink(theta)
##
## Log-likelihood: -31.7701 on 35 degrees of freedom
##
## Number of Fisher scoring iterations: 9
##
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):1'wilcox.test(OE.CPF$x~OE.CPF$Group.2)##
## Wilcoxon rank sum test with continuity correction
##
## data: OE.CPF$x by OE.CPF$Group.2
## W = 17.5, p-value = 0.02442
## alternative hypothesis: true location shift is not equal to 0tapply(OE.CPF$x,OE.CPF$Group.2,mean)## Control Bacteria
## 3.444444 5.000000tapply(OE.CPF$x,OE.CPF$Group.2,se)## Control Bacteria
## 0.4444444 0.44721365.3.4) Number of eggs per clutch
How many eggs do the females lay in each clutch?
OE.model.RP.NrEgg.Full <- glmer(Number_eggs~Treatment*Clutch_rank0+(1|Female_id),
data=OE.RP,family=poisson())
## Overdispersion
overdisp_fun(OE.model.RP.NrEgg.Full)## chisq ratio rdf p
## 1.148701e+03 1.491819e+01 7.700000e+01 1.361454e-190## Linearity
plot(resid(OE.model.RP.NrEgg.Full),fitted(OE.model.RP.NrEgg.Full))
## Homogeneity
plot(OE.model.RP.NrEgg.Full)
## Influence
infl <- influence(OE.model.RP.NrEgg.Full, obs = TRUE)
plot(infl, which = "cook")
The data is overdispersed. Use a negative binomial distribution.
OE.model.RP.NrEgg.Full <- glmer.nb(Number_eggs~Treatment*Clutch_rank0+
(1|Female_id),data=OE.RP)
## Overdispersion
overdisp_fun(OE.model.RP.NrEgg.Full)## chisq ratio rdf p
## 74.0007356 0.9736939 76.0000000 0.5435491## Linearity
plot(resid(OE.model.RP.NrEgg.Full),fitted(OE.model.RP.NrEgg.Full))
## Homogeneity
plot(OE.model.RP.NrEgg.Full)
## Influence
infl <- influence(OE.model.RP.NrEgg.Full, obs = TRUE)
plot(infl, which = "cook")
- Overdispersion: fixed.
- Linearity is fine.
- Homogeneity is fine.
- No obvious outliers.
Get the results
Display the model summary.
Get the relationship between clutch rank and number eggs for each treatment separately.
Display the mean number of eggs per clutch for each treatment.
summary(OE.model.RP.NrEgg.Full)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(7.7983) ( log )
## Formula: Number_eggs ~ Treatment * Clutch_rank0 + (1 | Female_id)
## Data: OE.RP
##
## AIC BIC logLik deviance df.resid
## 884.2 898.6 -436.1 872.2 76
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.98510 -0.57425 -0.00273 0.57214 2.80221
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female_id (Intercept) 2.297e-13 4.792e-07
## Number of obs: 82, groups: Female_id, 19
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.31488 0.09986 53.222 <2e-16 ***
## TreatmentBacteria 0.05585 0.13440 0.416 0.6777
## Clutch_rank0 -0.11608 0.05097 -2.277 0.0228 *
## TreatmentBacteria:Clutch_rank0 -0.14701 0.06138 -2.395 0.0166 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmnB Cltc_0
## TretmntBctr -0.743
## Clutch_rnk0 -0.761 0.566
## TrtmntB:C_0 0.632 -0.773 -0.830
## convergence code: 0
## boundary (singular) fit: see ?isSingularemtrends(OE.model.RP.NrEgg.Full,pairwise~Treatment,var="Clutch_rank0")## $emtrends
## Treatment Clutch_rank0.trend SE df asymp.LCL asymp.UCL
## Control -0.116 0.0510 Inf -0.216 -0.0162
## Bacteria -0.263 0.0342 Inf -0.330 -0.1961
##
## Trends are based on the log (transformed) scale
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## Control - Bacteria 0.147 0.0614 Inf 2.395 0.0166aggregate(OE.RP$Number_eggs, by=list(OE.RP$Treatment,OE.RP$Clutch_rank),FUN=mean)## Group.1 Group.2 x
## 1 Control 1 220.55556
## 2 Bacteria 1 179.90000
## 3 Control 2 149.37500
## 4 Bacteria 2 168.30000
## 5 Control 3 173.28571
## 6 Bacteria 3 150.44444
## 7 Control 4 145.83333
## 8 Bacteria 4 113.66667
## 9 Control 5 129.00000
## 10 Bacteria 5 64.71429
## 11 Bacteria 6 57.50000
## 12 Bacteria 7 11.00000aggregate(OE.RP$Number_eggs, by=list(OE.RP$Treatment,OE.RP$Clutch_rank),FUN=se)## Group.1 Group.2 x
## 1 Control 1 16.149571
## 2 Bacteria 1 18.274725
## 3 Control 2 16.439214
## 4 Bacteria 2 19.728744
## 5 Control 3 29.585803
## 6 Bacteria 3 7.717425
## 7 Control 4 28.917603
## 8 Bacteria 4 10.873004
## 9 Control 5 26.000000
## 10 Bacteria 5 9.430736
## 11 Bacteria 6 15.755951
## 12 Bacteria 7 NAPlot
Plot the data and add a prediction line including 95% confidence interval.
newdat <- data.frame(Clutch_rank0=rep(0:6,2),Treatment=rep(c("Control","Bacteria"),each=7))
pfun <- function(fit) {
predict(fit,newdata=newdat,re.form=NA,type="response")
}
prediction <- cbind(newdat,pred=predict(OE.model.RP.NrEgg.Full,newdata=newdat,re.form=NA,type="response"))
prediction <- cbind(prediction,confint(OE.model.RP.NrEgg.Full,method="boot",FUN=pfun))
colnames(prediction)[c(4,5)] <- c("lwr","upr")#svg("ResultsV3/OE.RP.NrEgg.Plot.svg",width=10,height=10)
ylimi <- range(OE.RP$Number_eggs,na.rm=TRUE)
plot(OE.RP$Number_eggs~jitter(OE.RP$Clutch_rank0,0.6),subset=OE.RP$Treatment=="Bacteria",pch=16,col=2,xlab="Clutch number",ylab="Number eggs",las=1)
points(OE.RP$Number_eggs~jitter(OE.RP$Clutch_rank0,0.6),subset=OE.RP$Treatment=="Control",pch=16,col=1)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Bacteria",col=2)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Control",col=1)
with(prediction[prediction$Treatment=="Control",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("black", alpha.f = 0.10), border = NA))
with(prediction[prediction$Treatment=="Bacteria",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("red", alpha.f = 0.10), border = NA))
legend("topright",pch=16,col=c(1,2),legend=c("Control","5mg/ml"))
#dev.off()5.3.5) Lifetime egg production
How does lifetime egg production compare between the two treatments?
Given the low sample size, the outcome will not be very reliable.
tapply(OE.LEP$x,OE.LEP$Group.2,mean)## Control Bacteria
## 614.0 655.3tapply(OE.LEP$x,OE.LEP$Group.2,median)## Control Bacteria
## 660.0 670.5tapply(OE.LEP$x,OE.LEP$Group.2,se)## Control Bacteria
## 83.89974 35.66014ci.calc <- function(x){boot.ci(boot(x,function(x,i) median(x[i]), R=1000))}
tapply(OE.LEP$x,OE.LEP$Group.2,ci.calc)## $Control
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot(x, function(x, i) median(x[i]), R = 1000))
##
## Intervals :
## Level Normal Basic
## 95% ( 483.0, 886.1 ) ( 399.0, 1003.0 )
##
## Level Percentile BCa
## 95% (317, 921 ) (317, 753 )
## Calculations and Intervals on Original Scale
## Some BCa intervals may be unstable
##
## $Bacteria
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot(x, function(x, i) median(x[i]), R = 1000))
##
## Intervals :
## Level Normal Basic
## 95% (605.1, 734.6 ) (600.0, 745.0 )
##
## Level Percentile BCa
## 95% (596.0, 741.0 ) (576.0, 735.2 )
## Calculations and Intervals on Original Scalewilcox.test(OE.LEP$x~OE.LEP$Group.2)##
## Wilcoxon rank sum test
##
## data: OE.LEP$x by OE.LEP$Group.2
## W = 40, p-value = 0.7197
## alternative hypothesis: true location shift is not equal to 05.3.5) Average egg weigt
What’s the average weight ef the eggs in each clutch and treatment?
OE.model.RP.Weight.Full <- lmer(Weight_per_egg~Treatment*Clutch_rank0+
(1|Female_id),data=OE.RP)
## Linearity
plot(resid(OE.model.RP.Weight.Full),fitted(OE.model.RP.Weight.Full))
## Homogeneity
plot(OE.model.RP.Weight.Full)
## Normality
qqnorm(residuals(OE.model.RP.Weight.Full))
## Influence
infl <- influence(OE.model.RP.Weight.Full, obs = TRUE)
plot(infl, which = "cook")
- Linearity is fine, but Outliers?
- Homogeneity could be better. Outlier?
- The qq-line is straight except for two points. Outlier?
- One possible outlier, but the influence is acceptable (Cook’s D < 1).
Check
OE.RP[c(24),]## Clutch_number Female_id Mating_date Treatment Laying_date Clutch_rank
## 24 24 101 2016-05-01 Bacteria 2016-05-01 1
## Number_eggs Weight_mg Hatching_date Number_hatched Time_until_hatching
## 24 90 0.0125 2016-05-13 80 12
## Hatching_percentage Time_until_laying Time_between_clutches
## 24 88.89 0 NA
## Weight_per_egg Notes Clutch_rank0
## 24 0.139 0OE.RP.Red <- OE.RP[-c(24),]
OE.model.RP.Weight.Red <- lmer(Weight_per_egg~Treatment*Clutch_rank0+
(1|Female_id),data=OE.RP.Red)
## Linearity
plot(resid(OE.model.RP.Weight.Red),fitted(OE.model.RP.Weight.Red))
## Homogeneity
plot(OE.model.RP.Weight.Red)
## Normality
qqnorm(residuals(OE.model.RP.Weight.Red))
## Influence
infl <- influence(OE.model.RP.Weight.Red, obs = TRUE)
plot(infl, which = "cook")
Removal of the outlier does not severly affect the model outcome.
Keep the full model, as there is no rason to remove outliers.
Get the results
Display the model summary.
Get the relationship between clutch rank and average egg weight for each treatment separately.
Display the mean egg weight per clutch for each treatment.
summary(OE.model.RP.Weight.Full)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Weight_per_egg ~ Treatment * Clutch_rank0 + (1 | Female_id)
## Data: OE.RP
##
## REML criterion at convergence: -490.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8644 -0.5614 0.0076 0.4864 3.7212
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female_id (Intercept) 4.752e-05 0.006894
## Residual 6.599e-05 0.008124
## Number of obs: 82, groups: Female_id, 19
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.090747 0.003226 33.967511 28.130
## TreatmentBacteria 0.005418 0.004337 31.363376 1.249
## Clutch_rank0 -0.001965 0.001232 68.881409 -1.596
## TreatmentBacteria:Clutch_rank0 -0.001970 0.001430 68.093592 -1.378
## Pr(>|t|)
## (Intercept) <2e-16 ***
## TreatmentBacteria 0.221
## Clutch_rank0 0.115
## TreatmentBacteria:Clutch_rank0 0.173
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmnB Cltc_0
## TretmntBctr -0.744
## Clutch_rnk0 -0.513 0.382
## TrtmntB:C_0 0.442 -0.503 -0.862emtrends(OE.model.RP.Weight.Full,pairwise~Treatment,var="Clutch_rank0")## $emtrends
## Treatment Clutch_rank0.trend SE df lower.CL upper.CL
## Control -0.00197 0.001232 68.9 -0.00442 0.000492
## Bacteria -0.00394 0.000726 65.7 -0.00538 -0.002487
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Control - Bacteria 0.00197 0.00143 68.1 1.378 0.1727aggregate(OE.RP$Weight_per_egg, by=list(OE.RP$Treatment,OE.RP$Clutch_rank),FUN=mean)## Group.1 Group.2 x
## 1 Control 1 0.09411111
## 2 Bacteria 1 0.09530000
## 3 Control 2 0.08175000
## 4 Bacteria 2 0.09140000
## 5 Control 3 0.08742857
## 6 Bacteria 3 0.08866667
## 7 Control 4 0.08883333
## 8 Bacteria 4 0.08700000
## 9 Control 5 0.08450000
## 10 Bacteria 5 0.08442857
## 11 Bacteria 6 0.07250000
## 12 Bacteria 7 0.07300000aggregate(OE.RP$Weight_per_egg, by=list(OE.RP$Treatment,OE.RP$Clutch_rank),FUN=se)## Group.1 Group.2 x
## 1 Control 1 0.002642833
## 2 Bacteria 1 0.006055392
## 3 Control 2 0.002788689
## 4 Bacteria 2 0.003070288
## 5 Control 3 0.002644465
## 6 Bacteria 3 0.002166667
## 7 Control 4 0.003145544
## 8 Bacteria 4 0.002236068
## 9 Control 5 0.003500000
## 10 Bacteria 5 0.005558593
## 11 Bacteria 6 0.003617089
## 12 Bacteria 7 NAPlot
Plot the data and add a prediction line including 95% confidence interval.
newdat <- data.frame(Clutch_rank0=rep(0:6,2),Treatment=rep(c("Control","Bacteria"),each=7))
pfun <- function(fit) {
predict(fit,newdata=newdat,re.form=NA,type="response")
}
prediction <- cbind(newdat,pred=predict(OE.model.RP.Weight.Full,newdata=newdat,re.form=NA,type="response"))
prediction <- cbind(prediction,confint(OE.model.RP.Weight.Full,method="boot",FUN=pfun))
colnames(prediction)[c(4,5)] <- c("lwr","upr")#svg("ResultsV3/OE.RP.Weight.Plot.svg",width=10,height=10)
ylimi <- range(OE.RP$Weight_per_egg,na.rm=TRUE)
plot(OE.RP$Weight_per_egg~jitter(OE.RP$Clutch_rank0,0.6),subset=OE.RP$Treatment=="Bacteria",ylim=c(0.06,0.14),pch=16,col=2,xlab="Clutch number",ylab="Average egg weight (ng)",las=1)
points(OE.RP$Weight_per_egg~jitter(OE.RP$Clutch_rank0,0.6),subset=OE.RP$Treatment=="Control",pch=16,col=1)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Bacteria",col=2)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Control",col=1)
with(prediction[prediction$Treatment=="Control",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("black", alpha.f = 0.10), border = NA))
with(prediction[prediction$Treatment=="Bacteria",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("red", alpha.f = 0.10), border = NA))
legend("topright",pch=16,col=c(1,2),legend=c("Control","5mg/ml"))
#dev.off()5.3.6) Hatching success
How well do the larvae of eatch clutch hatch in each treatment?
OE.RP.NoNA <- na.omit(OE.RP[,c(2,4,7,10,12,17)]) # there are no NAs, this is just for the prediction
OE.model.RP.Hatch.Full <- glmer(cbind(Number_hatched,(Number_eggs-Number_hatched))~
Treatment*Clutch_rank0+(1|Female_id),
data=OE.RP.NoNA,family=binomial())
## Linearity
plot(resid(OE.model.RP.Hatch.Full),fitted(OE.model.RP.Hatch.Full))
## Homogeneity
plot(OE.model.RP.Hatch.Full)
## Normality
qqnorm(residuals(OE.model.RP.Hatch.Full))
## Outlier
plot(hatvalues(OE.model.RP.Hatch.Full)~residuals(OE.model.RP.Hatch.Full))
## Influence
infl <- influence(OE.model.RP.Hatch.Full, obs = TRUE)
plot(infl, which = "cook")
- Ignore Linearity as it is not important for binomial regression.
- Homogeneity is fine.
- The qq line is fine.
- Four outliers. Cook’s D above 1.
Check
OE.RP.NoNA[c(49,52,66),]## Female_id Treatment Number_eggs Number_hatched Hatching_percentage
## 49 90 Control 180 110 61.11
## 52 85 Control 290 48 16.55
## 66 105 Control 273 212 77.66
## Clutch_rank0
## 49 0
## 52 0
## 66 2OE.RP.NoNA.Red <- OE.RP.NoNA[-c(49,52,66),]
OE.model.RP.Hatch.Red <- glmer(cbind(Number_hatched,(Number_eggs-Number_hatched))~
Treatment*Clutch_rank0+(1|Female_id),
data=OE.RP.NoNA.Red,family=binomial())
## Influence
infl <- influence(OE.model.RP.Hatch.Red, obs = TRUE)
plot(infl, which = "cook")
After removal of the outliers there are new values indicated as outliers. Since these three outliers don’t affect the model too strongly, and the whole process may run into a loop with constant new outliers, I’ll stick to the full model.
Get the results
Display the model summary.
Get the relationship between clutch rank and hatching success for each treatment separately.
Display the mean hatching success per clutch for each treatment.
summary(OE.model.RP.Hatch.Full)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## cbind(Number_hatched, (Number_eggs - Number_hatched)) ~ Treatment *
## Clutch_rank0 + (1 | Female_id)
## Data: OE.RP.NoNA
##
## AIC BIC logLik deviance df.resid
## 1209.4 1221.4 -599.7 1199.4 76
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.2389 -1.6669 0.3321 2.0003 6.9639
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female_id (Intercept) 0.6872 0.829
## Number of obs: 81, groups: Female_id, 19
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.78725 0.28085 2.803 0.00506 **
## TreatmentBacteria 1.52904 0.38938 3.927 8.61e-05 ***
## Clutch_rank0 0.20803 0.03311 6.283 3.33e-10 ***
## TreatmentBacteria:Clutch_rank0 -0.40669 0.04295 -9.469 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmnB Cltc_0
## TretmntBctr -0.721
## Clutch_rnk0 -0.111 0.080
## TrtmntB:C_0 0.085 -0.142 -0.771emtrends(OE.model.RP.Hatch.Full,pairwise~Treatment,var="Clutch_rank0",transform="response")## $emtrends
## Treatment Clutch_rank0.trend SE df asymp.LCL asymp.UCL
## Control 0.0373 0.00594 Inf 0.0257 0.0489
## Bacteria -0.0219 0.00301 Inf -0.0278 -0.0160
##
## Trends are obtained after back-transforming from the logit scale
## Results are given on the ( (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## Control - Bacteria 0.0592 0.00666 Inf 8.890 <.0001
##
## Results are given on the ( (not the response) scale.aggregate(OE.RP.NoNA$Hatching_percentage, by=list(OE.RP.NoNA$Clutch_rank,OE.RP.NoNA$Treatment),FUN=mean)## Group.1 Group.2 x
## 1 0 Control 65.37222
## 2 1 Control 70.05500
## 3 2 Control 78.30833
## 4 3 Control 81.93833
## 5 4 Control 76.45000
## 6 0 Bacteria 90.37800
## 7 1 Bacteria 87.83800
## 8 2 Bacteria 88.32111
## 9 3 Bacteria 84.69222
## 10 4 Bacteria 78.72000
## 11 5 Bacteria 80.94000
## 12 6 Bacteria 63.64000aggregate(OE.RP.NoNA$Hatching_percentage, by=list(OE.RP.NoNA$Clutch_rank,OE.RP.NoNA$Treatment),FUN=se)## Group.1 Group.2 x
## 1 0 Control 10.190074
## 2 1 Control 6.782383
## 3 2 Control 9.633534
## 4 3 Control 5.547081
## 5 4 Control 21.610000
## 6 0 Bacteria 1.739877
## 7 1 Bacteria 2.947389
## 8 2 Bacteria 2.541699
## 9 3 Bacteria 2.230284
## 10 4 Bacteria 5.212932
## 11 5 Bacteria 9.035608
## 12 6 Bacteria NAPlot
Plot the data and add a prediction line including 95% confidence interval.
newdat <- data.frame(Clutch_rank0=rep(0:6,2),Treatment=rep(c("Control","Bacteria"),each=7))
pfun <- function(fit) {
predict(fit,newdata=newdat,re.form=NA,type="response")
}
prediction <- cbind(newdat,pred=predict(OE.model.RP.Hatch.Full,newdata=newdat,re.form=NA,type="response"))
prediction <- cbind(prediction,confint(OE.model.RP.Hatch.Full,method="boot",FUN=pfun))
colnames(prediction)[c(4,5)] <- c("lwr","upr")
prediction[,3:5] <- prediction[,3:5]*100#svg("ResultsV3/OE.RP.Hatch.Plot.svg",width=10,height=10)
ylimi <- range(OE.RP.NoNA$Hatching_percentage,na.rm=TRUE)
plot(OE.RP.NoNA$Hatching_percentage~jitter(OE.RP.NoNA$Clutch_rank0,0.6),subset=OE.RP.NoNA$Treatment=="Bacteria",pch=16,col=2,ylim=ylimi,xlim=c(0,6),xlab="Clutch number",ylab="Hatching success (%)",las=1)
points(OE.RP.NoNA$Hatching_percentage~jitter(OE.RP.NoNA$Clutch_rank0,0.6),subset=OE.RP.NoNA$Treatment=="Control",pch=16,col=1)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Bacteria",col=2)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Control",col=1)
with(prediction[prediction$Treatment=="Control",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("black", alpha.f = 0.10), border = NA))
with(prediction[prediction$Treatment=="Bacteria",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("red", alpha.f = 0.10), border = NA))
legend("bottomright",pch=16,col=c(1,2),legend=c("Control","5mg/ml"))
#dev.off()5.4) Lifespan
5.4.1) Analyse
Check the proportional hazard assumption.
OE.model.LS.Full <- coxph(Surv(Day,Censor)~Feeding_treatment+frailty(Family),data=OE.LS.Fem)
cox.zph(OE.model.LS.Full)## rho chisq p
## Feeding_treatmentc5 -0.201 2.69 0.101
## Feeding_treatmentc10 -0.166 1.77 0.183
## GLOBAL NA 3.06 0.216Run the model and get the results.
OE.model.LS.Fem <- coxme(Surv(Day,Censor)~Feeding_treatment+(1|Family),data=OE.LS.Fem)
summary(OE.model.LS.Fem)## Cox mixed-effects model fit by maximum likelihood
## Data: OE.LS.Fem
## events, n = 68, 68
## Iterations= 5 23
## NULL Integrated Fitted
## Log-likelihood -221.9564 -219.7168 -215.3995
##
## Chisq df p AIC BIC
## Integrated loglik 4.48 3.00 0.214140 -1.52 -8.18
## Penalized loglik 13.11 5.43 0.029453 2.25 -9.80
##
## Model: Surv(Day, Censor) ~ Feeding_treatment + (1 | Family)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## Feeding_treatmentc5 0.4524649 1.572183 0.3287817 1.38 0.170
## Feeding_treatmentc10 0.6212852 1.861319 0.3156290 1.97 0.049
##
## Random effects
## Group Variable Std Dev Variance
## Family Intercept 0.3548425 0.1259132Plot
#svg("ResultsV3/OE.LS.Plot.svg",width=10,height=10)
xlimi <- c(0,max(OE.LS.Fem$Day,na.rm=TRUE))
modsum <- summary(survfit(Surv(Day,Censor)~Feeding_treatment+(1|Family),data=OE.LS.Fem),times=0:50)
with(modsum,plot(time,surv,type="n",xlim=xlimi,ylim=c(0,1),las=1,xlab="Days",ylab="Cumulative survival"))
polygon(c(modsum[["time"]][1:34],rev(modsum[["time"]][1:34])),c(modsum[["lower"]][1:34],rev(modsum[["upper"]][1:34])),
col = adjustcolor("black", alpha.f = 0.10), border = FALSE)
polygon(c(modsum[["time"]][36:84],rev(modsum[["time"]][36:84])),c(modsum[["lower"]][36:84],rev(modsum[["upper"]][36:84])),
col = adjustcolor("red", alpha.f = 0.10), border = FALSE)
polygon(c(modsum[["time"]][86:117],rev(modsum[["time"]][86:117])),c(modsum[["lower"]][86:117],rev(modsum[["upper"]][86:117])),
col = adjustcolor("blue", alpha.f = 0.10), border = FALSE)
lines(modsum[["time"]][1:35],modsum[["surv"]][1:35],type="s")
lines(modsum[["time"]][36:85],modsum[["surv"]][36:85],type="s",col="red")
lines(modsum[["time"]][86:118],modsum[["surv"]][86:118],type="s",col="blue")
legend("topright",legend=c("Control","5mg/ml","10mg/ml"),pch=c(15),col=c("black","red","blue"),cex=1)
#dev.off()6) Haemocoelic exposure
6.1) Gene expression
6.1.1) Assess PCA
## Target gene columns: 3,4,6,8,9,10,12
cortest.bartlett(HE.GE.Fem.Mat[,c(3,4,6,8,9,10,12)])## $chisq
## [1] 122.4893
##
## $p.value
## [1] 2.490741e-16
##
## $df
## [1] 21det(cor(as.matrix(HE.GE.Fem.Mat[,c(3,4,6,8,9,10,12)])))## [1] 0.1071155KMO(HE.GE.Fem.Mat[,c(3,4,6,8,9,10,12)])## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = HE.GE.Fem.Mat[, c(3, 4, 6, 8, 9, 10, 12)])
## Overall MSA = 0.46
## MSA for each item =
## Attacin B1.3g Lysozyme Pelle PGRP proPO Serpin
## 0.43 0.31 0.48 0.65 0.51 0.55 0.38Colinearity is okay, but the sampling adequacy is not. Removal of genes may help, but with so few genes the remaining PCA would not tell much more than a gene-by-gene analysis.
6.1.2) Gene by gene; ad hoc mixed linear models
Only the model specification stepp is shown here. These models are primarily used for the interaction effects, as the pairwise comparisons are of primary interest.
HE.model.adhoc.Attacin.Full <- lmer(Attacin~(Treatment+Timepoint)^2+(1|Family),
data=HE.GE.Fem.Mat)
HE.model.adhoc.B1.3g.Full <- lmer(B1.3g~(Treatment+Timepoint)^2+(1|Family),
data=HE.GE.Fem.Mat)
HE.model.adhoc.Lysozyme.Full <- lmer(Lysozyme~(Treatment+Timepoint)^2+(1|Family),
data=HE.GE.Fem.Mat)
HE.model.adhoc.Pelle.Full <- lmer(Pelle~(Treatment+Timepoint)^2+(1|Family),
data=HE.GE.Fem.Mat)
HE.model.adhoc.PGRP.Full <- lmer(PGRP~(Treatment+Timepoint)^2+(1|Family),
data=HE.GE.Fem.Mat)
HE.model.adhoc.proPO.Full <- lmer(proPO~(Treatment+Timepoint)^2+(1|Family),
data=HE.GE.Fem.Mat)
HE.model.adhoc.Serpin.Full <- lmer(Serpin~(Treatment+Timepoint)^2+(1|Family),
data=HE.GE.Fem.Mat)Adjust p-values
Adjust all p-values for the comparison of seven genes using false discovery rates.
pvals.adho <- rbind(Attacin=coef(summary(HE.model.adhoc.Attacin.Full))[,5],
B1.3g=coef(summary(HE.model.adhoc.B1.3g.Full))[,5],
Lysozyme=coef(summary(HE.model.adhoc.Lysozyme.Full))[,5],
Pelle=coef(summary(HE.model.adhoc.Pelle.Full))[,5],
PGRP=coef(summary(HE.model.adhoc.PGRP.Full))[,5],
proPO=coef(summary(HE.model.adhoc.proPO.Full))[,5],
Serpin=coef(summary(HE.model.adhoc.Serpin.Full))[,5])
qvals.adho<-apply(pvals.adho,2,p.adjust,method="fdr")Validate the models
1) Attacin:
## Attacin
## Linearity
plot(resid(HE.model.adhoc.Attacin.Full),fitted(HE.model.adhoc.Attacin.Full))
## Homogeneity
plot(HE.model.adhoc.Attacin.Full)
## Normality
qqnorm(residuals(HE.model.adhoc.Attacin.Full))
## Outlier
plot(hatvalues(HE.model.adhoc.Attacin.Full)~residuals(HE.model.adhoc.Attacin.Full))
- Linearity is fine.
- Homogeneity is fine.
- The qq-line looks okay for analysis.
- No notable outliers.
## B1.3g
## Linearity
plot(resid(HE.model.adhoc.B1.3g.Full),fitted(HE.model.adhoc.B1.3g.Full))
## Homogeneity
plot(HE.model.adhoc.B1.3g.Full)
## Normality
qqnorm(residuals(HE.model.adhoc.B1.3g.Full))
## Outlier
plot(hatvalues(HE.model.adhoc.B1.3g.Full)~residuals(HE.model.adhoc.B1.3g.Full))
- Linearity is fine.
- Homogeneity is fine.
- The qq-line looks okay for analysis.
- No notable outliers.
## Lysozyme
## Linearity
plot(resid(HE.model.adhoc.Lysozyme.Full),fitted(HE.model.adhoc.Lysozyme.Full))
## Homogeneity
plot(HE.model.adhoc.Lysozyme.Full)
## Normality
qqnorm(residuals(HE.model.adhoc.Lysozyme.Full))
## Outlier
plot(hatvalues(HE.model.adhoc.Lysozyme.Full)~residuals(HE.model.adhoc.Lysozyme.Full))
- Linearity is fine.
- Homogeneity is fine.
- The qq-line looks okay for analysis.
- No notable outliers.
## Pelle
## Linearity
plot(resid(HE.model.adhoc.Pelle.Full),fitted(HE.model.adhoc.Pelle.Full))
## Homogeneity
plot(HE.model.adhoc.Pelle.Full)
## Normality
qqnorm(residuals(HE.model.adhoc.Pelle.Full))
## Outlier
plot(hatvalues(HE.model.adhoc.Pelle.Full)~residuals(HE.model.adhoc.Pelle.Full))
- Linearity is fine.
- Homogeneity looks good enough.
- The qq-line looks okay for analysis.
- No notable outliers.
## PGRP
## Linearity
plot(resid(HE.model.adhoc.PGRP.Full),fitted(HE.model.adhoc.PGRP.Full))
## Homogeneity
plot(HE.model.adhoc.PGRP.Full)
## Normality
qqnorm(residuals(HE.model.adhoc.PGRP.Full))
## Outlier
plot(hatvalues(HE.model.adhoc.PGRP.Full)~residuals(HE.model.adhoc.PGRP.Full))
- Linearity is fine.
- Homogeneity is fine.
- The qq-line looks okay for analysis.
- No notable outliers.
## proPO
## Linearity
plot(resid(HE.model.adhoc.proPO.Full),fitted(HE.model.adhoc.proPO.Full))
## Homogeneity
plot(HE.model.adhoc.proPO.Full)
## Normality
qqnorm(residuals(HE.model.adhoc.proPO.Full))
## Outlier
plot(hatvalues(HE.model.adhoc.proPO.Full)~residuals(HE.model.adhoc.proPO.Full))
- Linearity is fine.
- Homogeneity looks good enough.
- The qq-line looks okay for analysis.
- No notable outliers.
## Serpin
## Linearity
plot(resid(HE.model.adhoc.Serpin.Full),fitted(HE.model.adhoc.Serpin.Full))
## Homogeneity
plot(HE.model.adhoc.Serpin.Full)
## Normality
qqnorm(residuals(HE.model.adhoc.Serpin.Full))
## Outlier
plot(hatvalues(HE.model.adhoc.Serpin.Full)~residuals(HE.model.adhoc.Serpin.Full))
- Linearity is fine.
- Homogeneity is slightly off. But everthing else is fine.
- The qq-line looks okay for analysis.
- No notable outliers.
6.1.3) Gene by gene; post hoc pairwise comparisons
Create post hoc multiple comparisons for a subset of pre-specified comparisons.
Adjust the p-values for an analysis of seven genes, using the fdr method.
comps<-c(
"t24.naive - t24.PBS = 0","t24.naive - t24.bacteria = 0",
"t72.naive - t72.PBS = 0","t72.naive - t72.bacteria = 0",
"t24.naive - t72.naive = 0","t24.PBS - t72.PBS = 0",
"t24.bacteria - t72.bacteria = 0"
)
HE.model.posthoc.Attacin <- lmer(Attacin~FullTreatment+(1|Family),data=HE.GE.Fem.Mat)
summ.Att <- summary(glht(HE.model.posthoc.Attacin,linfct=mcp(FullTreatment=comps)),test=adjusted(type="none"))
HE.model.posthoc.B1.3g <- lmer(B1.3g~FullTreatment+(1|Family),data=HE.GE.Fem.Mat)
summ.B13 <- summary(glht(HE.model.posthoc.B1.3g,linfct=mcp(FullTreatment=comps)),test=adjusted(type="none"))
HE.model.posthoc.Lysozyme <- lmer(Lysozyme~FullTreatment+(1|Family),data=HE.GE.Fem.Mat)
summ.Lys <- summary(glht(HE.model.posthoc.Lysozyme,linfct=mcp(FullTreatment=comps)),test=adjusted(type="none"))
HE.model.posthoc.Pelle <- lmer(Pelle~FullTreatment+(1|Family),data=HE.GE.Fem.Mat)
summ.Pel <- summary(glht(HE.model.posthoc.Pelle,linfct=mcp(FullTreatment=comps)),test=adjusted(type="none"))
HE.model.posthoc.PGRP <- lmer(PGRP~FullTreatment+(1|Family),data=HE.GE.Fem.Mat)
summ.PGR <- summary(glht(HE.model.posthoc.PGRP,linfct=mcp(FullTreatment=comps)),test=adjusted(type="none"))
HE.model.posthoc.proPO <- lmer(proPO~FullTreatment+(1|Family),data=HE.GE.Fem.Mat)
summ.pro <- summary(glht(HE.model.posthoc.proPO,linfct=mcp(FullTreatment=comps)),test=adjusted(type="none"))
HE.model.posthoc.Serpin <- lmer(Serpin~FullTreatment+(1|Family),data=HE.GE.Fem.Mat)
summ.Ser <- summary(glht(HE.model.posthoc.Serpin,linfct=mcp(FullTreatment=comps)),test=adjusted(type="none"))
pvals.poho <- rbind(Attacin=summ.Att$test$pvalues,
B1.3g=summ.B13$test$pvalues,
Lysozyme=summ.Lys$test$pvalues,
Pelle=summ.Pel$test$pvalues,
PGRP=summ.PGR$test$pvalues,
proPO=summ.pro$test$pvalues,
Serpin=summ.Ser$test$pvalues)
qvals.poho<-t(apply(pvals.poho,2,p.adjust,method="fdr"))Create table
Create a list with the results and write it to a table.
HE.GE.Fem.GxG.PH.Results <- list(Attacin=
summ.Att$test[c("coefficients","sigma","tstat","pvalues")] %>%
{cbind(estimate=
{paste(format(round(.$coefficients,digits=2),nsmall=2),
format(round(.$sigma,digits=2),nsmall=2),sep="±")},
tvalue=format(round(.$tstat,digits=4),nsmall=4),
pvalue=format(round(.$pvalues,digits=4),nsmall=4),
qvalue=format(round(qvals.poho[,"Attacin"],digits=4),nsmall=4))},
B1.3g=
summ.B13$test[c("coefficients","sigma","tstat","pvalues")] %>%
{cbind(estimate=
{paste(format(round(.$coefficients,digits=2),nsmall=2),
format(round(.$sigma,digits=2),nsmall=2),sep="±")},
tvalue=format(round(.$tstat,digits=4),nsmall=4),
pvalue=format(round(.$pvalues,digits=4),nsmall=4),
qvalue=format(round(qvals.poho[,"B1.3g"],digits=4),nsmall=4))},
Lysozyme=
summ.Lys$test[c("coefficients","sigma","tstat","pvalues")] %>%
{cbind(estimate=
{paste(format(round(.$coefficients,digits=2),nsmall=2),
format(round(.$sigma,digits=2),nsmall=2),sep="±")},
tvalue=format(round(.$tstat,digits=4),nsmall=4),
pvalue=format(round(.$pvalues,digits=4),nsmall=4),
qvalue=format(round(qvals.poho[,"Lysozyme"],digits=4),nsmall=4))},
Pelle=
summ.Pel$test[c("coefficients","sigma","tstat","pvalues")] %>%
{cbind(estimate=
{paste(format(round(.$coefficients,digits=2),nsmall=2),
format(round(.$sigma,digits=2),nsmall=2),sep="±")},
tvalue=format(round(.$tstat,digits=4),nsmall=4),
pvalue=format(round(.$pvalues,digits=4),nsmall=4),
qvalue=format(round(qvals.poho[,"Pelle"],digits=4),nsmall=4))},
PGRP=
summ.PGR$test[c("coefficients","sigma","tstat","pvalues")] %>%
{cbind(estimate=
{paste(format(round(.$coefficients,digits=2),nsmall=2),
format(round(.$sigma,digits=2),nsmall=2),sep="±")},
tvalue=format(round(.$tstat,digits=4),nsmall=4),
pvalue=format(round(.$pvalues,digits=4),nsmall=4),
qvalue=format(round(qvals.poho[,"PGRP"],digits=4),nsmall=4))},
proPO=
summ.pro$test[c("coefficients","sigma","tstat","pvalues")] %>%
{cbind(estimate=
{paste(format(round(.$coefficients,digits=2),nsmall=2),
format(round(.$sigma,digits=2),nsmall=2),sep="±")},
tvalue=format(round(.$tstat,digits=4),nsmall=4),
pvalue=format(round(.$pvalues,digits=4),nsmall=4),
qvalue=format(round(qvals.poho[,"proPO"],digits=4),nsmall=4))},
Serpin=
summ.Ser$test[c("coefficients","sigma","tstat","pvalues")] %>%
{cbind(estimate=
{paste(format(round(.$coefficients,digits=2),nsmall=2),
format(round(.$sigma,digits=2),nsmall=2),sep="±")},
tvalue=format(round(.$tstat,digits=4),nsmall=4),
pvalue=format(round(.$pvalues,digits=4),nsmall=4),
qvalue=format(round(qvals.poho[,"Serpin"],digits=4),nsmall=4))}
)
HE.GE.Fem.GxG.PH.Results## $Attacin
## estimate tvalue pvalue qvalue
## t24.naive - t24.PBS "-4.89±1.31" "-3.7305" "0.0002" "0.0013"
## t24.naive - t24.bacteria "-7.18±1.28" "-5.5923" "0.0000" "0.0000"
## t72.naive - t72.PBS "-4.41±1.33" "-3.3221" "0.0009" "0.0063"
## t72.naive - t72.bacteria "-3.07±1.32" "-2.3286" "0.0199" "0.0464"
## t24.naive - t72.naive "-0.81±1.32" "-0.6163" "0.5377" "0.7960"
## t24.PBS - t72.PBS "-0.32±1.30" "-0.2487" "0.8036" "0.8036"
## t24.bacteria - t72.bacteria " 3.30±1.29" " 2.5603" "0.0105" "0.0306"
##
## $B1.3g
## estimate tvalue pvalue qvalue
## t24.naive - t24.PBS "-1.79±0.62" "-2.8868" "0.0039" "0.0136"
## t24.naive - t24.bacteria "-2.24±0.62" "-3.6255" "0.0003" "0.0004"
## t72.naive - t72.PBS "-1.09±0.64" "-1.7093" "0.0874" "0.2563"
## t72.naive - t72.bacteria "-1.14±0.64" "-1.7905" "0.0734" "0.1027"
## t24.naive - t72.naive " 0.43±0.64" " 0.6759" "0.4991" "0.7960"
## t24.PBS - t72.PBS " 1.13±0.62" " 1.8251" "0.0680" "0.1586"
## t24.bacteria - t72.bacteria " 1.53±0.62" " 2.4804" "0.0131" "0.0306"
##
## $Lysozyme
## estimate tvalue pvalue qvalue
## t24.naive - t24.PBS " 0.04±0.39" " 0.1074" "0.9145" "0.9145"
## t24.naive - t24.bacteria " 0.23±0.37" " 0.6290" "0.5294" "0.6176"
## t72.naive - t72.PBS " 0.51±0.38" " 1.3314" "0.1831" "0.2563"
## t72.naive - t72.bacteria "-0.12±0.38" "-0.3142" "0.7533" "0.7533"
## t24.naive - t72.naive "-0.15±0.38" "-0.4094" "0.6823" "0.7960"
## t24.PBS - t72.PBS " 0.32±0.39" " 0.8040" "0.4214" "0.5899"
## t24.bacteria - t72.bacteria "-0.51±0.37" "-1.3563" "0.1750" "0.1750"
##
## $Pelle
## estimate tvalue pvalue qvalue
## t24.naive - t24.PBS "-1.42±0.53" "-2.6694" "0.0076" "0.0177"
## t24.naive - t24.bacteria "-2.75±0.53" "-5.1809" "0.0000" "0.0000"
## t72.naive - t72.PBS "-0.46±0.54" "-0.8494" "0.3956" "0.3956"
## t72.naive - t72.bacteria "-2.05±0.54" "-3.7579" "0.0002" "0.0012"
## t24.naive - t72.naive " 1.40±0.54" " 2.5696" "0.0102" "0.0713"
## t24.PBS - t72.PBS " 2.35±0.53" " 4.4367" "0.0000" "0.0001"
## t24.bacteria - t72.bacteria " 2.10±0.53" " 3.9600" "0.0001" "0.0005"
##
## $PGRP
## estimate tvalue pvalue qvalue
## t24.naive - t24.PBS "-0.21±0.47" "-0.4413" "0.6590" "0.7689"
## t24.naive - t24.bacteria "-0.14±0.47" "-0.2946" "0.7683" "0.7683"
## t72.naive - t72.PBS " 0.68±0.48" " 1.4175" "0.1563" "0.2563"
## t72.naive - t72.bacteria "-0.68±0.48" "-1.4262" "0.1538" "0.1795"
## t24.naive - t72.naive "-0.24±0.48" "-0.4940" "0.6213" "0.7960"
## t24.PBS - t72.PBS " 0.65±0.47" " 1.3905" "0.1644" "0.2877"
## t24.bacteria - t72.bacteria "-0.78±0.47" "-1.6750" "0.0939" "0.1096"
##
## $proPO
## estimate tvalue pvalue qvalue
## t24.naive - t24.PBS "-0.63±0.49" "-1.2897" "0.1972" "0.2760"
## t24.naive - t24.bacteria "-1.77±0.46" "-3.8250" "0.0001" "0.0002"
## t72.naive - t72.PBS "-0.75±0.48" "-1.5655" "0.1175" "0.2563"
## t72.naive - t72.bacteria "-0.98±0.47" "-2.1019" "0.0356" "0.0622"
## t24.naive - t72.naive " 0.32±0.47" " 0.6721" "0.5015" "0.7960"
## t24.PBS - t72.PBS " 0.20±0.49" " 0.4006" "0.6887" "0.8035"
## t24.bacteria - t72.bacteria " 1.10±0.46" " 2.3747" "0.0176" "0.0307"
##
## $Serpin
## estimate tvalue pvalue qvalue
## t24.naive - t24.PBS " 0.84±0.38" " 2.2039" "0.0275" "0.0482"
## t24.naive - t24.bacteria "-1.76±0.38" "-4.6109" "0.0000" "0.0000"
## t72.naive - t72.PBS "-0.40±0.39" "-1.0224" "0.3066" "0.3577"
## t72.naive - t72.bacteria "-1.18±0.39" "-3.0078" "0.0026" "0.0092"
## t24.naive - t72.naive " 0.06±0.39" " 0.1528" "0.8786" "0.8786"
## t24.PBS - t72.PBS "-1.18±0.38" "-3.0974" "0.0020" "0.0068"
## t24.bacteria - t72.bacteria " 0.64±0.38" " 1.6776" "0.0934" "0.1096"#lapply(HE.GE.Fem.GxG.PH.Results, function(x) write.table( data.frame(x),
# "ResultsV3/HE.GE.Fem.GxG.PH.Results.csv", append=T, sep=',' ))6.1.4) Gene by gene; Plots
head(HE.GE.Fem.Mat)## ID FullTreatment Attacin B1.3g Histone Lysozyme
## 1 L0116-202.1 t24.bacteria -17.81554 -26.04675 -23.95141 -24.57171
## 2 L0116-220.1 t24.bacteria -19.21890 -24.95628 -24.69162 -23.24640
## 3 L0116-224.3 t24.bacteria -20.57175 -27.64199 -22.99736 -25.75324
## 4 L0116-290.3 t24.bacteria -16.89788 -25.59975 -24.16983 -24.61726
## 5 L0116-256.4 t24.bacteria -17.91622 -24.92581 -23.99452 -24.23819
## 6 L0116-311.5 t24.bacteria -17.60451 -23.76316 -24.58395 -24.76010
## mrpL37 Pelle PGRP proPO S28rpS24 Serpin Larval_id
## 1 -23.49135 -25.42495 -28.00394 -30.45286 -23.80937 -21.59232 L0116-202
## 2 -24.02148 -25.14932 -28.55690 -30.15535 -24.08870 -20.72778 L0116-220
## 3 -23.13258 -23.04757 -27.29577 -29.41526 -23.33178 -21.26012 L0116-224
## 4 -24.52375 -25.86108 -28.09212 -29.28932 -24.80637 -21.46396 L0116-290
## 5 -24.00403 -26.77086 -27.84653 -29.72143 -24.24882 -21.34987 L0116-256
## 6 -25.42397 -23.95025 -28.24660 -28.09622 -26.17418 -22.20284 L0116-311
## Family Timepoint Treatment
## 1 1 t24 bacteria
## 2 1 t24 bacteria
## 3 3 t24 bacteria
## 4 3 t24 bacteria
## 5 4 t24 bacteria
## 6 5 t24 bacteriaHE.GE.Fem.Res.means <- aggregate(HE.GE.Fem.Mat[,c(3,4,6,8,9,10,12)],by=list(HE.GE.Fem.Mat$Treatment,
HE.GE.Fem.Mat$Timepoint),FUN=mean,na.rm=TRUE)
HE.GE.Fem.Res.CIs <- aggregate(HE.GE.Fem.Mat[,c(3,4,6,8,9,10,12)],by=list(HE.GE.Fem.Mat$Treatment,
HE.GE.Fem.Mat$Timepoint),FUN=CI)
xpos <- seq(0,2,1)+rep(c(0,4),each=3)+0.5
#svg("ResultsV3/HE.GE.Fem.Plot.svg",width=20,height=10)
par(mfrow=c(3,3))
par(mar=c(2,6,1,1))
for(i in 3:9){
ylimi <- round(range(c(floor(HE.GE.Fem.Res.means[,i]+HE.GE.Fem.Res.CIs[,i]*-1),ceil(HE.GE.Fem.Res.means[,i]+HE.GE.Fem.Res.CIs[,i]))))
plot(NA,xlim=c(0,7),ylim=ylimi,axes=FALSE,xlab="",cex.lab=5,ylab="")
axis(side=1,at=c(-2,8),labels=c("",""),font=2)
axis(side=2,at=seq(ylimi[1]-2,ylimi[2]+1,1),labels=seq(ylimi[1]-2,ylimi[2]+1,1),
las=2,font=2,cex.axis=3)
axis(side=3,at=c(-2,8))
axis(side=4,at=c(ylimi[1]-2,ylimi[2]+1))
errbar(x=xpos,y=HE.GE.Fem.Res.means[,i],yplus=HE.GE.Fem.Res.means[,i]+HE.GE.Fem.Res.CIs[,i],yminus=HE.GE.Fem.Res.means[,i]-HE.GE.Fem.Res.CIs[,i],lwd=4,pch="",add=TRUE)
points(x=xpos,y=HE.GE.Fem.Res.means[,i],pch=c(21,22,23),col=c("black"),cex=4,lwd=4,bg=rep(c("black","grey"),each=3))
text(0,ylimi[1],labels=colnames(HE.GE.Fem.Res.means)[i],cex=2,adj=0)
}
#dev.off()
Version two; add all points
head(HE.GE.Fem.Mat)## ID FullTreatment Attacin B1.3g Histone Lysozyme
## 1 L0116-202.1 t24.bacteria -17.81554 -26.04675 -23.95141 -24.57171
## 2 L0116-220.1 t24.bacteria -19.21890 -24.95628 -24.69162 -23.24640
## 3 L0116-224.3 t24.bacteria -20.57175 -27.64199 -22.99736 -25.75324
## 4 L0116-290.3 t24.bacteria -16.89788 -25.59975 -24.16983 -24.61726
## 5 L0116-256.4 t24.bacteria -17.91622 -24.92581 -23.99452 -24.23819
## 6 L0116-311.5 t24.bacteria -17.60451 -23.76316 -24.58395 -24.76010
## mrpL37 Pelle PGRP proPO S28rpS24 Serpin Larval_id
## 1 -23.49135 -25.42495 -28.00394 -30.45286 -23.80937 -21.59232 L0116-202
## 2 -24.02148 -25.14932 -28.55690 -30.15535 -24.08870 -20.72778 L0116-220
## 3 -23.13258 -23.04757 -27.29577 -29.41526 -23.33178 -21.26012 L0116-224
## 4 -24.52375 -25.86108 -28.09212 -29.28932 -24.80637 -21.46396 L0116-290
## 5 -24.00403 -26.77086 -27.84653 -29.72143 -24.24882 -21.34987 L0116-256
## 6 -25.42397 -23.95025 -28.24660 -28.09622 -26.17418 -22.20284 L0116-311
## Family Timepoint Treatment
## 1 1 t24 bacteria
## 2 1 t24 bacteria
## 3 3 t24 bacteria
## 4 3 t24 bacteria
## 5 4 t24 bacteria
## 6 5 t24 bacteriaGeneCols <-c(3,4,6,8,9,10,12)
HE.GE.Fem.Res.means <- aggregate(HE.GE.Fem.Mat[,c(3,4,6,8,9,10,12)],by=list(HE.GE.Fem.Mat$Treatment,
HE.GE.Fem.Mat$Timepoint),FUN=mean,na.rm=TRUE)
HE.GE.Fem.Res.CIs <- aggregate(HE.GE.Fem.Mat[,c(3,4,6,8,9,10,12)],by=list(HE.GE.Fem.Mat$Treatment,
HE.GE.Fem.Mat$Timepoint),FUN=CI)
xpos <- seq(0,2,1)+rep(c(0,4),each=3)+0.5
#svg("ResultsV3/HE.GE.Fem.Plot_Points.svg",width=20,height=10)
par(mfrow=c(3,3))
par(mar=c(2,6,1,1))
for(i in 3:9){
ylimi <- round(range(c(floor(HE.GE.Fem.Mat[,GeneCols[i-2]]),ceiling(HE.GE.Fem.Mat[,GeneCols[i-2]]))))
plot(NA,xlim=c(0,7),ylim=ylimi,axes=FALSE,xlab="",cex.lab=5,ylab="")
axis(side=1,at=c(-2,8),labels=c("",""),font=2)
axis(side=2,at=seq(ylimi[1]-2,ylimi[2]+1,1),labels=seq(ylimi[1]-2,ylimi[2]+1,1),
las=2,font=2,cex.axis=3)
axis(side=3,at=c(-2,8))
axis(side=4,at=c(ylimi[1]-2,ylimi[2]+1))
errbar(x=xpos,y=HE.GE.Fem.Res.means[,i],yplus=HE.GE.Fem.Res.means[,i]+HE.GE.Fem.Res.CIs[,i],yminus=HE.GE.Fem.Res.means[,i]-HE.GE.Fem.Res.CIs[,i],lwd=4,pch="",add=TRUE)
points(x=xpos,y=HE.GE.Fem.Res.means[,i],pch=c(21,22,23),col=c("black"),cex=4,lwd=4,bg=rep(c("black","gray"),each=3))
text(0,ylimi[1],labels=colnames(HE.GE.Fem.Res.means)[i],cex=2,adj=0)
points(x=rep(xpos,tapply(HE.GE.Fem.Mat[,GeneCols[i-2]],list(HE.GE.Fem.Mat$Treatment,HE.GE.Fem.Mat$Timepoint),length)),HE.GE.Fem.Mat[order(HE.GE.Fem.Mat$Timepoint,HE.GE.Fem.Mat$Treatment),GeneCols[i-2]],pch=4,col="gray32",cex=1.8,lwd=1.3)
}
#dev.off()
6.2) Encapsulation
6.2.2) Analysis
HE.model.EN.Full <- lmer(Corr_encaps~Treatment+(1|Family),data=HE.EN.Fem)
## Linearity
plot(resid(HE.model.EN.Full),fitted(HE.model.EN.Full))
## Homogeneity
plot(HE.model.EN.Full)
## Normality
qqnorm(residuals(HE.model.EN.Full))
## Influence
infl <- influence(HE.model.EN.Full, obs = TRUE)
plot(infl, which = "cook")
- Linearity is split, but otherwise fine.
- Homogeneity is fine.
- The qq-line is fine.
- No notable outliers.
Get the results
Display the model summary, and calculate the means of the gray values for each treatment.
As we are not interested in the comparison between the concentrations, the model summary provides sufficient information.
summary(HE.model.EN.Full)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Corr_encaps ~ Treatment + (1 | Family)
## Data: HE.EN.Fem
##
## REML criterion at convergence: 628.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.39103 -0.57882 -0.00857 0.70898 1.93836
##
## Random effects:
## Groups Name Variance Std.Dev.
## Family (Intercept) 6.949 2.636
## Residual 105.123 10.253
## Number of obs: 85, groups: Family, 11
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.5455 2.1600 22.0199 2.567 0.0176 *
## TreatmentPBS 0.7261 2.7727 78.0388 0.262 0.7941
## TreatmentBacterial -15.9691 2.7464 78.6447 -5.815 1.24e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtPBS
## TreatmntPBS -0.630
## TrtmntBctrl -0.657 0.501tapply(HE.EN.Fem$Corr_encaps,HE.EN.Fem$Treatment,mean,na.rm=TRUE)## Control PBS Bacterial
## 5.164186 6.238182 -10.550403tapply(HE.EN.Fem$Corr_encaps,HE.EN.Fem$Treatment,se)## Control PBS Bacterial
## 1.690743 1.934847 2.2458036.3) Life histories
6.3.1) Analysis
Prepare
Set the first clutch to be represented as zero. This way the intercept can be interpreted as the average difference in each life history trait at the first clutch.
Only PBS has clutch numbers larger than 4 (i.e. more than five clutches). I remove these, as PBS by itself is not of primary interest and the additional clutches are useless in the comparison to the control or bacteria treatments.
HE.RP$Clutch_rank0 <- HE.RP$Clutch_rank-1
HE.RP.Red <- HE.RP[-which(HE.RP$Clutch_rank0>=4),]6.3.2) Laying time
How does tratement and clutch number affect the timing of the clutches?
HE.model.RP.LayTime.Full <- lmer(yday(Laying_date)~Treatment*Clutch_rank0+
(1|Female_id),data=HE.RP.Red)
## Linearity
plot(resid(HE.model.RP.LayTime.Full),fitted(HE.model.RP.LayTime.Full))
## Homogeneity
plot(HE.model.RP.LayTime.Full)
## Normality
qqnorm(residuals(HE.model.RP.LayTime.Full))
## Influence
infl <- influence(HE.model.RP.LayTime.Full, obs = TRUE)
plot(infl, which = "cook")
- Linearity is fine.
- Homogeneity is slightly off.
- The qq-line is fine.
- No notable outliers.
Get the results
Display the model summary.
Get the relationship between clutch rank and egg laying time for each treatment separately.
summary(HE.model.RP.LayTime.Full)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: yday(Laying_date) ~ Treatment * Clutch_rank0 + (1 | Female_id)
## Data: HE.RP.Red
##
## REML criterion at convergence: 606.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8487 -0.5392 -0.1198 0.5262 2.1505
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female_id (Intercept) 43.566 6.600
## Residual 7.147 2.673
## Number of obs: 107, groups: Female_id, 40
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 84.37401 2.12226 37.71016 39.757
## TreatmentPBS 0.05302 2.71952 37.48006 0.019
## TreatmentBacteria 0.54894 2.93593 37.59424 0.187
## Clutch_rank0 6.75805 0.56840 64.77619 11.890
## TreatmentPBS:Clutch_rank0 -1.84567 0.70127 64.42882 -2.632
## TreatmentBacteria:Clutch_rank0 -2.32760 0.74042 64.50649 -3.144
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## TreatmentPBS 0.98455
## TreatmentBacteria 0.85269
## Clutch_rank0 < 2e-16 ***
## TreatmentPBS:Clutch_rank0 0.01061 *
## TreatmentBacteria:Clutch_rank0 0.00252 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtPBS TrtmnB Cltc_0 TPBS:C
## TreatmntPBS -0.780
## TretmntBctr -0.723 0.564
## Clutch_rnk0 -0.203 0.158 0.146
## TrtmPBS:C_0 0.164 -0.208 -0.119 -0.811
## TrtmntB:C_0 0.156 -0.121 -0.206 -0.768 0.622emtrends(HE.model.RP.LayTime.Full,pairwise~Treatment,var="Clutch_rank0")## $emtrends
## Treatment Clutch_rank0.trend SE df lower.CL upper.CL
## Control 6.76 0.568 64.8 5.62 7.89
## PBS 4.91 0.411 63.8 4.09 5.73
## Bacteria 4.43 0.474 64.1 3.48 5.38
##
## Trends are based on the yday (transformed) scale
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Control - PBS 1.846 0.701 64.4 2.632 0.0282
## Control - Bacteria 2.328 0.740 64.5 3.144 0.0070
## PBS - Bacteria 0.482 0.628 64.0 0.768 0.7239
##
## P value adjustment: tukey method for comparing a family of 3 estimatesaggregate(yday(HE.RP.Red$Laying_date), by=list(HE.RP.Red$Clutch_rank,HE.RP.Red$Treatment),FUN=mean)## Group.1 Group.2 x
## 1 1 Control 84.45455
## 2 2 Control 87.57143
## 3 3 Control 95.50000
## 4 4 Control 98.33333
## 5 1 PBS 84.29412
## 6 2 PBS 87.92857
## 7 3 PBS 94.33333
## 8 4 PBS 95.00000
## 9 1 Bacteria 85.16667
## 10 2 Bacteria 84.87500
## 11 3 Bacteria 89.62500
## 12 4 Bacteria 94.00000aggregate(yday(HE.RP.Red$Laying_date), by=list(HE.RP.Red$Clutch_rank,HE.RP.Red$Treatment),FUN=se)## Group.1 Group.2 x
## 1 1 Control 2.1547564
## 2 2 Control 1.7300859
## 3 3 Control 1.8574176
## 4 4 Control 2.3333333
## 5 1 PBS 1.5547950
## 6 2 PBS 1.0967699
## 7 3 PBS 2.5331140
## 8 4 PBS 2.3502786
## 9 1 Bacteria 2.1979099
## 10 2 Bacteria 0.8543481
## 11 3 Bacteria 1.6468314
## 12 4 Bacteria 3.0166206Plot
Plot the data and add a prediction line including 95% confidence interval.
newdat <- data.frame(Clutch_rank0=rep(0:3,3),Treatment=rep(c("Control","PBS","Bacteria"),each=4))
pfun <- function(fit) {
predict(fit,newdata=newdat,re.form=NA)
}
prediction <- cbind(newdat,pred=predict(HE.model.RP.LayTime.Full,newdata=newdat,re.form=NA))
prediction <- cbind(prediction,confint(HE.model.RP.LayTime.Full,method="boot",FUN=pfun))
colnames(prediction)[c(4,5)] <- c("lwr","upr")#svg("ResultsV3/HE.RP.LayDay.Plot.svg",width=10,height=10)
ylimi <- range(yday(HE.RP.Red$Laying_date),na.rm=TRUE)
plot(yday(HE.RP.Red$Laying_date)~jitter(HE.RP.Red$Clutch_rank0,0.6),subset=HE.RP.Red$Treatment=="Bacteria",ylim=ylimi,pch=16,col=2,xlab="Clutch number",ylab="Laying day",las=1,xaxt="n")
points(yday(HE.RP.Red$Laying_date)~jitter(HE.RP.Red$Clutch_rank0,0.6),subset=HE.RP.Red$Treatment=="Control",pch=16,col=1)
points(yday(HE.RP.Red$Laying_date)~jitter(HE.RP.Red$Clutch_rank0,0.6),subset=HE.RP.Red$Treatment=="PBS",pch=16,col=4)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Bacteria",col=2)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Control",col=1)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="PBS",col=4)
with(prediction[prediction$Treatment=="Control",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("black", alpha.f = 0.10), border = NA))
with(prediction[prediction$Treatment=="Bacteria",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("red", alpha.f = 0.10), border = NA))
with(prediction[prediction$Treatment=="PBS",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("blue", alpha.f = 0.10), border = NA))
legend("bottomright",pch=16,col=c(1,2,4),legend=c("Control","PBS","5mg/ml"))
axis(1,at=c(-1:5))
#dev.off()6.3.3) Number of clutches
How many clutches do the females of the different treatments lay in total?
HE.model.RP.NrClutch.Full <- glm(x~Group.2,data=HE.CPF,family=poisson())
## Overdispersion
overdisp_fun(HE.model.RP.NrClutch.Full)## chisq ratio rdf p
## 27.5117845 0.7435617 37.0000000 0.8718356## Linearity
plot(resid(HE.model.RP.NrClutch.Full),fitted(HE.model.RP.NrClutch.Full))
## Homogeneity
plot(HE.model.RP.NrClutch.Full)
## Normality
qqnorm(residuals(HE.model.RP.NrClutch.Full))
## Influence
#infl <- influence(HE.model.RP.NrClutch.Full, obs = TRUE)
#plot(infl, which = "cook")- Dispersion is acceptable.
- Linearity is fine.
- Homogeneity is fine.
- The qq-line looks slightly off. But it’s okay for analysis.
- No notable outliers.
Get the results
Display the model summary.
Get the average number of clutches produced in each treatment.
summary(HE.model.RP.NrClutch.Full)##
## Call:
## glm(formula = x ~ Group.2, family = poisson(), data = HE.CPF)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3427 -1.0551 0.1485 0.5881 1.5224
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.8979 0.1924 4.666 3.07e-06 ***
## Group.2PBS 0.2007 0.2380 0.843 0.399
## Group.2Bacteria 0.1137 0.2595 0.438 0.661
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 29.961 on 39 degrees of freedom
## Residual deviance: 29.233 on 37 degrees of freedom
## AIC: 146.62
##
## Number of Fisher Scoring iterations: 4tapply(HE.CPF$x,HE.CPF$Group.2,mean)## Control PBS Bacteria
## 2.454545 3.000000 2.750000tapply(HE.CPF$x,HE.CPF$Group.2,se)## Control PBS Bacteria
## 0.3899566 0.3834825 0.39167475.3.4) Number of eggs per clutch
How many eggs do the females lay in each clutch?
HE.model.RP.NrEgg.Full <- glmer(Number_eggs~Treatment*Clutch_rank0+(1|Female_id),
data=HE.RP.Red,family=poisson())
## Overdispersion
overdisp_fun(HE.model.RP.NrEgg.Full)## chisq ratio rdf p
## 2285.87231 22.85872 100.00000 0.00000## Linearity
plot(resid(HE.model.RP.NrEgg.Full),fitted(HE.model.RP.NrEgg.Full))
## Homogeneity
plot(HE.model.RP.NrEgg.Full)
## Normality
qqnorm(residuals(HE.model.RP.NrEgg.Full))
## Influence
infl <- influence(HE.model.RP.NrEgg.Full, obs = TRUE)
plot(infl, which = "cook")
The data is overdispersed. Use a negative binomial distribution.
HE.model.RP.NrEgg.Full <- glmer.nb(Number_eggs~Treatment*Clutch_rank0+
(1|Female_id),data=HE.RP.Red)
## Overdispersion
overdisp_fun(HE.model.RP.NrEgg.Full)## chisq ratio rdf p
## 81.0942441 0.8191338 99.0000000 0.9049601## Linearity
plot(resid(HE.model.RP.NrEgg.Full),fitted(HE.model.RP.NrEgg.Full))
## Homogeneity
plot(HE.model.RP.NrEgg.Full)
## Normality
qqnorm(residuals(HE.model.RP.NrEgg.Full))
## Influence
infl <- influence(HE.model.RP.NrEgg.Full, obs = TRUE)
plot(infl, which = "cook")
- Overdispersion: fixed.
- Homogeneity is fine.
- No notable outliers.
Get the results
Display the model summary.
Get the relationship between clutch rank and number eggs for each treatment separately.
Display the mean number of eggs per clutch for each treatment.
summary(HE.model.RP.NrEgg.Full)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(2.6044) ( log )
## Formula: Number_eggs ~ Treatment * Clutch_rank0 + (1 | Female_id)
## Data: HE.RP.Red
##
## AIC BIC logLik deviance df.resid
## 1218.0 1239.4 -601.0 1202.0 99
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.54258 -0.68058 -0.00046 0.56370 2.71638
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female_id (Intercept) 2.566e-14 1.602e-07
## Number of obs: 107, groups: Female_id, 40
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.93554 0.17585 28.067 <2e-16 ***
## TreatmentPBS 0.09020 0.22242 0.406 0.685
## TreatmentBacteria -0.12234 0.24163 -0.506 0.613
## Clutch_rank0 -0.03176 0.12366 -0.257 0.797
## TreatmentPBS:Clutch_rank0 -0.21323 0.15293 -1.394 0.163
## TreatmentBacteria:Clutch_rank0 0.01950 0.16263 0.120 0.905
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtPBS TrtmnB Cltc_0 TPBS:C
## TreatmntPBS -0.791
## TretmntBctr -0.728 0.575
## Clutch_rnk0 -0.729 0.576 0.530
## TrtmPBS:C_0 0.589 -0.733 -0.429 -0.809
## TrtmntB:C_0 0.554 -0.438 -0.739 -0.760 0.615
## convergence code: 0
## boundary (singular) fit: see ?isSingularemtrends(HE.model.RP.NrEgg.Full,pairwise~Treatment,var="Clutch_rank0")## $emtrends
## Treatment Clutch_rank0.trend SE df asymp.LCL asymp.UCL
## Control -0.0318 0.124 Inf -0.274 0.2106
## PBS -0.2450 0.090 Inf -0.421 -0.0686
## Bacteria -0.0123 0.106 Inf -0.219 0.1948
##
## Trends are based on the log (transformed) scale
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## Control - PBS 0.2132 0.153 Inf 1.394 0.3439
## Control - Bacteria -0.0195 0.163 Inf -0.120 0.9921
## PBS - Bacteria -0.2327 0.139 Inf -1.677 0.2140
##
## P value adjustment: tukey method for comparing a family of 3 estimatesaggregate(HE.RP.Red$Number_eggs, by=list(HE.RP.Red$Clutch_rank,HE.RP.Red$Treatment),FUN=mean)## Group.1 Group.2 x
## 1 1 Control 125.54545
## 2 2 Control 164.00000
## 3 3 Control 134.83333
## 4 4 Control 99.66667
## 5 1 PBS 137.52941
## 6 2 PBS 141.71429
## 7 3 PBS 89.88889
## 8 4 PBS 66.00000
## 9 1 Bacteria 111.00000
## 10 2 Bacteria 146.37500
## 11 3 Bacteria 124.50000
## 12 4 Bacteria 101.20000aggregate(HE.RP.Red$Number_eggs, by=list(HE.RP.Red$Clutch_rank,HE.RP.Red$Treatment),FUN=se)## Group.1 Group.2 x
## 1 1 Control 16.35591
## 2 2 Control 23.94040
## 3 3 Control 32.38046
## 4 4 Control 53.02934
## 5 1 PBS 16.49936
## 6 2 PBS 23.82363
## 7 3 PBS 19.07110
## 8 4 PBS 21.24909
## 9 1 Bacteria 14.74377
## 10 2 Bacteria 19.62409
## 11 3 Bacteria 24.75451
## 12 4 Bacteria 32.21397Plot
Plot the data and add a prediction line including 95% confidence interval.
newdat <- data.frame(Clutch_rank0=rep(0:3,3),Treatment=rep(c("Control","PBS","Bacteria"),each=4))
pfun <- function(fit) {
predict(fit,newdata=newdat,re.form=NA,type="response")
}
prediction <- cbind(newdat,pred=predict(HE.model.RP.NrEgg.Full,newdata=newdat,re.form=NA,type="response"))
prediction <- cbind(prediction,confint(HE.model.RP.NrEgg.Full,method="boot",FUN=pfun))
colnames(prediction)[c(4,5)] <- c("lwr","upr")#svg("ResultsV3/HE.RP.NrEgg.Plot.svg",width=10,height=10)
ylimi <- range(HE.RP.Red$Number_eggs,na.rm=TRUE)
plot(HE.RP.Red$Number_eggs~jitter(HE.RP.Red$Clutch_rank0,0.6),subset=HE.RP.Red$Treatment=="Bacteria",pch=16,col=2,ylim=ylimi,xlab="Clutch number",ylab="Number eggs",las=1,xaxt="n")
points(HE.RP.Red$Number_eggs~jitter(HE.RP.Red$Clutch_rank0,0.6),subset=HE.RP.Red$Treatment=="Control",pch=16,col=1)
points(HE.RP.Red$Number_eggs~jitter(HE.RP.Red$Clutch_rank0,0.6),subset=HE.RP.Red$Treatment=="PBS",pch=16,col=4)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Bacteria",col=2)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Control",col=1)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="PBS",col=4)
with(prediction[prediction$Treatment=="Control",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("black", alpha.f = 0.10), border = NA))
with(prediction[prediction$Treatment=="Bacteria",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("red", alpha.f = 0.10), border = NA))
with(prediction[prediction$Treatment=="PBS",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("blue", alpha.f = 0.10), border = NA))
legend("topright",pch=16,col=c(1,4,2),legend=c("Control","PBS","5mg/ml"))
axis(1,at=c(-1:5))
#dev.off()6.3.5) Lifetime egg production
How does lifetime egg production compare between the two treatments?
Given the low sample size, the outcome will not be very reliable.
tapply(HE.LEP$x,HE.LEP$Group.2,mean)## Control PBS Bacteria
## 330.6364 352.4118 333.7500tapply(HE.LEP$x,HE.LEP$Group.2,median)## Control PBS Bacteria
## 300.0 303.0 350.5tapply(HE.LEP$x,HE.LEP$Group.2,se)## Control PBS Bacteria
## 73.70473 47.66257 66.64278ci.calc <- function(x){boot.ci(boot(x,function(x,i) median(x[i]), R=1000))}
tapply(HE.LEP$x,HE.LEP$Group.2,ci.calc)## $Control
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot(x, function(x, i) median(x[i]), R = 1000))
##
## Intervals :
## Level Normal Basic
## 95% (137.9, 490.7 ) (102.0, 507.0 )
##
## Level Percentile BCa
## 95% ( 93, 498 ) ( 84, 383 )
## Calculations and Intervals on Original Scale
## Some BCa intervals may be unstable
##
## $PBS
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot(x, function(x, i) median(x[i]), R = 1000))
##
## Intervals :
## Level Normal Basic
## 95% (111.3, 418.9 ) ( 90.0, 367.0 )
##
## Level Percentile BCa
## 95% (239, 516 ) (237, 450 )
## Calculations and Intervals on Original Scale
## Some BCa intervals may be unstable
##
## $Bacteria
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot(x, function(x, i) median(x[i]), R = 1000))
##
## Intervals :
## Level Normal Basic
## 95% (163.2, 558.3 ) (226.5, 580.5 )
##
## Level Percentile BCa
## 95% (120.5, 474.5 ) ( 83.0, 466.0 )
## Calculations and Intervals on Original Scalesummary.aov(lm(HE.LEP$x~HE.LEP$Group.2))## Df Sum Sq Mean Sq F value Pr(>F)
## HE.LEP$Group.2 2 4025 2012 0.041 0.96
## Residuals 37 1801717 48695plot(HE.LEP$x~HE.LEP$Group.2)
6.3.5) Average egg weigt
What’s the average weight ef the eggs in each clutch and treatment?
HE.model.RP.Weight.Full <- lmer(Weight_per_egg~Treatment*Clutch_rank0+
(1|Female_id),data=HE.RP.Red)
## Linearity
plot(resid(HE.model.RP.Weight.Full),fitted(HE.model.RP.Weight.Full))
## Homogeneity
plot(HE.model.RP.Weight.Full)
## Normality
qqnorm(residuals(HE.model.RP.Weight.Full))
## Influence
infl <- influence(HE.model.RP.Weight.Full, obs = TRUE)
plot(infl, which = "cook")
- Linearity could be better, possible outliers?.
- Homogeneity could be better.
- The qq-line is straight except for two extreme points.
- Three possible outliers, but the influence is acceptable (Cook’s D < 1).
This is borderline. I think the violations are just acceptable for analysis. Keep it as it is.
summary(HE.model.RP.Weight.Full)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Weight_per_egg ~ Treatment * Clutch_rank0 + (1 | Female_id)
## Data: HE.RP.Red
##
## REML criterion at convergence: -519.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.99515 -0.44380 0.05967 0.59661 1.74575
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female_id (Intercept) 5.589e-05 0.007476
## Residual 5.407e-05 0.007353
## Number of obs: 89, groups: Female_id, 38
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.085532 0.003160 47.590881 27.068
## TreatmentPBS 0.004768 0.004033 45.246724 1.182
## TreatmentBacteria 0.005064 0.004366 44.523148 1.160
## Clutch_rank0 -0.005111 0.001693 56.294187 -3.019
## TreatmentPBS:Clutch_rank0 -0.001752 0.002075 55.560595 -0.844
## TreatmentBacteria:Clutch_rank0 -0.001082 0.002230 54.049153 -0.485
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## TreatmentPBS 0.24329
## TreatmentBacteria 0.25232
## Clutch_rank0 0.00381 **
## TreatmentPBS:Clutch_rank0 0.40216
## TreatmentBacteria:Clutch_rank0 0.62954
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtPBS TrtmnB Cltc_0 TPBS:C
## TreatmntPBS -0.783
## TretmntBctr -0.724 0.567
## Clutch_rnk0 -0.443 0.347 0.321
## TrtmPBS:C_0 0.361 -0.430 -0.262 -0.816
## TrtmntB:C_0 0.336 -0.263 -0.441 -0.759 0.619emtrends(HE.model.RP.Weight.Full,pairwise~Treatment,var="Clutch_rank0")## $emtrends
## Treatment Clutch_rank0.trend SE df lower.CL upper.CL
## Control -0.00511 0.00169 56.3 -0.00850 -0.00172
## PBS -0.00686 0.00120 54.1 -0.00927 -0.00446
## Bacteria -0.00619 0.00145 51.1 -0.00911 -0.00328
##
## Degrees-of-freedom method: satterthwaite
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Control - PBS 0.00175 0.00208 55.6 0.844 0.6773
## Control - Bacteria 0.00108 0.00223 54.0 0.485 0.8786
## PBS - Bacteria -0.00067 0.00188 52.3 -0.356 0.9327
##
## P value adjustment: tukey method for comparing a family of 3 estimatesaggregate(HE.RP.Red$Weight_per_egg, by=list(HE.RP.Red$Clutch_rank,HE.RP.Red$Treatment),FUN=mean,na.rm=TRUE)## Group.1 Group.2 x
## 1 1 Control 0.08347796
## 2 2 Control 0.08179785
## 3 3 Control 0.08166831
## 4 4 Control 0.06360635
## 5 1 PBS 0.09049094
## 6 2 PBS 0.08552898
## 7 3 PBS 0.07360705
## 8 4 PBS 0.07556703
## 9 1 Bacteria 0.09064904
## 10 2 Bacteria 0.08501371
## 11 3 Bacteria 0.07942519
## 12 4 Bacteria 0.07544719aggregate(HE.RP.Red$Weight_per_egg, by=list(HE.RP.Red$Clutch_rank,HE.RP.Red$Treatment),FUN=se)## Group.1 Group.2 x
## 1 1 Control 0.0028780725
## 2 2 Control 0.0007386683
## 3 3 Control 0.0058234331
## 4 4 Control 0.0077998992
## 5 1 PBS 0.0025443736
## 6 2 PBS 0.0021443211
## 7 3 PBS 0.0072844559
## 8 4 PBS 0.0064722856
## 9 1 Bacteria 0.0026831925
## 10 2 Bacteria 0.0021943711
## 11 3 Bacteria 0.0025631482
## 12 4 Bacteria 0.0055903937Plot
Plot the data and add a prediction line including 95% confidence interval.
newdat <- data.frame(Clutch_rank0=rep(0:6,3),Treatment=rep(c("Control","PBS","Bacteria"),each=7))
pfun <- function(fit) {
predict(fit,newdata=newdat,re.form=NA,type="response")
}
prediction <- cbind(newdat,pred=predict(HE.model.RP.Weight.Full,newdata=newdat,re.form=NA,type="response"))
prediction <- cbind(prediction,confint(HE.model.RP.Weight.Full,method="boot",FUN=pfun))
colnames(prediction)[c(4,5)] <- c("lwr","upr")#svg("ResultsV3/HE.RP.Weight.Plot.svg",width=10,height=10)
ylimi <- range(HE.RP.Red$Weight_per_egg,na.rm=TRUE)
plot(HE.RP.Red$Weight_per_egg~jitter(HE.RP.Red$Clutch_rank0,0.6),subset=HE.RP.Red$Treatment=="Bacteria",pch=16,col=2,ylim=ylimi,xlab="Clutch number",ylab="Average egg weight (ng)",las=1,xaxt="n")
points(HE.RP.Red$Weight_per_egg~jitter(HE.RP.Red$Clutch_rank0,0.6),subset=HE.RP.Red$Treatment=="Control",pch=16,col=1)
points(HE.RP.Red$Weight_per_egg~jitter(HE.RP.Red$Clutch_rank0,0.6),subset=HE.RP.Red$Treatment=="PBS",pch=16,col=4)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Bacteria",col=2)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Control",col=1)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="PBS",col=4)
with(prediction[prediction$Treatment=="Control",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("black", alpha.f = 0.10), border = NA))
with(prediction[prediction$Treatment=="Bacteria",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("red", alpha.f = 0.10), border = NA))
with(prediction[prediction$Treatment=="PBS",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("blue", alpha.f = 0.10), border = NA))
legend("topright",pch=16,col=c(1,2,4),legend=c("Control","PBS","5mg/ml"))
axis(1,at=c(-1:5))
#dev.off()5.3.6) Hatching success
How well do the larvae of eatch clutch hatch in each treatment?
HE.RP.Red.NoNA <- na.omit(HE.RP.Red[,c(2,5,8,9,14,15)]) # remove the NAs for the prediction
HE.model.RP.Hatch.Full <- glmer(cbind(Number_hatched,(Number_eggs-Number_hatched))~
Treatment*Clutch_rank0+(1|Female_id),
data=HE.RP.Red.NoNA,family=binomial())
## Homogeneity
plot(HE.model.RP.Hatch.Full)
## Normality
qqnorm(residuals(HE.model.RP.Hatch.Full))
## Influence
infl <- influence(HE.model.RP.Hatch.Full, obs = TRUE)
plot(infl, which = "cook")
- Homogeneity is fine.
- The qq line is fine.
- One severe outlier, others possible.
Check
HE.RP.Red.NoNA[c(30,59,62),]## Female_id Number_eggs Number_hatched Treatment Hatching_percentage
## 44 96 204 74 Control 36.27451
## 83 155 208 148 Bacteria 71.15385
## 87 158 237 231 PBS 97.46835
## Clutch_rank0
## 44 3
## 83 0
## 87 0None of the data points have a specific note. Keep them in.
Get the results
Display the model summary.
Get the relationship between clutch rank and hatching success for each treatment separately.
Display the mean hatching success per clutch for each treatment.
summary(HE.model.RP.Hatch.Full)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## cbind(Number_hatched, (Number_eggs - Number_hatched)) ~ Treatment *
## Clutch_rank0 + (1 | Female_id)
## Data: HE.RP.Red.NoNA
##
## AIC BIC logLik deviance df.resid
## 1442.3 1459.0 -714.2 1428.3 73
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.6478 -1.9149 -0.0575 2.1125 10.3542
##
## Random effects:
## Groups Name Variance Std.Dev.
## Female_id (Intercept) 2.844 1.687
## Number of obs: 80, groups: Female_id, 34
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.87364 0.53817 1.623 0.10452
## TreatmentPBS -0.23707 0.69471 -0.341 0.73292
## TreatmentBacteria 0.28702 0.78931 0.364 0.71614
## Clutch_rank0 -0.77569 0.06074 -12.770 < 2e-16 ***
## TreatmentPBS:Clutch_rank0 0.24088 0.07953 3.029 0.00246 **
## TreatmentBacteria:Clutch_rank0 -0.25140 0.08075 -3.113 0.00185 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtPBS TrtmnB Cltc_0 TPBS:C
## TreatmntPBS -0.775
## TretmntBctr -0.682 0.528
## Clutch_rnk0 -0.066 0.051 0.045
## TrtmPBS:C_0 0.050 -0.070 -0.034 -0.764
## TrtmntB:C_0 0.049 -0.038 -0.090 -0.752 0.574emtrends(HE.model.RP.Hatch.Full,pairwise~Treatment,var="Clutch_rank0",transform="response")## $emtrends
## Treatment Clutch_rank0.trend SE df asymp.LCL asymp.UCL
## Control -0.194 0.0152 Inf -0.224 -0.164
## PBS -0.134 0.0128 Inf -0.159 -0.108
## Bacteria -0.257 0.0133 Inf -0.283 -0.230
##
## Trends are obtained after back-transforming from the logit scale
## Results are given on the ( (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## Control - PBS -0.0603 0.0199 Inf -3.033 0.0068
## Control - Bacteria 0.0627 0.0202 Inf 3.109 0.0053
## PBS - Bacteria 0.1230 0.0185 Inf 6.658 <.0001
##
## Results are given on the ( (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimatesaggregate(HE.RP.Red.NoNA$Hatching_percentage, by=list(HE.RP.Red.NoNA$Clutch_rank,HE.RP.Red.NoNA$Treatment),FUN=mean, na.rm=TRUE)## Group.1 Group.2 x
## 1 0 Control 69.08693
## 2 1 Control 53.38117
## 3 2 Control 52.61996
## 4 3 Control 31.04048
## 5 0 PBS 63.42643
## 6 1 PBS 49.43677
## 7 2 PBS 67.55406
## 8 3 PBS 34.24426
## 9 0 Bacteria 75.33376
## 10 1 Bacteria 74.45500
## 11 2 Bacteria 41.81914
## 12 3 Bacteria 61.22144aggregate(HE.RP.Red.NoNA$Hatching_percentage, by=list(HE.RP.Red.NoNA$Clutch_rank,HE.RP.Red.NoNA$Treatment),FUN=se)## Group.1 Group.2 x
## 1 0 Control 6.934919
## 2 1 Control 15.362355
## 3 2 Control 10.628423
## 4 3 Control 5.234029
## 5 0 PBS 9.376832
## 6 1 PBS 10.141073
## 7 2 PBS 11.048047
## 8 3 PBS 9.379078
## 9 0 Bacteria 10.336117
## 10 1 Bacteria 9.879232
## 11 2 Bacteria 13.708610
## 12 3 Bacteria 22.528986Plot
Plot the data and add a prediction line including 95% confidence interval.
newdat <- data.frame(Clutch_rank0=rep(0:6,3),Treatment=rep(c("Control","PBS","Bacteria"),each=7))
pfun <- function(fit) {
predict(fit,newdata=newdat,re.form=NA,type="response")
}
prediction <- cbind(newdat,pred=predict(HE.model.RP.Hatch.Full,newdata=newdat,re.form=NA,type="response"))
prediction <- cbind(prediction,confint(HE.model.RP.Hatch.Full,method="boot",FUN=pfun))
colnames(prediction)[c(4,5)] <- c("lwr","upr")
prediction[,3:5] <- prediction[,3:5]*100#svg("ResultsV3/HE.RP.Hatch.Plot.svg",width=10,height=10)
ylimi <- range(HE.RP.Red$Hatching_percentage,na.rm=TRUE)
plot(HE.RP.Red.NoNA$Hatching_percentage~jitter(HE.RP.Red.NoNA$Clutch_rank0,0.6),subset=HE.RP.Red.NoNA$Treatment=="Bacteria",pch=16,col=2,ylim=ylimi,xlab="Clutch number",ylab="Hatching success (%)",las=1,xaxt="n")
points(HE.RP.Red.NoNA$Hatching_percentage~jitter(HE.RP.Red.NoNA$Clutch_rank0,0.6),subset=HE.RP.Red.NoNA$Treatment=="Control",pch=16,col=1)
points(HE.RP.Red.NoNA$Hatching_percentage~jitter(HE.RP.Red.NoNA$Clutch_rank0,0.6),subset=HE.RP.Red.NoNA$Treatment=="PBS",pch=16,col=4)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Bacteria",col=2)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="Control",col=1)
lines(prediction$pred~prediction$Clutch_rank0,subset=prediction$Treatment=="PBS",col=4)
with(prediction[prediction$Treatment=="Control",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("black", alpha.f = 0.10), border = NA))
with(prediction[prediction$Treatment=="Bacteria",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("red", alpha.f = 0.10), border = NA))
with(prediction[prediction$Treatment=="PBS",],
polygon(x=c(Clutch_rank0,rev(Clutch_rank0)),
y=c(lwr, rev(upr)),col = adjustcolor("blue", alpha.f = 0.10), border = NA))
legend("bottom",pch=16,col=c(1,2,4),legend=c("Control","PBS","5mg/ml"))
axis(1,at=c(-1:5))
#dev.off()6.4) Lifespan
6.4.1) Analyse
Check the proportional hazard assumption.
HE.model.LS.Full <- coxph(Surv(Day,Censor)~Injection_treatment+frailty(Family),data=HE.LS.Fem)
cox.zph(HE.model.LS.Full)## rho chisq p
## Injection_treatmentinjection -0.0963 0.416 0.519Run the model and get the results.
HE.model.LS.Full <- coxme(Surv(Day,Censor)~Injection_treatment+(1|Family),data=HE.LS.Fem)
summary(HE.model.LS.Full)## Cox mixed-effects model fit by maximum likelihood
## Data: HE.LS.Fem
## events, n = 48, 48 (8 observations deleted due to missingness)
## Iterations= 2 12
## NULL Integrated Fitted
## Log-likelihood -140.6739 -136.6063 -136.5894
##
## Chisq df p AIC BIC
## Integrated loglik 8.14 2.00 0.0171180 4.14 0.39
## Penalized loglik 8.17 1.02 0.0043869 6.14 4.23
##
## Model: Surv(Day, Censor) ~ Injection_treatment + (1 | Family)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## Injection_treatmentinjection 0.8908534 2.437209 0.3096367 2.88 0.004
##
## Random effects
## Group Variable Std Dev Variance
## Family Intercept 0.0199995211 0.0003999808Plot
#svg("ResultsV3/HE.LS.Plot.svg",width=10,height=10)
xlimi <- c(0,max(HE.LS.Fem$Day,na.rm=TRUE))
modsum <- summary(survfit(Surv(Day,Censor)~Injection_treatment+(1|Family),data=HE.LS.Fem),times=0:50)
with(modsum,plot(time,surv,type="n",xlim=xlimi,ylim=c(0,1),las=1,xlab="Days",ylab="Cumulative survival"))
polygon(c(modsum[["time"]][1:42],rev(modsum[["time"]][1:42])),c(modsum[["lower"]][1:42],rev(modsum[["upper"]][1:42])),
col = adjustcolor("black", alpha.f = 0.10), border = FALSE)
polygon(c(modsum[["time"]][44:76],rev(modsum[["time"]][44:76])),c(modsum[["lower"]][44:76],rev(modsum[["upper"]][44:76])),
col = adjustcolor("red", alpha.f = 0.10), border = FALSE)
lines(modsum[["time"]][1:43],modsum[["surv"]][1:43],type="s")
lines(modsum[["time"]][44:77],modsum[["surv"]][44:77],type="s",col="red")
legend("topright",legend=c("Control","5mg/ml"),pch=c(15),col=c("black","red"),cex=1)
#dev.off()