################################################################### # # Code to produce Figure 2 in: Burgess, Sander, Bueno (xxxx) # How relatedness between mates influences reproductive success: # an experimental analysis of self-fertilization and biparental # inbreeding in a marine bryozoan. Ecology and Evolution #: ##-## # #################################################################### # setwd("") library(lme4) # version 1.1-21 dat <- read.csv("Data for figure 2 and 3.csv",header=T,colClasses=c(rep("factor",6),rep("numeric",4))) ## Make sure to check 30.10 imports correctly (not as 30.1) # Calculate juveniles per ovicell for plotting dat$Prop.Sett <- dat$Settlers / dat$Ovicells dat$Prop.Sett <- ifelse(is.nan(dat$Prop.Sett),0,dat$Prop.Sett) # Change order of treatment levels for plotting dat$Treatment <- factor(dat$Treatment, levels=c("Within","Sib","Alone")) # Panel A (Site 1, Dog Island) # Fit models d <- dat[dat$Site=="Dog Island" & dat$Treatment!="Alone",] m1 <- glmer(Settlers ~ Treatment + (1|GrandMother.ID),data=d,family="poisson") m2 <- glmer(Settlers ~ 1 + (1|GrandMother.ID),data=d,family="poisson") # Get p-values pvalA <- anova(m1,m2,test="Chisq")$'Pr(>Chisq)'[2] pvalA <- ifelse(pvalA>0.001,round(pvalA,3),"<0.001") # Get fitted values m1 <- glmer(Settlers ~ Treatment-1 + (1|GrandMother.ID),data=d,family="poisson") newdat <- with(d, expand.grid(Treatment=c("Within","Sib"),Settlers=0)) newdat$pred <- predict(m1,newdat,re.form=NA) mm <- model.matrix(terms(m1),newdat)[,1:2] # Uncertainty in fixed effects newdat$V <- diag(mm %*% tcrossprod(vcov(m1),mm)) newdat$plo <- newdat$pred - 1.96*sqrt(newdat$V) newdat$phi <- newdat$pred + 1.96*sqrt(newdat$V) # Uncertainty in fixed effects and random effect variance newdat$Vr <- newdat$V + VarCorr(m1)$GrandMother.ID[1] newdat$plor <- newdat$pred - 1.96*sqrt(newdat$Vr) newdat$phir <- newdat$pred + 1.96*sqrt(newdat$Vr) mnA <- exp(newdat$pred) lwrA <- exp(newdat$plor) uprA <- exp(newdat$phir) # Panel B (Site 2, Marine Lab) d <- dat[dat$Site=="MLW" & dat$Treatment!="Alone",] m1 <- glmer(Settlers ~ Treatment + (1|GrandMother.ID),data=d,family="poisson") m2 <- glmer(Settlers ~ 1 + (1|GrandMother.ID),data=d,family="poisson") # Get p-values pvalB <- anova(m1,m2,test="Chisq")$'Pr(>Chisq)'[2] pvalB <- ifelse(pvalB>0.001,round(pvalB,3),"<0.001") # Get fitted values m1 <- glmer(Settlers ~ Treatment-1 + (1|GrandMother.ID),data=d,family="poisson") newdat <- with(d, expand.grid(Treatment=c("Within","Sib"),Settlers=0)) newdat$pred <- predict(m1,newdat,re.form=NA) mm <- model.matrix(terms(m1),newdat)[,1:2] # Uncertainty in fixed effects newdat$V <- diag(mm %*% tcrossprod(vcov(m1),mm)) newdat$plo <- newdat$pred - 1.96*sqrt(newdat$V) newdat$phi <- newdat$pred + 1.96*sqrt(newdat$V) # Uncertainty in fixed effects and random effect variance newdat$Vr <- newdat$V + VarCorr(m1)$GrandMother.ID[1] newdat$plor <- newdat$pred - 1.96*sqrt(newdat$Vr) newdat$phir <- newdat$pred + 1.96*sqrt(newdat$Vr) mnB <- exp(newdat$pred) lwrB <- exp(newdat$plor) uprB <- exp(newdat$phir) # Panel C (Site 1, Dog Island) d <- dat[dat$Site=="Dog Island" & dat$Treatment!="Alone",] m1 <- glmer(cbind(Settlers,Ovicells) ~ Treatment + (1|GrandMother.ID),data=d,family="binomial") m2 <- glmer(cbind(Settlers,Ovicells) ~ 1 + (1|GrandMother.ID),data=d,family="binomial") # Get p-values pvalC <- anova(m1,m2,test="Chisq")$'Pr(>Chisq)'[2] pvalC <- ifelse(pvalC>0.001,round(pvalC,3),"<0.001") # Get fitted values m1 <- glmer(cbind(Settlers,Ovicells) ~ Treatment-1 + (1|GrandMother.ID),data=d,family="binomial") newdat <- with(d, expand.grid(Treatment=c("Within","Sib"),Settlers=0,Ovicells=0)) newdat$pred <- predict(m1,newdat,re.form=NA) mm <- model.matrix(terms(m1),newdat)[,1:2] # Uncertainty in fixed effects newdat$V <- diag(mm %*% tcrossprod(vcov(m1),mm)) newdat$plo <- newdat$pred - 1.96*sqrt(newdat$V) newdat$phi <- newdat$pred + 1.96*sqrt(newdat$V) # Uncertainty in fixed effects and random effect variance newdat$Vr <- newdat$V + VarCorr(m1)$GrandMother.ID[1] newdat$plor <- newdat$pred - 1.96*sqrt(newdat$Vr) newdat$phir <- newdat$pred + 1.96*sqrt(newdat$Vr) mnC <- plogis(newdat$pred) lwrC <- plogis(newdat$plor) uprC <- plogis(newdat$phir) # Panel D (Site 2, Marine Lab) d <- dat[dat$Site=="MLW" & dat$Treatment!="Alone",] m1 <- glmer(cbind(Settlers,Ovicells) ~ Treatment + (1|GrandMother.ID),data=d,family="binomial") m2 <- glmer(cbind(Settlers,Ovicells) ~ 1 + (1|GrandMother.ID),data=d,family="binomial") # Get p-values pvalD <- anova(m1,m2,test="Chisq")$'Pr(>Chisq)'[2] pvalD <- ifelse(pvalD>0.001,round(pvalD,3),"<0.001") # Get fitted values m1 <- glmer(cbind(Settlers,Ovicells) ~ Treatment-1 + (1|GrandMother.ID),data=d,family="binomial") newdat <- with(d, expand.grid(Treatment=c("Within","Sib"),Settlers=0,Ovicells=0)) newdat$pred <- predict(m1,newdat,re.form=NA) mm <- model.matrix(terms(m1),newdat)[,1:2] # Uncertainty in fixed effects newdat$V <- diag(mm %*% tcrossprod(vcov(m1),mm)) newdat$plo <- newdat$pred - 1.96*sqrt(newdat$V) newdat$phi <- newdat$pred + 1.96*sqrt(newdat$V) # Uncertainty in fixed effects and random effect variance newdat$Vr <- newdat$V + VarCorr(m1)$GrandMother.ID[1] newdat$plor <- newdat$pred - 1.96*sqrt(newdat$Vr) newdat$phir <- newdat$pred + 1.96*sqrt(newdat$Vr) mnD <- plogis(newdat$pred) lwrD <- plogis(newdat$plor) uprD <- plogis(newdat$phir) # Set plotting parameters col1 <- adjustcolor("#D55E00",alpha.f=0.5) # orange (within) col2 <- adjustcolor("#0072B2",alpha.f=0.5) # blue (sib) col3 <- adjustcolor("#999999",alpha.f=0.5) # grey (alone) pchs <- 0:3 cexs <- 2 cex.ax <- 1.5 xlabs <- c("Non-sib\nmating","Sib\nmating","Selfing") # Plot Figure 2 quartz(width=7,height=5) par(mfrow=c(2,2),mar=c(0,1,1,1),oma=c(5.2,4.5,1,0.5)) # Pabel A d <- dat[dat$Site=="Dog Island",] ind <- unique(d$GrandMother.ID) plot(c(0.7,3.2),c(0,250),type="n",las=1,xaxt="n",yaxt="n",ylab="",xlab="",bty="l") offset <- -0.1 for(i in 1:length(ind)){ with(d[d$GrandMother.ID==ind[i],], points(as.numeric(Treatment)+offset,Settlers, pch=pchs[i],cex=cexs,lwd=2, col=ifelse(Treatment=="Within",col1,ifelse(Treatment=="Sib",col2,col3)))) offset <- offset + 0.05 } axis(side=2,at=seq(0,250,50),labels=seq(0,250,50),cex.axis=cex.ax,las=1) mtext("Site 1 (Dog Island)",side=3,adj=0.5,cex=1.2) mtext("Number of juveniles",side=2,line=4,cex=1.2) legend(0.3,250,legend=eval(paste("A) p =",pvalA)),cex=cex.ax,bty="n") points(1:3,c(mnA,0),cex=3,pch=19,col=adjustcolor("black",alpha.f=0.5)) segments(1:2,lwrA,1:2,uprA,lwd=2) # Panel B d <- dat[dat$Site=="MLW",] ind <- unique(d$GrandMother.ID) plot(c(0.7,3.2),c(0,250),type="n",las=1,xaxt="n",yaxt="n",ylab="",xlab="",bty="l") offset <- -0.1 for(i in 1:length(ind)){ with(d[d$GrandMother.ID==ind[i],], points(as.numeric(Treatment)+offset,Settlers, pch=pchs[i],cex=cexs,lwd=2, col=ifelse(Treatment=="Within",col1,ifelse(Treatment=="Sib",col2,col3)))) offset <- offset + 0.05 } axis(side=2,at=seq(0,250,50),labels=seq(0,250,50),cex.axis=cex.ax,las=1) mtext("Site 2 (Marine Lab)",side=3,adj=0.5,cex=1.2) legend(0.3,250,legend=eval(paste("B) p =",pvalB)),cex=cex.ax,bty="n") points(1:3,c(mnB,0),cex=3,pch=19,col=adjustcolor("black",alpha.f=0.5)) segments(1:2,lwrB,1:2,uprB,lwd=2) # Panel C d <- dat[dat$Site=="Dog Island",] ind <- unique(d$GrandMother.ID) plot(c(0.7,3.2),c(0,0.65),type="n",las=1,xaxt="n",yaxt="n",ylab="",xlab="",bty="l") offset <- -0.1 for(i in 1:length(ind)){ with(d[d$GrandMother.ID==ind[i],], points(as.numeric(Treatment)+offset,Prop.Sett, pch=pchs[i],cex=cexs,lwd=2, col=ifelse(Treatment=="Within",col1,ifelse(Treatment=="Sib",col2,col3)))) offset <- offset + 0.05 } axis(side=2,at=seq(0,0.6,0.1),labels=seq(0,0.6,0.1),cex.axis=cex.ax,las=1) mtext("Juveniles per ovicell",side=2,line=4,cex=1.2) legend(0.3,0.65,legend=eval(paste("C) p =",pvalC)),cex=cex.ax,bty="n") mtext(xlabs,side=1,line=2,at=c(1,2,2.9)) points(1:3,c(mnC,0),cex=3,pch=19,col=adjustcolor("black",alpha.f=0.5)) segments(1:2,lwrC,1:2,uprC,lwd=2) # D d <- dat[dat$Site=="MLW",] ind <- unique(d$GrandMother.ID) plot(c(0.7,3.2),c(0,0.65),type="n",las=1,xaxt="n",yaxt="n",ylab="",xlab="",bty="l") offset <- -0.1 for(i in 1:length(ind)){ with(d[d$GrandMother.ID==ind[i],], points(as.numeric(Treatment)+offset,Prop.Sett, pch=pchs[i],cex=cexs,lwd=2, col=ifelse(Treatment=="Within",col1,ifelse(Treatment=="Sib",col2,col3)))) offset <- offset + 0.05 } axis(side=2,at=seq(0,0.6,0.1),labels=seq(0,0.6,0.1),cex.axis=cex.ax,las=1) legend(0.3,0.65,legend=eval(paste("D) p =",pvalD)),cex=cex.ax,bty="n") mtext(xlabs,side=1,line=2,at=c(1,2,2.9)) points(1:3,c(mnD,0),cex=3,pch=19,col=adjustcolor("black",alpha.f=0.5)) segments(1:2,lwrD,1:2,uprD,lwd=2) mtext("Mating partner treatment",side=1,outer=T,line=4,cex=1.2)