#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