Figure 1D - Leslie fly fitness matrix

library(popbio)

# FITNESS MATRIX WITH VARIABLE DEVELOPMENT TIME PER VIAL (100x) 
fitdat<-read.csv(FlyFitnessRawData.csv↩, header=TRUE, sep=",") 
devdat<-read.csv(DevelopmentData.csv↩, header=TRUE, sep=",") 

uniq <- unique(unlist(fitdat$treat))
uniq1 <- unique(unlist(fitdat$vial))
fit_rep<-vector("numeric", 100)
fit_vial<-vector("numeric", 5)
d<-matrix(0, nrow=5, length(uniq), byrow=F)
d<-as.data.frame(d)

# iterate through treatments 
for(m in 1:length(uniq)){

trt<-subset(fitdat, treat == uniq[m])

# iterate through vials 
for(n in 1:length(uniq1)){

	vial<-subset(trt, vial == uniq1[n])
	fecund<-vector("numeric", 39)
	surv<-vector("numeric", 39)

# select development time within treatment to use 	
	for (i in 1:100){	
	
		trt2<-subset(devdat, treatment == uniq1[n])
		uniq2<-unique(unlist(trt2$replicate))
		r<-sample(uniq2, 1, replace=TRUE)
		rep<-subset(trt2, replicate==r)
		dev<-rep$development

		fecund<-vial$fecundity_adj
		f_add<-rep(0,dev+6)
		fecundity<-c(f_add, fecund)
		time<-seq(0, length(fecundity)-1, by=1)
	
		surv<-vial$survival
		s_add<-rep(1,dev+5)
		survival<-c(s_add, surv)

#create matrix of all zeros
		fitM<-matrix(0, nrow=length(fecundity), length(fecundity), byrow=T)

#populate matirx with fecundity and survival values
		for (k in 1:length(fecundity)){
				
		fitM[1,k]<-fecundity[k]

		for (l in 1:length(survival)){
			
			fitM[l+1,l] <- survival[l]
				
			}
		}		
		
#calculate dominant eigenvalue
	fit<-lambda(fitM)

#store results in vector
	fit_rep[i]<-fit
	
	}
	
fit_vial[n]<-mean(fit_rep)

}

#store results in dataframe	
d[,m]<-fit_vial

}

d

check <- read.csv(FlyFitness.csv↩)

# PLOT FLY FITNESS (from fitness dataframe above (d) modified in Excel)

mydata<-read.csv(FlyFitness.csv↩, header=TRUE, sep=",")
data<-subset(mydata, treat != "33")
attach(data)

st.err <- function(x) {
		sd(x)/sqrt(length(x))
		}

stderr<-aggregate(fitness, by=list(treat), FUN="st.err")
meanFit<-aggregate(fitness, by=list(treat), FUN="mean")


spacing<-c(0.92, 1.08, 0.85, 1, 1.15, 1.62, 1.59, 2.28, 1.75, 2.19, 2.05, 1.95, 2.11, 1.71, 1.87, 2.78, 2.94, 3.02, 3.1, 3.2, 3.28, 2.66, 3.375, 2.56, 3.55, 3.9, 3.99, 4.08, 4.17, 4.33, 5, 0)
points<-c(1:32)
points<-as.character(points)

par(mar=c(5,5,2,2))
plot(spacing, meanFit$x, ylim=c(1.15, 1.25), ylab="Fitness", xlab="Number of species", pch=21, cex=3.5, cex.axis=1.4, cex.lab=1.6)
text(spacing, meanFit$x, labels=points)
arrows(spacing,meanFit$x-stderr$x,spacing,meanFit$x+stderr$x, code=3, length=0.02, angle = 90)