#Males Harm Females Less when Competing with Familiar Relatives #Supplementary Material: R Code #Proceedings of the Royal Society of London B: Biological Sciences #Samuel J Lymbery and Leigh W Simmons #Load Libraries library(lme4) library(glmmML) library(pscl) library(AER) library(MASS) library(RVAideMemoire) #Load data data1=read.csv("Lymbery & Simmons 2017 Raw Data.csv") data1$M1P=as.factor(data1$M1P) data1$M2P=as.factor(data1$M2P) data1$M3P=as.factor(data1$M3P) data1$M1P2=as.factor(data1$M1P2) data1$M1P3=as.factor(data1$M3P) data1$ID=as.factor(data1$ID) data1$FP=as.factor(data1$FP) #Subset to remove females laying less than 15 eggs #This includes 1RF, 2RUF and 2URF data2=subset(data1,Eggs2>=15) head(data2) #Subset by treatment R=subset(data2,Relatedness=="R") RF=subset(R,Familiarity=="F") RUF=subset(R,Familiarity=="UF") UR=subset(data2,Relatedness=="UR") URF=subset(UR,Familiarity=="F") URUF=subset(UR,Familiarity=="UF") length(RF$Female) length(RUF$Female) length(URF$Female) length(URUF$Female) #Outlier detection using boxplots (more than 1.5 times interquartile range) #Eggs boxplot(Eggs2~Treatment, names=c("Familiar","Unfamiliar","Familiar","Unfamiliar"), data=data2, ylab=c("Total Number of Eggs Produced per Female"), xlab=c("Related Unrelated"), ylim=c(0,100)) #Offspring boxplot(Offspring~Treatment, names=c("Familiar","Unfamiliar","Familiar","Unfamiliar"), data=data2, ylab=c("Total Number of Offspring Produced per Female"), xlab=c("Related Unrelated"), ylim=c(0,100)) #From the boxplots, groups with ID 72(RF), 80(RF), 171(RUF), 100(URF), 200(URUF), #204(URUF) are outliers for both #eggs and offspring newdata=read.csv("Lymbery & Simmons Raw Data Minus Outliers.csv") newdata$M1P=as.factor(newdata$M1P) newdata$M2P=as.factor(newdata$M2P) newdata$M3P=as.factor(newdata$M3P) newdata$M1P2=as.factor(newdata$M1P2) newdata$M1P3=as.factor(newdata$M3P) newdata$ID=as.factor(newdata$ID) newdata$FP=as.factor(newdata$FP) #So for egg and offspring analysis we use newdata #Subset this by treatment R2=subset(newdata,Relatedness=="R") RF2=subset(R2,Familiarity=="F") RUF2=subset(R2,Familiarity=="UF") UR2=subset(newdata,Relatedness=="UR") URF2=subset(UR2,Familiarity=="F") URUF2=subset(UR2,Familiarity=="UF") length(RF2$Female) length(RUF2$Female) length(URF2$Female) length(URUF2$Female) ###LIFESPAN### #For this we use data2 rather than newdata #Boxplot par(bty='l') boxplot(Lifespan.2~Treatment, names=c("Familiar","Unfamiliar","Familiar","Unfamiliar"), data=data2, ylab=c("Female Lifespan Post-Grouping (Days)"), xlab=c("Related Unrelated"), ylim=c(0,10)) #Interaction Plot meansEb=by(data2$Lifespan.2,list(data2$Familiarity2,data2$Relatedness2),mean) seb=function(x) sqrt(var(x)/length(x)) sesEb=by(data2$Lifespan.2,list(data2$Familiarity2,data2$Relatedness2),seb) pdf(file="Figure 1.pdf", width = 3, height = 3, family = "Helvetica") par(bty='l') par(family="serif",font=1) LSIP=interaction.plot(x.factor=data2$Relatedness2, trace.factor=data2$Familiarity2, response=data2$Lifespan.2, fun=mean, type="b", col=c("black","black"), pch = c(1,19), xlab=NA, ylab=c("female lifespan post-grouping (days)"), ylim=c(4.5,6), legend=FALSE, family="serif", font=1, font.lab=1, cex.axis=0.75, cex.lab=0.75, cex.main=0.75, cex=0.75, las=1) lines(c(1,1),c(meansEb[1]-sesEb[1],meansEb[1]+sesEb[1])) lines(c(1,1),c(meansEb[2]-sesEb[2],meansEb[2]+sesEb[2])) lines(c(2,2),c(meansEb[3]-sesEb[3],meansEb[3]+sesEb[3])) lines(c(2,2),c(meansEb[4]-sesEb[4],meansEb[4]+sesEb[4])) dev.off() #Model without blocking LS1=glm(Lifespan.2~Relatedness*Familiarity,data=data2,family=poisson(link="log")) summary(LS1) 40.331/200 #Ratio=0.20 #It is not overdispersed #Add FP LS2=glmer(Lifespan.2~Relatedness*Familiarity+(1|FP),data=data2,family=poisson(link="log")) overdisp.glmer(LS2) anova(LS2,LS1,test="Chisq") #FP does not explain sig variation so drop it #Now test sig of fixed effects LS3=glm(Lifespan.2~Relatedness+Familiarity,data=data2,family=poisson(link="log")) anova(LS1,LS3,test="Chisq") #p=0.46 #Interaction is not significant so drop it LS4=glm(Lifespan.2~Relatedness,data=data2,family=poisson(link="log")) anova(LS3,LS4,test="Chisq") #p=0.86 #Familiarity is not significant LS5=glm(Lifespan.2~Familiarity,data=data2,family=poisson(link="log")) anova(LS3,LS5,test="Chisq") #p=0.39 #Relatednesss is not significant #Nothing is significant, so we don't need to block ###EGGS### #Boxplot par(bty='l') boxplot(Eggs2~Treatment, names=c("Familiar","Unfamiliar","Familiar","Unfamiliar"), data=newdata, ylab=c("Total Number of Eggs Produced per Female"), xlab=c("Related Unrelated"), ylim=c(0,100)) #Interaction Plot meansEb=by(newdata$Eggs2,list(newdata$Familiarity2,newdata$Relatedness2),mean) seb=function(x) sqrt(var(x)/length(x)) sesEb=by(newdata$Eggs2,list(newdata$Familiarity2,newdata$Relatedness2),seb) pdf(file="Figure 2b.pdf", width = 3, height = 3, family = "Helvetica") par(bty='l') par(family="serif",font=1) EggsIP=interaction.plot(x.factor=newdata$Relatedness2, trace.factor=newdata$Familiarity2, response=newdata$Eggs2, fun=mean, type="b", col=c("black","black"), pch = c(1,19), xlab=NA, ylab=c("number of eggs per female"), ylim=c(60,75), legend=FALSE, family="serif", font=1, font.lab=1, cex.axis=0.75, cex.lab=0.75, cex.main=0.75, cex=0.75, las=1) lines(c(1,1),c(meansEb[1]-sesEb[1],meansEb[1]+sesEb[1])) lines(c(1,1),c(meansEb[2]-sesEb[2],meansEb[2]+sesEb[2])) lines(c(2,2),c(meansEb[3]-sesEb[3],meansEb[3]+sesEb[3])) lines(c(2,2),c(meansEb[4]-sesEb[4],meansEb[4]+sesEb[4])) dev.off() #Models without MP Eggs1=glm(Eggs2~Relatedness*Familiarity,newdata,family=poisson(link="log")) summary(Eggs1) 297.23/194 #Ratio=1.53, so it's not overdispersed #Add FP and test for sig Eggs2=glmer(Eggs2~Relatedness*Familiarity+(1|FP),newdata,family=poisson(link="log")) overdisp.glmer(Eggs2) #Ratio=1.16, so it's not overdispersed anova(Eggs2,Eggs1,test="Chisq") #chisq=16.90,p<0.001 #Keep FP #Test for sig of interaction Eggs3=glmer(Eggs2~Relatedness+Familiarity+(1|FP),data=newdata,family=poisson(link="log")) anova(Eggs2,Eggs3,test="Chisq") #p=0.028 #Interaction is significant #Try Blocking Method One, with full dataset and MP Eggs4=glmer(Eggs2~Relatedness*Familiarity+(1|FP)+(1|M1P2),data=newdata, family=poisson(link="log")) overdisp.glmer(Eggs4) #Ratio is 1.04 so it's not overdispersed #Test sig of MP with FP Eggs5=glmer(Eggs2~Relatedness*Familiarity+(1|FP),data=newdata, family=poisson(link="log")) anova(Eggs4,Eggs5,test="Chisq") #Chisq=1.58, p=0.21 #MP does not explain sig variance when we have FP #Now test sig of FP with MP Eggs6=glmer(Eggs2~Relatedness*Familiarity+(1|M1P2),data=newdata, family=poisson(link="log")) anova(Eggs4,Eggs6,test="Chisq") #Chisq=13.60, p<0.001 #So FP does explain significant variance even when we have MP #The model with FP but not MP may be most parsimonious, but for now we will test the model with both #as well i.e. Eggs4 overdisp.glmer(Eggs4) #Ratio=1.04, not overdispersed #Test sig of interaction Eggs7=glmer(Eggs2~Relatedness+Familiarity+(1|FP)+(1|M1P2),data=newdata, family=poisson(link="log")) anova(Eggs4,Eggs7) #Chisq=3.43, p=0.064 #Interaction is marginally non-sig #Drop it to test for main effects anyway but we may be better off keeping it overdisp.glmer(Eggs7) #Test MP Eggs7b=glmer(Eggs2~Relatedness+Familiarity+(1|FP),data=newdata, family=poisson(link="log")) anova(Eggs7,Eggs7b) #And FP Eggs7c=glmer(Eggs2~Relatedness+Familiarity+(1|M1P2),data=newdata, family=poisson(link="log")) anova(Eggs7,Eggs7c) #Test Familiarity Eggs8=glmer(Eggs2~Relatedness+(1|FP)+(1|M1P2),data=newdata, family=poisson(link="log")) anova(Eggs7,Eggs8) #Chisq=1.22, p=0.27 #Familiarity is not sig Eggs9=glmer(Eggs2~Familiarity+(1|FP)+(1|M1P2),data=newdata, family=poisson(link="log")) anova(Eggs7,Eggs9) #Chisq=0.15, p=0.70 #Relatedness is not sig #Blocking Method Two: Subset to RF and RUF EggsR=glm(Eggs2~Familiarity,data=R2,family=poisson(link="log")) summary(EggsR) 142.2/97 #1.47 so it's not overdispersed #Add FP and test sig EggsR2=glmer(Eggs2~Familiarity+(1|FP),data=R2,family=poisson(link="log")) overdisp.glmer(EggsR2) #Ratio=0.97, not overdispersed anova(EggsR2,EggsR) #Chisq=7.48, p=0.0062 #Keep FP #Sig of familiarity EggsR3=glmer(Eggs2~1+(1|FP),data=R2,family=poisson(link="log")) anova(EggsR2,EggsR3,test="Chisq") #Chisq=6.30, p=0.012 #This says it's significant so continue with blocking EggsR4=glmer(Eggs2~Familiarity+(1|FP)+(1|M1P2),data=R2,family=poisson(link="log")) overdisp.glmer(EggsR4) #Ratio 0.95, so it's not overdispersed anymore #Test for significance of MP with FP anova(EggsR4,EggsR2) #Chisq=0.059, p=081 #MP dpesn't explain sig variance once we have FP #Now try testing for sig of FP with MP EggsR5=glmer(Eggs2~Familiarity+(1|M1P2),data=R2,family=poisson(link="log")) anova(EggsR4,EggsR5) #Chisq=6.83, p=0.0089 #FP still sig #Our most parsimonious model has FP but not MP #We will also test the sig of a model with both (EggsR4), just to be sure #Test for effect of familiarity EggsR6=glmer(Eggs2~(1|FP)+(1|M1P2),data=R2,family=poisson(link="log")) anova(EggsR4,EggsR6) #Chisq=5.62, p=0.018 #Familiarity is significant #Blocking Method Two continued: Subset to URF and URUF EggsUR1=glm(Eggs2~Familiarity,data=UR2,family=poisson(link="log")) summary(EggsUR1) 155.03/97 #Ratio=1.60, so it's not overdispersed #Add FP and test for sig EggsUR2=glmer(Eggs2~Familiarity+(1|FP),data=UR2,family=poisson(link="log")) overdisp.glmer(EggsUR2) #Ratio=1.09, not overdispersed anova(EggsUR2,EggsUR1,test="Chisq") #Chisq=10.18, p=0.0014 #Keep FP #Test for sig of familiarity EggsUR3=glmer(Eggs2~(1|FP),data=UR2,family=poisson(link="log")) anova(EggsUR3,EggsUR2) #Chisq=0.27, p=0.61 #Not sig so no need to continue with blocking ###OFFSPRING### #Boxplot par(bty='l') boxplot(Offspring~Treatment, names=c("Familiar","Unfamiliar","Familiar","Unfamiliar"), data=newdata, ylab=c("Total Number of Offspring Produced per Female"), xlab=c("Related Unrelated"), ylim=c(0,100)) #Interaction Plot meansEb=by(newdata$Offspring,list(newdata$Familiarity2,newdata$Relatedness2),mean) seb=function(x) sqrt(var(x)/length(x)) sesEb=by(newdata$Offspring,list(newdata$Familiarity2,newdata$Relatedness2),seb) pdf(file="Figure 3b.pdf", width = 3, height = 3, family = "Helvetica") par(bty='l') par(family="serif",font=1) OffIP=interaction.plot(x.factor=newdata$Relatedness2, trace.factor=newdata$Familiarity2, response=newdata$Offspring, fun=mean, type="b", col=c("black","black"), pch = c(1,19), xlab=NA, ylab=c("number of offspring per female"), ylim=c(50,65), legend=FALSE, family="serif", font=1, font.lab=1, cex.axis=0.75, cex.lab=0.75, cex.main=0.75, cex=0.75, las=1) lines(c(1,1),c(meansEb[1]-sesEb[1],meansEb[1]+sesEb[1])) lines(c(1,1),c(meansEb[2]-sesEb[2],meansEb[2]+sesEb[2])) lines(c(2,2),c(meansEb[3]-sesEb[3],meansEb[3]+sesEb[3])) lines(c(2,2),c(meansEb[4]-sesEb[4],meansEb[4]+sesEb[4])) dev.off() #Models without MP Off1=glm(Offspring~Relatedness*Familiarity,newdata,family=poisson(link="log")) summary(Off1) 362.47/194 #Ratio=1.87, so it's not overdispersed #Add FP and test for sig Off2=glmer(Offspring~Relatedness*Familiarity+(1|FP),newdata,family=poisson(link="log")) overdisp.glmer(Off2) #Ratio=1.32, so it's not overdispersed anova(Off2,Off1,test="Chisq") #chisq=34.34,p<0.001 #Keep FP #Test for sig of interaction Off3=glmer(Offspring~Relatedness+Familiarity+(1|FP),data=newdata,family=poisson(link="log")) anova(Off2,Off3,test="Chisq") #Chisq=6.46,p=0.011 #Interaction is significant #Try Blocking Method One, with full dataset and MP Off4=glmer(Offspring~Relatedness*Familiarity+(1|FP)+(1|M1P2),data=newdata, family=poisson(link="log")) overdisp.glmer(Off4) #Ratio is 1.11 so it's not overdispersed #Test sig of MP with FP Off5=glmer(Offspring~Relatedness*Familiarity+(1|FP),data=newdata, family=poisson(link="log")) anova(Off4,Off5,test="Chisq") #Chisq=3.36, p=0.067 #MP is marginally non sig when we have FP #Now test sig of FP with MP Off6=glmer(Offspring~Relatedness*Familiarity+(1|M1P2),data=newdata, family=poisson(link="log")) anova(Off4,Off6,test="Chisq") #Chisq=26.45, p<0.001 #So FP does explain significant variance even when we have MP #The model with FP but not MP may be most parsimonious, but MP is only marginally non-sig #so for now we will test the model with both as well i.e. Off4 overdisp.glmer(Off4) #Ratio=1.11, not overdispersed #Test sig of interaction Off7=glmer(Offspring~Relatedness+Familiarity+(1|FP)+(1|M1P2),data=newdata, family=poisson(link="log")) anova(Off4,Off7) #Chisq=4.10, p=0.043 #Interaction is sig so we don't test main effects #Blocking Method Two: Subset to RF and RUF OffR=glm(Offspring~Familiarity,data=R2,family=poisson(link="log")) summary(OffR) 199.9/97 #2.06 so it's overdispersed #Add FP and test sig OffR2=glmer(Offspring~Familiarity+(1|FP),data=R2,family=poisson(link="log")) overdisp.glmer(OffR2) #Ratio=1.23, not overdispersed anymore anova(OffR2,OffR) #Chisq=20.34, p<0.001 #Keep FP #Sig of familiarity OffR3=glmer(Offspring~(1|FP),data=R2,family=poisson(link="log")) anova(OffR2,OffR3,test="Chisq") #Chisq=6.48, p=0.011 #This says it's significant so continue with blocking OffR4=glmer(Offspring~Familiarity+(1|FP)+(1|M1P2),data=R2,family=poisson(link="log")) overdisp.glmer(OffR4) #Ratio 1.002, so it's not overdispersed #Test for significance of MP with FP anova(OffR4,OffR2) #Chisq=3.40, p=0.065 #MP is marginally non sig when we have FP #Now try testing for sig of FP with MP OffR5=glmer(Offspring~Familiarity+(1|M1P2),data=R2,family=poisson(link="log")) anova(OffR4,OffR5) #Chisq=15.70, p<0.001 #FP still sig #Our most parsimonious model may have FP but not MP, but MP is only marginally non-sig #so we will also test the sig of a model with both (OffR4) #Test for effect of familiarity OffR6=glmer(Offspring~(1|FP)+(1|M1P2),data=R2,family=poisson(link="log")) anova(OffR4,OffR6) #Chisq=5.18, p=0.023 #Familiarity is still significant #Blocking Method Two continued: Subset to URF and URUF OffUR1=glm(Offspring~Familiarity,data=UR2,family=poisson(link="log")) summary(OffUR1) 162.56/97 #Ratio=1.68, so it's not overdispersed #Add FP and test for sig OffUR2=glmer(Offspring~Familiarity+(1|FP),data=UR2,family=poisson(link="log")) overdisp.glmer(OffUR2) #Ratio=1.12, not overdispersed anova(OffUR2,OffUR1,test="Chisq") #Chisq=14.14, p=0.00015 #Keep FP #Test for sig of familiarity OffUR3=glmer(Offspring~(1|FP),data=UR2,family=poisson(link="log")) anova(OffUR3,OffUR2) #Chisq=0.67, p=0.41 #Not sig so no need to continue with blocking ###DISTRIBUTION OF FP### library(aylmer) str(newdata$FP) summary(RF2$FP) summary(RUF2$FP) summary(URF2$FP) summary(URUF2$FP) FPTable=matrix(c(summary(RF2$FP),summary(RUF2$FP),summary(URF2$FP),summary(URUF2$FP)),byrow=TRUE,nrow=64) is.na(FPTable)=!FPTable FPTable View(FPTable) aylmer.test(FPTable,simulate.p.value=TRUE,B=100000) #P~1, so they are randomly distributed across treatments #We do not need to include FP