# You will need to alter the setwd to the directory where you have the needed data files. setwd("/Users/mmgdepartment/Dropbox/All Hands") # First we read in the raw data data.raw<-read.table("All.Hands.Concatenated.Data.csv", sep=",", header=T) # Then we make sure column types are correct. For example, Person should be a # categorical variable, not a continuous one. data.raw$Flask<-factor(data.raw$Flask) data.raw$Person<-factor(data.raw$Person) # Next we clean up the raw data. Any value of 0 will be converted to a missing value; # any row with missing values will be dropped, since we cannot calculated fitness with # missing or 0 values. data.trimmed<-data.raw data.trimmed[data.trimmed==0]<- NA data.trimmed<-na.omit(data.trimmed) # Next we need to calculate fitness. This is based on dividing the number of generations # in the evolved population by the number of generations in the reference clone. So, # first we calculate the number of generations of each population in the competition. data.trimmed$Gen.R<-log((data.trimmed$R.3*2*10^9)/(data.trimmed$R.0*2*10^3), base=2) data.trimmed$Gen.W<-log((data.trimmed$W.3*2*10^9)/(data.trimmed$W.0*2*10^3), base=2) # Then we subset the data differently depending on whether the reference population is # arabinose positive or negative. data.positive<-subset(data.trimmed, Arabinose=="Positive") data.negative<-subset(data.trimmed, Arabinose=="Negative") # Then we calculate fitness data.positive$Fitness<-data.positive$Gen.W/data.positive$Gen.R data.negative$Fitness<-data.negative$Gen.R/data.negative$Gen.W # We can also take this opportunity to create a new column, of what the evolved population # is. data.positive$Evolved.Population<-data.positive$White.Pop data.negative$Evolved.Population<-data.negative$Red.Pop # And then we combine the two data sets back together. data.combined<-rbind(data.positive, data.negative) summary(data.combined) # Now we log transform the fitness data data.combined$Log.Fitness<-log(data.combined$Fitness, base=10) #data.combined$Log.Fitness<-log(data.combined$Fitness, base=2.718281828459045) # Now let's sort the data file in a way to make it easily readable data.combined<-data.combined[order(data.combined$Replicate, data.combined$Person, data.combined$Flask), ] # Now we calculate the mean and standard deviation of each population-time point combination: data.transformed<-c() flask.list<-unique(data.combined$Flask) for (i in 1:length(flask.list) ) { current.flask<-flask.list[i] current.data<-subset(data.combined, Flask==current.flask) current.mean<-mean(current.data$Log.Fitness) current.sd<-sd(current.data$Log.Fitness) current.data$Flask.Mean<-current.mean current.data$Flask.sd<-current.sd data.transformed<-rbind(data.transformed, current.data) } summary(data.transformed) data.transformed<-data.transformed[order(data.transformed$Replicate, data.transformed$Person), ] # Now we calculate the Z score for each measured fitness: data.transformed$Flask.z<-(data.transformed$Flask.Mean-data.transformed$Log.Fitness)/data.transformed$Flask.sd # Now we want to know how many measurements fall within certain distances from their # relevant means. To do this, we don't care whether the measurement is above or below # its mean, merely the magnitude of its distance. So we take the absolute value of # the z score, and then we'll bin those results. data.transformed$Absolute.z<-abs(data.transformed$Flask.z) data.4.or.more<-subset(data.transformed, Absolute.z>=4) data.3.to.4<-subset(data.transformed, Absolute.z>=3 & Absolute.z<4) data.2.to.3<-subset(data.transformed, Absolute.z>=2 & Absolute.z<3) data.1.to.2<-subset(data.transformed, Absolute.z>=1 & Absolute.z<2) data.less.than.1<-subset(data.transformed, Absolute.z<1) length(data.4.or.more$Flask) length(data.3.to.4$Flask) length(data.2.to.3$Flask) # We will accept anything with an absolute z score of less than 3. data.accepted<-data.transformed[data.transformed$Absolute.z<3, ] # This data file now has 1113 entries, out of the 1126 measurements which produced # calculable fitness measurements ############################################### # Now we will create Table 1, showing estimates of mean fitness for each of these # populations, at each of the time points. Population.list<-unique(data.accepted$Evolved.Population) data.40k<-subset(data.accepted, Generation==40000) data.50k<-subset(data.accepted, Generation==50000) data.60k<-subset(data.accepted, Generation==60000) Mean.Fit.40k<-matrix(nrow=length(Population.list), ncol=4) Mean.Fit.50k<-matrix(nrow=length(Population.list), ncol=4) Mean.Fit.60k<-matrix(nrow=length(Population.list), ncol=4) for (i in 1:length(Population.list) ) { current.pop<-Population.list[i] current.data<-subset(data.40k, Evolved.Population==current.pop) N<-length(current.data$Flask) Mean.Fit<-mean(current.data$Fitness) SD.Fit<-sd(current.data$Fitness) Mean.Fit.40k[i, 1]<-as.character(current.pop) Mean.Fit.40k[i, 2]<-as.numeric(as.character(N)) Mean.Fit.40k[i, 3]<-as.numeric(as.character(Mean.Fit)) Mean.Fit.40k[i, 4]<-as.numeric(as.character(SD.Fit)) } for (i in 1:length(Population.list) ) { current.pop<-Population.list[i] current.data<-subset(data.50k, Evolved.Population==current.pop) N<-length(current.data$Flask) Mean.Fit<-mean(current.data$Fitness) SD.Fit<-sd(current.data$Fitness) Mean.Fit.50k[i, 1]<-as.character(current.pop) Mean.Fit.50k[i, 2]<-as.numeric(as.character(N)) Mean.Fit.50k[i, 3]<-as.numeric(as.character(Mean.Fit)) Mean.Fit.50k[i, 4]<-as.numeric(as.character(SD.Fit)) } for (i in 1:length(Population.list) ) { current.pop<-Population.list[i] current.data<-subset(data.60k, Evolved.Population==current.pop) N<-length(current.data$Flask) Mean.Fit<-mean(current.data$Fitness) SD.Fit<-sd(current.data$Fitness) Mean.Fit.60k[i, 1]<-as.character(current.pop) Mean.Fit.60k[i, 2]<-as.numeric(as.character(N)) Mean.Fit.60k[i, 3]<-as.numeric(as.character(Mean.Fit)) Mean.Fit.60k[i, 4]<-as.numeric(as.character(SD.Fit)) } colnames(Mean.Fit.40k)<-c("Population", "N", "Mean", "SD") colnames(Mean.Fit.50k)<-c("Population", "N", "Mean", "SD") colnames(Mean.Fit.60k)<-c("Population", "N", "Mean", "SD") Grand.Mean.40k<-matrix(nrow=1, ncol=4) Grand.Mean.50k<-matrix(nrow=1, ncol=4) Grand.Mean.60k<-matrix(nrow=1, ncol=4) Grand.Mean.40k[1, 1]<-as.character("Grand Mean") Grand.Mean.40k[1, 2]<-as.numeric(length(Mean.Fit.40k[,3])) Grand.Mean.40k[1, 3]<-as.numeric(mean(as.numeric(Mean.Fit.40k[,3]))) Grand.Mean.40k[1, 4]<-as.numeric(sd(as.numeric(Mean.Fit.40k[,3]))) Grand.Mean.50k[1, 1]<-as.character("Grand Mean") Grand.Mean.50k[1, 2]<-as.numeric(length(Mean.Fit.50k[,3])) Grand.Mean.50k[1, 3]<-as.numeric(mean(as.numeric(Mean.Fit.50k[,3]))) Grand.Mean.50k[1, 4]<-as.numeric(sd(as.numeric(Mean.Fit.50k[,3]))) Grand.Mean.60k[1, 1]<-as.character("Grand Mean") Grand.Mean.60k[1, 2]<-as.numeric(length(Mean.Fit.60k[,3])) Grand.Mean.60k[1, 3]<-as.numeric(mean(as.numeric(Mean.Fit.60k[,3]))) Grand.Mean.60k[1, 4]<-as.numeric(sd(as.numeric(Mean.Fit.60k[,3]))) colnames(Grand.Mean.40k)<-c("Population", "N", "Mean", "SD") colnames(Grand.Mean.50k)<-c("Population", "N", "Mean", "SD") colnames(Grand.Mean.60k)<-c("Population", "N", "Mean", "SD") Mean.Fit.40k<-rbind(Mean.Fit.40k, Grand.Mean.40k) Mean.Fit.50k<-rbind(Mean.Fit.50k, Grand.Mean.50k) Mean.Fit.60k<-rbind(Mean.Fit.60k, Grand.Mean.60k) # This portion rounds the data for better display Mean.Fit.40k[,3]<-round(as.numeric(Mean.Fit.40k[,3]), digits=4) Mean.Fit.40k[,4]<-round(as.numeric(Mean.Fit.40k[,4]), digits=4) Mean.Fit.50k[,3]<-round(as.numeric(Mean.Fit.50k[,3]), digits=4) Mean.Fit.50k[,4]<-round(as.numeric(Mean.Fit.50k[,4]), digits=4) Mean.Fit.60k[,3]<-round(as.numeric(Mean.Fit.60k[,3]), digits=4) Mean.Fit.60k[,4]<-round(as.numeric(Mean.Fit.60k[,4]), digits=4) Table.1<-cbind(Mean.Fit.40k, Mean.Fit.50k[, 2:4], Mean.Fit.60k[, 2:4] ) Table.1<-as.data.frame(Table.1) # This places the data in rows by population, and in columns of first all information # from generation 40,000; then all information from generation 50,000; then all information # from generation 60,000. I don't know how to create column spanning names in R, so that # can be done later in a spreadsheet application. # Tables were renumbered on submission, as many were moved to the Supplement. What used to be Table 1 became Table S2. write.csv(file="Table.S2.csv", Table.1) ############################################## # Now we will create Table 2. # The easiest way to do this is to create a new data file, which has 3 columns # of Fitness data for each replicate: 1 for 40,000 generations; 1 for 50,000 # generations, and 1 for 60,000 generations. # We will need to create a new column concatenating Person and Replicate, # so that the 22 replicates from Person 1 are treated independent of each other. data.accepted$Person.Plus.Replicate<-paste(as.character(data.accepted$Person), as.character(data.accepted$Replicate)) Population.list<-unique(data.accepted$Evolved.Population) Replicate.list<-unique(data.accepted$Person.Plus.Replicate) data.by.change<-matrix(nrow=9*length(Replicate.list), ncol=8) i<-1 for (i in 1:length(Population.list)) { current.population<-Population.list[i] current.data<-subset(data.accepted, Evolved.Population==current.population) j<-1 for (j in 1:length(Replicate.list)) { current.replicate<-Replicate.list[j] inner.data<-subset(current.data, Person.Plus.Replicate==current.replicate) Fit.40<-subset(inner.data, Generation==40000) Fit.50<-subset(inner.data, Generation==50000) Fit.60<-subset(inner.data, Generation==60000) data.by.change[(j-1)*9+i, 1]<-as.character(current.population) data.by.change[(j-1)*9+i, 2]<-as.character(current.replicate) if (length(Fit.40$Fitness)>0) { data.by.change[(j-1)*9+i, 3]<-as.numeric(as.character(Fit.40$Fitness)) } if (length(Fit.50$Fitness)>0) { data.by.change[(j-1)*9+i, 4]<-as.numeric(as.character(Fit.50$Fitness)) } if (length(Fit.60$Fitness)>0) { data.by.change[(j-1)*9+i, 5]<-as.numeric(as.character(Fit.60$Fitness)) } } } for (i in 1: length(data.by.change[,1])) { if (!is.na(data.by.change[i, 3])==TRUE) { if (!is.na(data.by.change[i, 5])==TRUE) { value<-as.numeric(as.character(data.by.change[i, 5]))-as.numeric(as.character(data.by.change[i, 3])) data.by.change[i, 6]<-value } if (!is.na(data.by.change[i, 4])==TRUE) { value<-as.numeric(as.character(data.by.change[i, 4]))-as.numeric(as.character(data.by.change[i, 3])) data.by.change[i, 7]<-value } } if (!is.na(data.by.change[i, 4])==TRUE) { if (!is.na(data.by.change[i, 5])==TRUE) { value<-as.numeric(as.character(data.by.change[i, 5]))-as.numeric(as.character(data.by.change[i, 4])) data.by.change[i, 8]<-value } } } colnames(data.by.change)<-c("Population", "Person.Plus.Replicate", "Fit.40k", "Fit.50k", "Fit.60k", "Fit.60.Minus.Fit.40", "Fit.50.Minus.Fit.40", "Fit.60.Minus.Fit.50") data.by.change<-as.data.frame(data.by.change) # Now that the data file has been created, we can use it to find how mean fitness has changed in each of the # the populations at each time range we care about. data.for.table.2<-na.omit(data.by.change) Table.2<-matrix(nrow=9, ncol=14) Population.list<-unique(data.for.table.2$Population) for (i in 1:length(Population.list)) { current.pop<-Population.list[i] current.data<-subset(data.for.table.2, Population==current.pop) N<-length(current.data$Fit.40k) dif.1<-as.numeric(as.character(current.data$Fit.60.Minus.Fit.40)) dif.2<-as.numeric(as.character(current.data$Fit.50.Minus.Fit.40)) dif.3<-as.numeric(as.character(current.data$Fit.60.Minus.Fit.50)) Mean.Fit.40.to.60<-mean(dif.1) SD.Fit.40.to.60<-sd(dif.1) Mean.Fit.40.to.50<-mean(dif.2) SD.Fit.40.to.50<-sd(dif.2) Mean.Fit.50.to.60<-mean(dif.3) SD.Fit.50.to.60<-sd(dif.3) values.1<-as.numeric(as.character(current.data$Fit.60.Minus.Fit.40)) values.2<-as.numeric(as.character(current.data$Fit.50.Minus.Fit.40)) values.3<-as.numeric(as.character(current.data$Fit.60.Minus.Fit.50)) t.40.to.60<-as.numeric(as.character(t.test(x=values.1, alternative="greater", mu=0)[1])) p.40.to.60<-as.numeric(as.character(t.test(x=values.1, alternative="greater", mu=0)[3])) t.40.to.50<-as.numeric(as.character(t.test(x=values.2, alternative="greater", mu=0)[1])) p.40.to.50<-as.numeric(as.character(t.test(x=values.2, alternative="greater", mu=0)[3])) t.50.to.60<-as.numeric(as.character(t.test(x=values.3, alternative="greater", mu=0)[1])) p.50.to.60<-as.numeric(as.character(t.test(x=values.3, alternative="greater", mu=0)[3])) Table.2[i, 1]<-as.character(current.pop) Table.2[i, 2]<-as.numeric(as.character(N)) Table.2[i, 3]<-as.numeric(as.character(Mean.Fit.40.to.60)) Table.2[i, 4]<-as.numeric(as.character(SD.Fit.40.to.60)) Table.2[i, 5]<-as.numeric(as.character(t.40.to.60)) Table.2[i, 6]<-as.numeric(as.character(p.40.to.60)) Table.2[i, 7]<-as.numeric(as.character(Mean.Fit.40.to.50)) Table.2[i, 8]<-as.numeric(as.character(SD.Fit.40.to.50)) Table.2[i, 9]<-as.numeric(as.character(t.40.to.50)) Table.2[i, 10]<-as.numeric(as.character(p.40.to.50)) Table.2[i, 11]<-as.numeric(as.character(Mean.Fit.50.to.60)) Table.2[i, 12]<-as.numeric(as.character(SD.Fit.50.to.60)) Table.2[i, 13]<-as.numeric(as.character(t.50.to.60)) Table.2[i, 14]<-as.numeric(as.character(p.50.to.60)) } # Now creating a version with extra digits in case I need more precision for # later analyses. Table.2[,3]<-round(as.numeric(as.character(Table.2[,3])), digits=7) Table.2[,4]<-round(as.numeric(as.character(Table.2[,4])), digits=7) Table.2[,5]<-round(as.numeric(as.character(Table.2[,5])), digits=7) Table.2[,6]<-round(as.numeric(as.character(Table.2[,6])), digits=18) Table.2[,7]<-round(as.numeric(as.character(Table.2[,7])), digits=7) Table.2[,8]<-round(as.numeric(as.character(Table.2[,8])), digits=7) Table.2[,9]<-round(as.numeric(as.character(Table.2[,9])), digits=7) Table.2[,10]<-round(as.numeric(as.character(Table.2[,10])), digits=18) Table.2[,11]<-round(as.numeric(as.character(Table.2[,11])), digits=7) Table.2[,12]<-round(as.numeric(as.character(Table.2[,12])), digits=7) Table.2[,13]<-round(as.numeric(as.character(Table.2[,13])), digits=7) Table.2[,14]<-round(as.numeric(as.character(Table.2[,14])), digits=18) Table.2.extra.digits<-Table.2 # Now we round these off for more useful display. Table.2[,3]<-round(as.numeric(as.character(Table.2[,3])), digits=4) Table.2[,4]<-round(as.numeric(as.character(Table.2[,4])), digits=4) Table.2[,5]<-round(as.numeric(as.character(Table.2[,5])), digits=4) Table.2[,6]<-round(as.numeric(as.character(Table.2[,6])), digits=18) Table.2[,7]<-round(as.numeric(as.character(Table.2[,7])), digits=4) Table.2[,8]<-round(as.numeric(as.character(Table.2[,8])), digits=4) Table.2[,9]<-round(as.numeric(as.character(Table.2[,9])), digits=4) Table.2[,10]<-round(as.numeric(as.character(Table.2[,10])), digits=18) Table.2[,11]<-round(as.numeric(as.character(Table.2[,11])), digits=4) Table.2[,12]<-round(as.numeric(as.character(Table.2[,12])), digits=4) Table.2[,13]<-round(as.numeric(as.character(Table.2[,13])), digits=4) Table.2[,14]<-round(as.numeric(as.character(Table.2[,14])), digits=18) # This gives us a table with uncorrected p values. We do have 27 tests in this table, # so now I'm going to do a sequential Bonferroni correction — take the most significant # value and adjust for the fact there are 27 tests; then the next most and adjust for # the 26 remaining, etc. # First, pull all the p values into one listing p.value.matrix<-matrix(nrow=27, ncol=4) for (i in 1:9) { p.value.matrix[i, 1]<-Table.2[i, 1] p.value.matrix[i, 2]<-as.numeric(40000) p.value.matrix[i, 3]<-Table.2[i, 6] p.value.matrix[i+9, 1]<-Table.2[i, 1] p.value.matrix[i+9, 2]<-as.numeric(50000) p.value.matrix[i+9, 3]<-Table.2[i, 10] p.value.matrix[i+18, 1]<-Table.2[i, 1] p.value.matrix[i+18, 2]<-as.numeric(60000) p.value.matrix[i+18, 3]<-Table.2[i, 14] } # Then sort them from smallest to largest p.value.matrix<-p.value.matrix[order(as.numeric(as.character(p.value.matrix[,3]))),] # Then create the sequential corrections for (i in 1:27) { current.p<-as.numeric(as.character(p.value.matrix[i, 3])) adjusted.p<-current.p*(28-i) p.value.matrix[i, 4]<-as.numeric(as.character(adjusted.p)) } # Then sort them back into the original ordering. p.value.matrix<-p.value.matrix[order(as.numeric(as.character(p.value.matrix[,2]))),] ordered.p.value.matrix<-matrix(nrow=0, ncol=4) for (i in 1:length(Population.list)) { current.pop<-Population.list[i] current.matrix<-subset(p.value.matrix, p.value.matrix[,1]==current.pop) ordered.p.value.matrix<-rbind(ordered.p.value.matrix, current.matrix) } ordered.p.value.matrix<-ordered.p.value.matrix[order(as.numeric(as.character(ordered.p.value.matrix[,2]))),] # Now we can create an alternate version of Table 2, corrected for multiple comparisons corrected.p.40<-ordered.p.value.matrix[1:9, 4] corrected.p.50<-ordered.p.value.matrix[10:18, 4] corrected.p.60<-ordered.p.value.matrix[19:27, 4] Table.2.alternate<-Table.2 Table.2.alternate[,6]<-corrected.p.40 Table.2.alternate[,10]<-corrected.p.50 Table.2.alternate[,14]<-corrected.p.60 # Round all of the p values to 4 decimal places: Table.2[,6]<-round(as.numeric(as.character(Table.2[,6])), digits=4) Table.2[,10]<-round(as.numeric(as.character(Table.2[,10])), digits=4) Table.2[,14]<-round(as.numeric(as.character(Table.2[,14])), digits=4) Table.2.alternate[,6]<-round(as.numeric(as.character(Table.2.alternate[,6])), digits=4) Table.2.alternate[,10]<-round(as.numeric(as.character(Table.2.alternate[,10])), digits=4) Table.2.alternate[,14]<-round(as.numeric(as.character(Table.2.alternate[,14])), digits=4) # Replace very small p values with a < 0.0001 marker for (i in 1: length(Table.2[,6])) { value.1<-as.numeric(as.character(Table.2[i,6])) value.2<-as.numeric(as.character(Table.2[i,10])) value.3<-as.numeric(as.character(Table.2[i,14])) if (value.1<0.0001) { Table.2[i, 6]<-"< 0.0001" } if (value.2<0.0001) { Table.2[i, 10]<-"< 0.0001" } if (value.3<0.0001) { Table.2[i, 14]<-"< 0.0001" } } colnames(Table.2)<-c("Population", "N", "Mean", "SD", "t", "p", "Mean", "SD", "t", "p", "Mean", "SD", "t", "p") Table.2<-as.data.frame(Table.2) for (i in 1: length(Table.2.alternate[,6])) { value.1<-as.numeric(as.character(Table.2.alternate[i,6])) value.2<-as.numeric(as.character(Table.2.alternate[i,10])) value.3<-as.numeric(as.character(Table.2.alternate[i,14])) if (value.1<0.0001) { Table.2.alternate[i, 6]<-"< 0.0001" } if (value.2<0.0001) { Table.2.alternate[i, 10]<-"< 0.0001" } if (value.3<0.0001) { Table.2.alternate[i, 14]<-"< 0.0001" } } colnames(Table.2.alternate)<-c("Population", "N", "Mean", "SD", "t", "p - adj", "Mean", "SD", "t", "p - adj", "Mean", "SD", "t", "p - adj") Table.2.alternate<-as.data.frame(Table.2.alternate) # Like with Table 1, I don't know how to make a column name span across columns, so I'll need to do that in # a spreadsheet application afterwards. # Tables were renumbered on submission, as many were moved to the Supplement. What used to be Table 2 is now Table 1. # We used the alternate version, which involves the sequential Bonferroni correction. #write.csv(file="Table.2.csv", Table.2) write.csv(file="Table.1.csv", Table.2.alternate) ########################## # Now for Table 3. Thankfully, we can pull some of this data straight from Table 2. Table.3<-matrix(ncol=4, nrow=9) Table.3[,1]<-as.character(Table.2[,1]) Table.3[,2]<-as.numeric(as.character(Table.2.extra.digits[,3])) Table.3[,3]<-as.numeric(as.character(Table.2.extra.digits[,7])) Table.3[,4]<-as.numeric(as.character(Table.2.extra.digits[,11])) # Find the means and standard deviations among population points Table.3.grand<-matrix(ncol=4, nrow=2) Table.3.grand[1, 1]<-"Grand Mean" Table.3.grand[1, 2]<-mean(as.numeric(as.character(Table.3[, 2]))) Table.3.grand[1, 3]<-mean(as.numeric(as.character(Table.3[, 3]))) Table.3.grand[1, 4]<-mean(as.numeric(as.character(Table.3[, 4]))) Table.3.grand[2, 1]<-"SD" Table.3.grand[2, 2]<-sd(as.numeric(as.character(Table.3[, 2]))) Table.3.grand[2, 3]<-sd(as.numeric(as.character(Table.3[, 3]))) Table.3.grand[2, 4]<-sd(as.numeric(as.character(Table.3[, 4]))) # Find the exptected values from the hyperbolic model Table.3.predicted.hyper<-matrix(nrow=1, ncol=4) Table.3.predicted.hyper[1,1]<-"Predicted by hypberbolic model" Table.3.predicted.hyper[1,2]<-as.numeric((1.6525/1.6308)-1) Table.3.predicted.hyper[1,3]<-as.numeric((1.6436/1.6308)-1) Table.3.predicted.hyper[1,4]<-as.numeric((1.6525/1.6437)-1) # Find whether values are different from hyperbolic model Table.3.results.hyper<-matrix(nrow=1, ncol=4) Table.3.results.hyper[1, 1]<-"Two-tailed p" for (i in 2:4) { test.data<-as.numeric(as.character(Table.3[,i])) prediction<-as.numeric(as.character(Table.3.predicted.hyper[1, i])) result<-as.numeric(as.character(t.test(test.data, mu=prediction, alternative="two.sided")[3])) Table.3.results.hyper[1, i]<-as.numeric(as.character(result)) } # Find predictions of power law model Table.3.predicted.power<-matrix(nrow=1, ncol=4) Table.3.predicted.power[1,1]<-"Predicted by power-law model" Table.3.predicted.power[1,2]<-as.numeric((1.7246/1.6597)-1) Table.3.predicted.power[1,3]<-as.numeric((1.6951/1.6597)-1) Table.3.predicted.power[1,4]<-as.numeric((1.7246/1.6951)-1) # Find whether results are different from power law model Table.3.results.power<-matrix(nrow=1, ncol=4) Table.3.results.power[1, 1]<-"Two-tailed p" for (i in 2:4) { test.data<-as.numeric(as.character(Table.3[,i])) prediction<-as.numeric(as.character(Table.3.predicted.power[1, i])) result<-as.numeric(as.character(t.test(test.data, mu=prediction, alternative="two.sided")[3])) Table.3.results.power[1, i]<-as.numeric(as.character(result)) } # Stitch components together Table.3<-rbind(Table.3, Table.3.grand, Table.3.predicted.hyper, Table.3.results.hyper, Table.3.predicted.power, Table.3.results.power) # Round values for convenience in display Table.3[,2]<-round(as.numeric(as.character(Table.3[,2])), digits=4) Table.3[,3]<-round(as.numeric(as.character(Table.3[,3])), digits=4) Table.3[,4]<-round(as.numeric(as.character(Table.3[,4])), digits=4) # Replace extremely small p values with < 0.0001 row.list<-c(13, 15) for (i in 1:2) { j<-row.list[i] value.1<-as.numeric(as.character(Table.3[j,2])) value.2<-as.numeric(as.character(Table.3[j,3])) value.3<-as.numeric(as.character(Table.3[j,4])) if (value.1<0.0001) { Table.3[j, 2]<-"< 0.0001" } if (value.2<0.0001) { Table.3[j, 3]<-"< 0.0001" } if (value.3<0.0001) { Table.3[j, 4]<-"< 0.0001" } } # Name the columns: colnames(Table.3)<-c("Population", "40,000 to 60,000", "40,000 to 50,000", "50,000 to 60,000") Table.3<-as.data.frame(Table.3) # Tables were renumbered on submission, as many were moved to the Supplement. What used to be at the bottom of Table 3 is now # Table 2. #write.csv(file="Table.2.csv", Table.3) Alt.Table.2<-Table.3 Alt.Table.2<-Alt.Table.2[10:15,] write.csv(file="Table.2.csv", Alt.Table.2) # The following lines are used to generate a statistic reported in the main text: observed.early<-Table.3[1:9, 3] observed.late<-Table.3[1:9, 4] observed.early<-as.numeric(as.character(observed.early)) observed.late<-as.numeric(as.character(observed.late)) expected.late<-observed.early-0.0039 diff.late<-expected.late-observed.late t.test(diff.late, alternative=c("two.sided")) diff.2.late<-observed.early-observed.late t.test(diff.2.late, mu=0.0039, alternative=c("two.sided")) t.test(diff.late, alternative=c("two.sided"), mu=0.0039) ############################### # For Tables 4-6, I need a dummy variable for replicate which counts each replicate differently, regardless # of which person did it. This is essentially Person.Plus.Replicate from data.accepted, but as a factor. ###################### # Now I will try to implement the Friedman method nonparametric analog to ANOVAs for tables 4, 5, and 6. # Let's try this as a function, so I can just call it when needed. Friedman.p<-function(data.set) { Rep.list<-unique(data.set$Dummy.Replicate) Pop.list<-unique(data.set$Evolved.Population) rank.matrix<-matrix(nrow=length(Pop.list), ncol=length(Rep.list)) for (i in 1:length(Rep.list)) { current.rep<-Rep.list[i] current.data<-subset(data.set, Dummy.Replicate==current.rep) Flask.list<-current.data$Flask Fit.list<-current.data$Fitness current.matrix<-matrix(ncol=3, nrow=length(Pop.list)) current.matrix[,1]<-as.numeric(as.character(Flask.list)) current.matrix[,2]<-as.numeric(as.character(Fit.list)) current.matrix[,2]<-as.numeric(as.character(current.matrix[,2]))+1000 current.matrix<-current.matrix[order(as.numeric(as.character(current.matrix[,2]))),] current.matrix[,3]<-seq(1:length(Pop.list)) current.matrix<-current.matrix[order(current.matrix[,1]),] current.matrix[,2]<-as.numeric(as.character(current.matrix[,2]))-1000 rank.matrix[,i]<-as.numeric(as.character(current.matrix[,3])) } rank.sum<-matrix(nrow=length(Pop.list), ncol=1) for (i in 1:length(Pop.list)) { current.data<-as.numeric(as.character(rank.matrix[i,])) current.sum<-sum(as.numeric(as.character(current.data))) rank.sum[i, 1]<-current.sum } a<-length(Pop.list) b<-length(Rep.list) squared.sum.ranks<-0 for (i in 1:length(Pop.list)) { current.rank.sum<-as.numeric(as.character(rank.sum[i,1])) value<-as.numeric(as.character(current.rank.sum))^2 squared.sum.ranks<-as.numeric(as.character(squared.sum.ranks))+as.numeric(as.character(value)) } computed.value<-((12/(a*b*(a+1)))*(as.numeric(as.character(squared.sum.ranks)))) - (3*b*(a+1)) calc.p<-1-pchisq(computed.value, (a-1)) return(calc.p) } ### data.accepted$Dummy.Replicate<-factor(data.accepted$Person.Plus.Replicate) Dummy.list<-unique(data.accepted$Dummy.Replicate) # Next, we need to drop replicates where not all the populations had useable data at a given generation. # We do this separately for the three different generations. data.for.ANOVA.40<-c() for (i in 1:length(Dummy.list)) { current.rep<-Dummy.list[i] current.data<-subset(data.accepted, Dummy.Replicate==current.rep) current.data<-subset(current.data, Generation==40000) current.length<-length(current.data$Fitness) if (current.length==9) { data.for.ANOVA.40<-rbind(data.for.ANOVA.40, current.data) } } data.for.ANOVA.50<-c() for (i in 1:length(Dummy.list)) { current.rep<-Dummy.list[i] current.data<-subset(data.accepted, Dummy.Replicate==current.rep) current.data<-subset(current.data, Generation==50000) current.length<-length(current.data$Fitness) if (current.length==9) { data.for.ANOVA.50<-rbind(data.for.ANOVA.50, current.data) } } data.for.ANOVA.60<-c() for (i in 1:length(Dummy.list)) { current.rep<-Dummy.list[i] current.data<-subset(data.accepted, Dummy.Replicate==current.rep) current.data<-subset(current.data, Generation==60000) current.length<-length(current.data$Fitness) if (current.length==9) { data.for.ANOVA.60<-rbind(data.for.ANOVA.60, current.data) } } # Now we look at dropping certain populations: data.for.ANOVA.40.6<-subset(data.for.ANOVA.40, Evolved.Population=="Ara + 1" | Evolved.Population=="Ara + 2" | Evolved.Population=="Ara + 4" | Evolved.Population=="Ara + 5" | Evolved.Population=="Ara - 5" | Evolved.Population=="Ara - 6") data.for.ANOVA.40.5<-subset(data.for.ANOVA.40.6, Evolved.Population=="Ara + 2" | Evolved.Population=="Ara + 4" | Evolved.Population=="Ara + 5" | Evolved.Population=="Ara - 5" | Evolved.Population=="Ara - 6" ) data.for.ANOVA.50.6<-subset(data.for.ANOVA.50, Evolved.Population=="Ara + 1" | Evolved.Population=="Ara + 2" | Evolved.Population=="Ara + 4" | Evolved.Population=="Ara + 5" | Evolved.Population=="Ara - 5" | Evolved.Population=="Ara - 6") data.for.ANOVA.50.5<-subset(data.for.ANOVA.50.6, Evolved.Population=="Ara + 2" | Evolved.Population=="Ara + 4" | Evolved.Population=="Ara + 5" | Evolved.Population=="Ara - 5" | Evolved.Population=="Ara - 6" ) data.for.ANOVA.60.6<-subset(data.for.ANOVA.60, Evolved.Population=="Ara + 1" | Evolved.Population=="Ara + 2" | Evolved.Population=="Ara + 4" | Evolved.Population=="Ara + 5" | Evolved.Population=="Ara - 5" | Evolved.Population=="Ara - 6") data.for.ANOVA.60.5<-subset(data.for.ANOVA.60.6, Evolved.Population=="Ara + 2" | Evolved.Population=="Ara + 4" | Evolved.Population=="Ara + 5" | Evolved.Population=="Ara - 5" | Evolved.Population=="Ara - 6" ) # Now we can calculate the ANOVAs: anova.full.40<-aov(Fitness ~ Error(Evolved.Population*Dummy.Replicate), data=data.for.ANOVA.40) anova.6.40<-aov(Fitness ~ Error(Evolved.Population*Dummy.Replicate), data=data.for.ANOVA.40.6) anova.5.40<-aov(Fitness ~ Error(Evolved.Population*Dummy.Replicate), data=data.for.ANOVA.40.5) anova.full.50<-aov(Fitness ~ Error(Evolved.Population*Dummy.Replicate), data=data.for.ANOVA.50) anova.6.50<-aov(Fitness ~ Error(Evolved.Population*Dummy.Replicate), data=data.for.ANOVA.50.6) anova.5.50<-aov(Fitness ~ Error(Evolved.Population*Dummy.Replicate), data=data.for.ANOVA.50.5) anova.full.60<-aov(Fitness ~ Error(Evolved.Population*Dummy.Replicate), data=data.for.ANOVA.60) anova.6.60<-aov(Fitness ~ Error(Evolved.Population*Dummy.Replicate), data=data.for.ANOVA.60.6) anova.5.60<-aov(Fitness ~ Error(Evolved.Population*Dummy.Replicate), data=data.for.ANOVA.60.5) F.p.40.full<-Friedman.p(data.for.ANOVA.40) F.p.40.6<-Friedman.p(data.for.ANOVA.40.6) F.p.40.5<-Friedman.p(data.for.ANOVA.40.5) F.p.50.full<-Friedman.p(data.for.ANOVA.50) F.p.50.6<-Friedman.p(data.for.ANOVA.50.6) F.p.50.5<-Friedman.p(data.for.ANOVA.50.5) F.p.60.full<-Friedman.p(data.for.ANOVA.60) F.p.60.6<-Friedman.p(data.for.ANOVA.60.6) F.p.60.5<-Friedman.p(data.for.ANOVA.60.5) # Unfortunately, with random effects, R does not display the ANOVA tables in ways like # we're used to. Here I extract the relevant information of Sum of Squares, degrees of freedom # and Mean Square from what they give me; calculate the F values from those; and calculate # p values from the F statistics. # First the ANOVA table for all populations at 40,000 generations: anova.table.40.full<-matrix(nrow=6, ncol=7) anova.table.40.full[1,1]<-"All nine populations" anova.table.40.full[2,1]<-"Source" anova.table.40.full[2,2]<-"SS" anova.table.40.full[2,3]<-"df" anova.table.40.full[2,4]<-"MS" anova.table.40.full[2,5]<-"F" anova.table.40.full[2,6]<-"p" anova.table.40.full[2,7]<-"Friedman's p" anova.table.40.full[3,1]<-"Population" anova.table.40.full[3,2]<-round(as.numeric(as.character(summary(anova.full.40)[[1]][[1]][[2]])), digits=5) anova.table.40.full[3,3]<-as.numeric(as.character(summary(anova.full.40)[[1]][[1]][[1]])) anova.table.40.full[3,4]<-round(as.numeric(as.character(summary(anova.full.40)[[1]][[1]][[3]])), digits=5) anova.table.40.full[4,1]<-"Block" anova.table.40.full[4,2]<-round(as.numeric(as.character(summary(anova.full.40)[[2]][[1]][[2]])), digits=5) anova.table.40.full[4,3]<-as.numeric(as.character(summary(anova.full.40)[[2]][[1]][[1]])) anova.table.40.full[4,4]<-round(as.numeric(as.character(summary(anova.full.40)[[2]][[1]][[3]])), digits=5) anova.table.40.full[5,1]<-"Error" anova.table.40.full[5,2]<-round(as.numeric(as.character(summary(anova.full.40)[[3]][[1]][[2]])), digits=5) anova.table.40.full[5,3]<-as.numeric(as.character(summary(anova.full.40)[[3]][[1]][[1]])) anova.table.40.full[5,4]<-round(as.numeric(as.character(summary(anova.full.40)[[3]][[1]][[3]])), digits=5) anova.table.40.full[6,1]<-"Total" anova.table.40.full[3,5]<-round(as.numeric(as.character(summary(anova.full.40)[[1]][[1]][[3]]))/as.numeric(as.character(summary(anova.full.40)[[3]][[1]][[3]])), digits=5) anova.table.40.full[3,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.40.full[3,5])), as.numeric(as.character(summary(anova.full.40)[[1]][[1]][[1]])), as.numeric(as.character(summary(anova.full.40)[[3]][[1]][[1]])) ))), digits=4) anova.table.40.full[3,7]<-F.p.40.full anova.table.40.full[4,5]<-round(as.numeric(as.character(summary(anova.full.40)[[2]][[1]][[3]]))/as.numeric(as.character(summary(anova.full.40)[[3]][[1]][[3]])), digits=5) anova.table.40.full[4,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.40.full[4,5])), as.numeric(as.character(summary(anova.full.40)[[2]][[1]][[1]])), as.numeric(as.character(summary(anova.full.40)[[3]][[1]][[1]])) ))), digits=4) if(as.numeric(as.character(anova.table.40.full[3,6]))<0.0001) { anova.table.40.full[3,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.40.full[4,6]))<0.0001) { anova.table.40.full[4,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.40.full[3,7]))<0.0001) { anova.table.40.full[3,7]<-"< 0.0001"} anova.table.40.full[6,2]<-as.numeric(as.character(anova.table.40.full[3,2]))+as.numeric(as.character(anova.table.40.full[4,2]))+as.numeric(as.character(anova.table.40.full[5,2])) anova.table.40.full[6,3]<-as.numeric(as.character(anova.table.40.full[3,3]))+as.numeric(as.character(anova.table.40.full[4,3]))+as.numeric(as.character(anova.table.40.full[5,3])) anova.table.40.full[6,4]<-round(as.numeric(as.character(anova.table.40.full[6,2]))/as.numeric(as.character(anova.table.40.full[6,3])), digits=5) # Then for just the 6 non-hypermutator populations: anova.table.40.6<-matrix(nrow=6, ncol=7) anova.table.40.6[1,1]<-"Excluding three hypermutator populations" anova.table.40.6[2,1]<-"Source" anova.table.40.6[2,2]<-"SS" anova.table.40.6[2,3]<-"df" anova.table.40.6[2,4]<-"MS" anova.table.40.6[2,5]<-"F" anova.table.40.6[2,6]<-"p" anova.table.40.6[2,7]<-"Friedman's p" anova.table.40.6[3,1]<-"Population" anova.table.40.6[3,2]<-round(as.numeric(as.character(summary(anova.6.40)[[1]][[1]][[2]])), digits=5) anova.table.40.6[3,3]<-as.numeric(as.character(summary(anova.6.40)[[1]][[1]][[1]])) anova.table.40.6[3,4]<-round(as.numeric(as.character(summary(anova.6.40)[[1]][[1]][[3]])), digits=5) anova.table.40.6[4,1]<-"Block" anova.table.40.6[4,2]<-round(as.numeric(as.character(summary(anova.6.40)[[2]][[1]][[2]])), digits=5) anova.table.40.6[4,3]<-as.numeric(as.character(summary(anova.6.40)[[2]][[1]][[1]])) anova.table.40.6[4,4]<-round(as.numeric(as.character(summary(anova.6.40)[[2]][[1]][[3]])), digits=5) anova.table.40.6[5,1]<-"Error" anova.table.40.6[5,2]<-round(as.numeric(as.character(summary(anova.6.40)[[3]][[1]][[2]])), digits=5) anova.table.40.6[5,3]<-as.numeric(as.character(summary(anova.6.40)[[3]][[1]][[1]])) anova.table.40.6[5,4]<-round(as.numeric(as.character(summary(anova.6.40)[[3]][[1]][[3]])), digits=5) anova.table.40.6[6,1]<-"Total" anova.table.40.6[3,5]<-round(as.numeric(as.character(summary(anova.6.40)[[1]][[1]][[3]]))/as.numeric(as.character(summary(anova.6.40)[[3]][[1]][[3]])), digits=5) anova.table.40.6[3,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.40.6[3,5])), as.numeric(as.character(summary(anova.6.40)[[1]][[1]][[1]])), as.numeric(as.character(summary(anova.6.40)[[3]][[1]][[1]])) ))), digits=4) anova.table.40.6[3,7]<-F.p.40.6 anova.table.40.6[4,5]<-round(as.numeric(as.character(summary(anova.6.40)[[2]][[1]][[3]]))/as.numeric(as.character(summary(anova.6.40)[[3]][[1]][[3]])), digits=5) anova.table.40.6[4,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.40.6[4,5])), as.numeric(as.character(summary(anova.6.40)[[2]][[1]][[1]])), as.numeric(as.character(summary(anova.6.40)[[3]][[1]][[1]])) ))), digits=4) if(as.numeric(as.character(anova.table.40.6[3,6]))<0.0001) { anova.table.40.6[3,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.40.6[4,6]))<0.0001) { anova.table.40.6[4,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.40.6[3,7]))<0.0001) { anova.table.40.6[3,7]<-"< 0.0001"} anova.table.40.6[6,2]<-as.numeric(as.character(anova.table.40.6[3,2]))+as.numeric(as.character(anova.table.40.6[4,2]))+as.numeric(as.character(anova.table.40.6[5,2])) anova.table.40.6[6,3]<-as.numeric(as.character(anova.table.40.6[3,3]))+as.numeric(as.character(anova.table.40.6[4,3]))+as.numeric(as.character(anova.table.40.6[5,3])) anova.table.40.6[6,4]<-round(as.numeric(as.character(anova.table.40.6[6,2]))/as.numeric(as.character(anova.table.40.6[6,3])), digits=5) # And then again dropping Ara + 1 anova.table.40.5<-matrix(nrow=6, ncol=7) anova.table.40.5[1,1]<-"Also exclusing population Ara+1" anova.table.40.5[2,1]<-"Source" anova.table.40.5[2,2]<-"SS" anova.table.40.5[2,3]<-"df" anova.table.40.5[2,4]<-"MS" anova.table.40.5[2,5]<-"F" anova.table.40.5[2,6]<-"p" anova.table.40.5[2,7]<-"Friedman.p" anova.table.40.5[3,1]<-"Population" anova.table.40.5[3,2]<-round(as.numeric(as.character(summary(anova.5.40)[[1]][[1]][[2]])), digits=5) anova.table.40.5[3,3]<-as.numeric(as.character(summary(anova.5.40)[[1]][[1]][[1]])) anova.table.40.5[3,4]<-round(as.numeric(as.character(summary(anova.5.40)[[1]][[1]][[3]])), digits=5) anova.table.40.5[4,1]<-"Block" anova.table.40.5[4,2]<-round(as.numeric(as.character(summary(anova.5.40)[[2]][[1]][[2]])), digits=5) anova.table.40.5[4,3]<-as.numeric(as.character(summary(anova.5.40)[[2]][[1]][[1]])) anova.table.40.5[4,4]<-round(as.numeric(as.character(summary(anova.5.40)[[2]][[1]][[3]])), digits=5) anova.table.40.5[5,1]<-"Error" anova.table.40.5[5,2]<-round(as.numeric(as.character(summary(anova.5.40)[[3]][[1]][[2]])), digits=5) anova.table.40.5[5,3]<-as.numeric(as.character(summary(anova.5.40)[[3]][[1]][[1]])) anova.table.40.5[5,4]<-round(as.numeric(as.character(summary(anova.5.40)[[3]][[1]][[3]])), digits=5) anova.table.40.5[6,1]<-"Total" anova.table.40.5[3,5]<-round(as.numeric(as.character(summary(anova.5.40)[[1]][[1]][[3]]))/as.numeric(as.character(summary(anova.5.40)[[3]][[1]][[3]])), digits=5) anova.table.40.5[3,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.40.5[3,5])), as.numeric(as.character(summary(anova.5.40)[[1]][[1]][[1]])), as.numeric(as.character(summary(anova.5.40)[[3]][[1]][[1]])) ))), digits=4) anova.table.40.5[3,7]<-F.p.40.5 anova.table.40.5[4,5]<-round(as.numeric(as.character(summary(anova.5.40)[[2]][[1]][[3]]))/as.numeric(as.character(summary(anova.5.40)[[3]][[1]][[3]])), digits=5) anova.table.40.5[4,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.40.5[4,5])), as.numeric(as.character(summary(anova.5.40)[[2]][[1]][[1]])), as.numeric(as.character(summary(anova.5.40)[[3]][[1]][[1]])) ))), digits=4) if(as.numeric(as.character(anova.table.40.5[3,6]))<0.0001) { anova.table.40.5[3,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.40.5[4,6]))<0.0001) { anova.table.40.5[4,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.40.5[3,7]))<0.0001) { anova.table.40.5[3,7]<-"< 0.0001"} anova.table.40.5[6,2]<-as.numeric(as.character(anova.table.40.5[3,2]))+as.numeric(as.character(anova.table.40.5[4,2]))+as.numeric(as.character(anova.table.40.5[5,2])) anova.table.40.5[6,3]<-as.numeric(as.character(anova.table.40.5[3,3]))+as.numeric(as.character(anova.table.40.5[4,3]))+as.numeric(as.character(anova.table.40.5[5,3])) anova.table.40.5[6,4]<-round(as.numeric(as.character(anova.table.40.5[6,2]))/as.numeric(as.character(anova.table.40.5[6,3])), digits=5) # And now we stitch these together to form Table 4: Table.4<-rbind(anova.table.40.full, anova.table.40.6, anova.table.40.5) Table.4<-as.data.frame(Table.4) # Tables were renumbered on submission, as many were moved to the Supplement. What used to be Table 4 is now Table S3. write.csv(file="Table.S3.csv", Table.4) ########################## # Table 5 repeats the procedure for Table 4, just on the data from 50,000 generations. anova.table.50.full<-matrix(nrow=6, ncol=7) anova.table.50.full[1,1]<-"All nine populations" anova.table.50.full[2,1]<-"Source" anova.table.50.full[2,2]<-"SS" anova.table.50.full[2,3]<-"df" anova.table.50.full[2,4]<-"MS" anova.table.50.full[2,5]<-"F" anova.table.50.full[2,6]<-"p" anova.table.50.full[2,7]<-"Friedman's p" anova.table.50.full[3,1]<-"Population" anova.table.50.full[3,2]<-round(as.numeric(as.character(summary(anova.full.50)[[1]][[1]][[2]])), digits=5) anova.table.50.full[3,3]<-as.numeric(as.character(summary(anova.full.50)[[1]][[1]][[1]])) anova.table.50.full[3,4]<-round(as.numeric(as.character(summary(anova.full.50)[[1]][[1]][[3]])), digits=5) anova.table.50.full[4,1]<-"Block" anova.table.50.full[4,2]<-round(as.numeric(as.character(summary(anova.full.50)[[2]][[1]][[2]])), digits=5) anova.table.50.full[4,3]<-as.numeric(as.character(summary(anova.full.50)[[2]][[1]][[1]])) anova.table.50.full[4,4]<-round(as.numeric(as.character(summary(anova.full.50)[[2]][[1]][[3]])), digits=5) anova.table.50.full[5,1]<-"Error" anova.table.50.full[5,2]<-round(as.numeric(as.character(summary(anova.full.50)[[3]][[1]][[2]])), digits=5) anova.table.50.full[5,3]<-as.numeric(as.character(summary(anova.full.50)[[3]][[1]][[1]])) anova.table.50.full[5,4]<-round(as.numeric(as.character(summary(anova.full.50)[[3]][[1]][[3]])), digits=5) anova.table.50.full[6,1]<-"Total" anova.table.50.full[3,5]<-round(as.numeric(as.character(summary(anova.full.50)[[1]][[1]][[3]]))/as.numeric(as.character(summary(anova.full.50)[[3]][[1]][[3]])), digits=5) anova.table.50.full[3,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.50.full[3,5])), as.numeric(as.character(summary(anova.full.50)[[1]][[1]][[1]])), as.numeric(as.character(summary(anova.full.50)[[3]][[1]][[1]])) ))), digits=4) anova.table.50.full[3,7]<-F.p.50.full anova.table.50.full[4,5]<-round(as.numeric(as.character(summary(anova.full.50)[[2]][[1]][[3]]))/as.numeric(as.character(summary(anova.full.50)[[3]][[1]][[3]])), digits=5) anova.table.50.full[4,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.50.full[4,5])), as.numeric(as.character(summary(anova.full.50)[[2]][[1]][[1]])), as.numeric(as.character(summary(anova.full.50)[[3]][[1]][[1]])) ))), digits=4) if(as.numeric(as.character(anova.table.50.full[3,6]))<0.0001) { anova.table.50.full[3,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.50.full[4,6]))<0.0001) { anova.table.50.full[4,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.50.full[3,7]))<0.0001) { anova.table.50.full[3,7]<-"< 0.0001"} anova.table.50.full[6,2]<-as.numeric(as.character(anova.table.50.full[3,2]))+as.numeric(as.character(anova.table.50.full[4,2]))+as.numeric(as.character(anova.table.50.full[5,2])) anova.table.50.full[6,3]<-as.numeric(as.character(anova.table.50.full[3,3]))+as.numeric(as.character(anova.table.50.full[4,3]))+as.numeric(as.character(anova.table.50.full[5,3])) anova.table.50.full[6,4]<-round(as.numeric(as.character(anova.table.50.full[6,2]))/as.numeric(as.character(anova.table.50.full[6,3])), digits=5) # Then for just the 6 non-hypermutator populations: anova.table.50.6<-matrix(nrow=6, ncol=7) anova.table.50.6[1,1]<-"Excluding three hypermutator populations" anova.table.50.6[2,1]<-"Source" anova.table.50.6[2,2]<-"SS" anova.table.50.6[2,3]<-"df" anova.table.50.6[2,4]<-"MS" anova.table.50.6[2,5]<-"F" anova.table.50.6[2,6]<-"p" anova.table.50.6[2,7]<-"Friedman's p" anova.table.50.6[3,1]<-"Population" anova.table.50.6[3,2]<-round(as.numeric(as.character(summary(anova.6.50)[[1]][[1]][[2]])), digits=5) anova.table.50.6[3,3]<-as.numeric(as.character(summary(anova.6.50)[[1]][[1]][[1]])) anova.table.50.6[3,4]<-round(as.numeric(as.character(summary(anova.6.50)[[1]][[1]][[3]])), digits=5) anova.table.50.6[4,1]<-"Block" anova.table.50.6[4,2]<-round(as.numeric(as.character(summary(anova.6.50)[[2]][[1]][[2]])), digits=5) anova.table.50.6[4,3]<-as.numeric(as.character(summary(anova.6.50)[[2]][[1]][[1]])) anova.table.50.6[4,4]<-round(as.numeric(as.character(summary(anova.6.50)[[2]][[1]][[3]])), digits=5) anova.table.50.6[5,1]<-"Error" anova.table.50.6[5,2]<-round(as.numeric(as.character(summary(anova.6.50)[[3]][[1]][[2]])), digits=5) anova.table.50.6[5,3]<-as.numeric(as.character(summary(anova.6.50)[[3]][[1]][[1]])) anova.table.50.6[5,4]<-round(as.numeric(as.character(summary(anova.6.50)[[3]][[1]][[3]])), digits=5) anova.table.50.6[6,1]<-"Total" anova.table.50.6[3,5]<-round(as.numeric(as.character(summary(anova.6.50)[[1]][[1]][[3]]))/as.numeric(as.character(summary(anova.6.50)[[3]][[1]][[3]])), digits=5) anova.table.50.6[3,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.50.6[3,5])), as.numeric(as.character(summary(anova.6.50)[[1]][[1]][[1]])), as.numeric(as.character(summary(anova.6.50)[[3]][[1]][[1]])) ))), digits=4) anova.table.50.6[3,7]<-F.p.50.6 anova.table.50.6[4,5]<-round(as.numeric(as.character(summary(anova.6.50)[[2]][[1]][[3]]))/as.numeric(as.character(summary(anova.6.50)[[3]][[1]][[3]])), digits=5) anova.table.50.6[4,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.50.6[4,5])), as.numeric(as.character(summary(anova.6.50)[[2]][[1]][[1]])), as.numeric(as.character(summary(anova.6.50)[[3]][[1]][[1]])) ))), digits=4) if(as.numeric(as.character(anova.table.50.6[3,6]))<0.0001) { anova.table.50.6[3,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.50.6[4,6]))<0.0001) { anova.table.50.6[4,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.50.6[3,7]))<0.0001) { anova.table.50.6[3,7]<-"< 0.0001"} anova.table.50.6[6,2]<-as.numeric(as.character(anova.table.50.6[3,2]))+as.numeric(as.character(anova.table.50.6[4,2]))+as.numeric(as.character(anova.table.50.6[5,2])) anova.table.50.6[6,3]<-as.numeric(as.character(anova.table.50.6[3,3]))+as.numeric(as.character(anova.table.50.6[4,3]))+as.numeric(as.character(anova.table.50.6[5,3])) anova.table.50.6[6,4]<-round(as.numeric(as.character(anova.table.50.6[6,2]))/as.numeric(as.character(anova.table.50.6[6,3])), digits=5) # And then again dropping Ara + 1 anova.table.50.5<-matrix(nrow=6, ncol=7) anova.table.50.5[1,1]<-"Also exclusing population Ara+1" anova.table.50.5[2,1]<-"Source" anova.table.50.5[2,2]<-"SS" anova.table.50.5[2,3]<-"df" anova.table.50.5[2,4]<-"MS" anova.table.50.5[2,5]<-"F" anova.table.50.5[2,6]<-"p" anova.table.50.5[3,1]<-"Population" anova.table.50.5[3,2]<-round(as.numeric(as.character(summary(anova.5.50)[[1]][[1]][[2]])), digits=5) anova.table.50.5[3,3]<-as.numeric(as.character(summary(anova.5.50)[[1]][[1]][[1]])) anova.table.50.5[3,4]<-round(as.numeric(as.character(summary(anova.5.50)[[1]][[1]][[3]])), digits=5) anova.table.50.5[4,1]<-"Block" anova.table.50.5[4,2]<-round(as.numeric(as.character(summary(anova.5.50)[[2]][[1]][[2]])), digits=5) anova.table.50.5[4,3]<-as.numeric(as.character(summary(anova.5.50)[[2]][[1]][[1]])) anova.table.50.5[4,4]<-round(as.numeric(as.character(summary(anova.5.50)[[2]][[1]][[3]])), digits=5) anova.table.50.5[5,1]<-"Error" anova.table.50.5[5,2]<-round(as.numeric(as.character(summary(anova.5.50)[[3]][[1]][[2]])), digits=5) anova.table.50.5[5,3]<-as.numeric(as.character(summary(anova.5.50)[[3]][[1]][[1]])) anova.table.50.5[5,4]<-round(as.numeric(as.character(summary(anova.5.50)[[3]][[1]][[3]])), digits=5) anova.table.50.5[6,1]<-"Total" anova.table.50.5[3,5]<-round(as.numeric(as.character(summary(anova.5.50)[[1]][[1]][[3]]))/as.numeric(as.character(summary(anova.5.50)[[3]][[1]][[3]])), digits=5) anova.table.50.5[3,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.50.5[3,5])), as.numeric(as.character(summary(anova.5.50)[[1]][[1]][[1]])), as.numeric(as.character(summary(anova.5.50)[[3]][[1]][[1]])) ))), digits=4) anova.table.50.5[3,7]<-F.p.50.5 anova.table.50.5[4,5]<-round(as.numeric(as.character(summary(anova.5.50)[[2]][[1]][[3]]))/as.numeric(as.character(summary(anova.5.50)[[3]][[1]][[3]])), digits=5) anova.table.50.5[4,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.50.5[4,5])), as.numeric(as.character(summary(anova.5.50)[[2]][[1]][[1]])), as.numeric(as.character(summary(anova.5.50)[[3]][[1]][[1]])) ))), digits=4) if(as.numeric(as.character(anova.table.50.5[3,6]))<0.0001) { anova.table.50.5[3,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.50.5[4,6]))<0.0001) { anova.table.50.5[4,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.50.5[3,7]))<0.0001) { anova.table.50.5[3,7]<-"< 0.0001"} anova.table.50.5[6,2]<-as.numeric(as.character(anova.table.50.5[3,2]))+as.numeric(as.character(anova.table.50.5[4,2]))+as.numeric(as.character(anova.table.50.5[5,2])) anova.table.50.5[6,3]<-as.numeric(as.character(anova.table.50.5[3,3]))+as.numeric(as.character(anova.table.50.5[4,3]))+as.numeric(as.character(anova.table.50.5[5,3])) anova.table.50.5[6,4]<-round(as.numeric(as.character(anova.table.50.5[6,2]))/as.numeric(as.character(anova.table.50.5[6,3])), digits=5) # And now we stitch these together to form Table 5: Table.5<-rbind(anova.table.50.full, anova.table.50.6, anova.table.50.5) Table.5<-as.data.frame(Table.5) # Tables were renumbered on submission, as many were moved to the Supplement. What used to be Table 5 is now Table S4 write.csv(file="Table.S4.csv", Table.5) #################### # And now Table 6, doing the same thing for the data from 60,000 generations: anova.table.60.full<-matrix(nrow=6, ncol=7) anova.table.60.full[1,1]<-"All nine populations" anova.table.60.full[2,1]<-"Source" anova.table.60.full[2,2]<-"SS" anova.table.60.full[2,3]<-"df" anova.table.60.full[2,4]<-"MS" anova.table.60.full[2,5]<-"F" anova.table.60.full[2,6]<-"p" anova.table.60.full[2,7]<-"Friedman's p" anova.table.60.full[3,1]<-"Population" anova.table.60.full[3,2]<-round(as.numeric(as.character(summary(anova.full.60)[[1]][[1]][[2]])), digits=5) anova.table.60.full[3,3]<-as.numeric(as.character(summary(anova.full.60)[[1]][[1]][[1]])) anova.table.60.full[3,4]<-round(as.numeric(as.character(summary(anova.full.60)[[1]][[1]][[3]])), digits=5) anova.table.60.full[4,1]<-"Block" anova.table.60.full[4,2]<-round(as.numeric(as.character(summary(anova.full.60)[[2]][[1]][[2]])), digits=5) anova.table.60.full[4,3]<-as.numeric(as.character(summary(anova.full.60)[[2]][[1]][[1]])) anova.table.60.full[4,4]<-round(as.numeric(as.character(summary(anova.full.60)[[2]][[1]][[3]])), digits=5) anova.table.60.full[5,1]<-"Error" anova.table.60.full[5,2]<-round(as.numeric(as.character(summary(anova.full.60)[[3]][[1]][[2]])), digits=5) anova.table.60.full[5,3]<-as.numeric(as.character(summary(anova.full.60)[[3]][[1]][[1]])) anova.table.60.full[5,4]<-round(as.numeric(as.character(summary(anova.full.60)[[3]][[1]][[3]])), digits=5) anova.table.60.full[6,1]<-"Total" anova.table.60.full[3,5]<-round(as.numeric(as.character(summary(anova.full.60)[[1]][[1]][[3]]))/as.numeric(as.character(summary(anova.full.60)[[3]][[1]][[3]])), digits=5) anova.table.60.full[3,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.60.full[3,5])), as.numeric(as.character(summary(anova.full.60)[[1]][[1]][[1]])), as.numeric(as.character(summary(anova.full.60)[[3]][[1]][[1]])) ))), digits=4) anova.table.60.full[3,7]<-F.p.60.full anova.table.60.full[4,5]<-round(as.numeric(as.character(summary(anova.full.60)[[2]][[1]][[3]]))/as.numeric(as.character(summary(anova.full.60)[[3]][[1]][[3]])), digits=5) anova.table.60.full[4,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.60.full[4,5])), as.numeric(as.character(summary(anova.full.60)[[2]][[1]][[1]])), as.numeric(as.character(summary(anova.full.60)[[3]][[1]][[1]])) ))), digits=4) if(as.numeric(as.character(anova.table.60.full[3,6]))<0.0001) { anova.table.60.full[3,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.60.full[4,6]))<0.0001) { anova.table.60.full[4,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.60.full[3,7]))<0.0001) { anova.table.60.full[3,7]<-"< 0.0001"} anova.table.60.full[6,2]<-as.numeric(as.character(anova.table.60.full[3,2]))+as.numeric(as.character(anova.table.60.full[4,2]))+as.numeric(as.character(anova.table.60.full[5,2])) anova.table.60.full[6,3]<-as.numeric(as.character(anova.table.60.full[3,3]))+as.numeric(as.character(anova.table.60.full[4,3]))+as.numeric(as.character(anova.table.60.full[5,3])) anova.table.60.full[6,4]<-round(as.numeric(as.character(anova.table.60.full[6,2]))/as.numeric(as.character(anova.table.60.full[6,3])), digits=5) # Then for just the 6 non-hypermutator populations: anova.table.60.6<-matrix(nrow=6, ncol=7) anova.table.60.6[1,1]<-"Excluding three hypermutator populations" anova.table.60.6[2,1]<-"Source" anova.table.60.6[2,2]<-"SS" anova.table.60.6[2,3]<-"df" anova.table.60.6[2,4]<-"MS" anova.table.60.6[2,5]<-"F" anova.table.60.6[2,6]<-"p" anova.table.60.6[2,7]<-"Friedman's p" anova.table.60.6[3,1]<-"Population" anova.table.60.6[3,2]<-round(as.numeric(as.character(summary(anova.6.60)[[1]][[1]][[2]])), digits=5) anova.table.60.6[3,3]<-as.numeric(as.character(summary(anova.6.60)[[1]][[1]][[1]])) anova.table.60.6[3,4]<-round(as.numeric(as.character(summary(anova.6.60)[[1]][[1]][[3]])), digits=5) anova.table.60.6[4,1]<-"Block" anova.table.60.6[4,2]<-round(as.numeric(as.character(summary(anova.6.60)[[2]][[1]][[2]])), digits=5) anova.table.60.6[4,3]<-as.numeric(as.character(summary(anova.6.60)[[2]][[1]][[1]])) anova.table.60.6[4,4]<-round(as.numeric(as.character(summary(anova.6.60)[[2]][[1]][[3]])), digits=5) anova.table.60.6[5,1]<-"Error" anova.table.60.6[5,2]<-round(as.numeric(as.character(summary(anova.6.60)[[3]][[1]][[2]])), digits=5) anova.table.60.6[5,3]<-as.numeric(as.character(summary(anova.6.60)[[3]][[1]][[1]])) anova.table.60.6[5,4]<-round(as.numeric(as.character(summary(anova.6.60)[[3]][[1]][[3]])), digits=5) anova.table.60.6[6,1]<-"Total" anova.table.60.6[3,5]<-round(as.numeric(as.character(summary(anova.6.60)[[1]][[1]][[3]]))/as.numeric(as.character(summary(anova.6.60)[[3]][[1]][[3]])), digits=5) anova.table.60.6[3,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.60.6[3,5])), as.numeric(as.character(summary(anova.6.60)[[1]][[1]][[1]])), as.numeric(as.character(summary(anova.6.60)[[3]][[1]][[1]])) ))), digits=4) anova.table.60.6[3,7]<-F.p.60.6 anova.table.60.6[4,5]<-round(as.numeric(as.character(summary(anova.6.60)[[2]][[1]][[3]]))/as.numeric(as.character(summary(anova.6.60)[[3]][[1]][[3]])), digits=5) anova.table.60.6[4,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.60.6[4,5])), as.numeric(as.character(summary(anova.6.60)[[2]][[1]][[1]])), as.numeric(as.character(summary(anova.6.60)[[3]][[1]][[1]])) ))), digits=4) if(as.numeric(as.character(anova.table.60.6[3,6]))<0.0001) { anova.table.60.6[3,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.60.6[4,6]))<0.0001) { anova.table.60.6[4,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.60.6[3,7]))<0.0001) { anova.table.60.6[3,7]<-"< 0.0001"} anova.table.60.6[6,2]<-as.numeric(as.character(anova.table.60.6[3,2]))+as.numeric(as.character(anova.table.60.6[4,2]))+as.numeric(as.character(anova.table.60.6[5,2])) anova.table.60.6[6,3]<-as.numeric(as.character(anova.table.60.6[3,3]))+as.numeric(as.character(anova.table.60.6[4,3]))+as.numeric(as.character(anova.table.60.6[5,3])) anova.table.60.6[6,4]<-round(as.numeric(as.character(anova.table.60.6[6,2]))/as.numeric(as.character(anova.table.60.6[6,3])), digits=5) # And then again dropping Ara + 1 anova.table.60.5<-matrix(nrow=6, ncol=7) anova.table.60.5[1,1]<-"Also exclusing population Ara+1" anova.table.60.5[2,1]<-"Source" anova.table.60.5[2,2]<-"SS" anova.table.60.5[2,3]<-"df" anova.table.60.5[2,4]<-"MS" anova.table.60.5[2,5]<-"F" anova.table.60.5[2,6]<-"p" anova.table.60.5[2,7]<-"Friedman's p" anova.table.60.5[3,1]<-"Population" anova.table.60.5[3,2]<-round(as.numeric(as.character(summary(anova.5.60)[[1]][[1]][[2]])), digits=5) anova.table.60.5[3,3]<-as.numeric(as.character(summary(anova.5.60)[[1]][[1]][[1]])) anova.table.60.5[3,4]<-round(as.numeric(as.character(summary(anova.5.60)[[1]][[1]][[3]])), digits=5) anova.table.60.5[4,1]<-"Block" anova.table.60.5[4,2]<-round(as.numeric(as.character(summary(anova.5.60)[[2]][[1]][[2]])), digits=5) anova.table.60.5[4,3]<-as.numeric(as.character(summary(anova.5.60)[[2]][[1]][[1]])) anova.table.60.5[4,4]<-round(as.numeric(as.character(summary(anova.5.60)[[2]][[1]][[3]])), digits=5) anova.table.60.5[5,1]<-"Error" anova.table.60.5[5,2]<-round(as.numeric(as.character(summary(anova.5.60)[[3]][[1]][[2]])), digits=5) anova.table.60.5[5,3]<-as.numeric(as.character(summary(anova.5.60)[[3]][[1]][[1]])) anova.table.60.5[5,4]<-round(as.numeric(as.character(summary(anova.5.60)[[3]][[1]][[3]])), digits=5) anova.table.60.5[6,1]<-"Total" anova.table.60.5[3,5]<-round(as.numeric(as.character(summary(anova.5.60)[[1]][[1]][[3]]))/as.numeric(as.character(summary(anova.5.60)[[3]][[1]][[3]])), digits=5) anova.table.60.5[3,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.60.5[3,5])), as.numeric(as.character(summary(anova.5.60)[[1]][[1]][[1]])), as.numeric(as.character(summary(anova.5.60)[[3]][[1]][[1]])) ))), digits=4) anova.table.60.5[3,7]<-F.p.50.5 anova.table.60.5[4,5]<-round(as.numeric(as.character(summary(anova.5.60)[[2]][[1]][[3]]))/as.numeric(as.character(summary(anova.5.60)[[3]][[1]][[3]])), digits=5) anova.table.60.5[4,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.60.5[4,5])), as.numeric(as.character(summary(anova.5.60)[[2]][[1]][[1]])), as.numeric(as.character(summary(anova.5.60)[[3]][[1]][[1]])) ))), digits=4) if(as.numeric(as.character(anova.table.60.5[3,6]))<0.0001) { anova.table.60.5[3,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.60.5[4,6]))<0.0001) { anova.table.60.5[4,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.60.5[3,7]))<0.0001) { anova.table.60.5[3,7]<-"< 0.0001"} anova.table.60.5[6,2]<-as.numeric(as.character(anova.table.60.5[3,2]))+as.numeric(as.character(anova.table.60.5[4,2]))+as.numeric(as.character(anova.table.60.5[5,2])) anova.table.60.5[6,3]<-as.numeric(as.character(anova.table.60.5[3,3]))+as.numeric(as.character(anova.table.60.5[4,3]))+as.numeric(as.character(anova.table.60.5[5,3])) anova.table.60.5[6,4]<-round(as.numeric(as.character(anova.table.60.5[6,2]))/as.numeric(as.character(anova.table.60.5[6,3])), digits=5) # And now we stitch these together to form Table 6: Table.6<-rbind(anova.table.60.full, anova.table.60.6, anova.table.60.5) Table.6<-as.data.frame(Table.6) # Tables were numbered on submission. What used to be Table 6 is now Table S5. write.csv(file="Table.S5.csv", Table.6) # Now to create Table 7. We have already separated the data into 3 subsets, one for each # generation. Those are data.40k, data.50k, and data.60k. # Now we create a linear model of each time point. lm.40k.9<-lm(Fitness ~ Evolved.Population, data=data.40k, na.action=na.omit, model=T) lm.50k.9<-lm(Fitness ~ Evolved.Population, data=data.50k, na.action=na.omit, model=T) lm.60k.9<-lm(Fitness ~ Evolved.Population, data=data.60k, na.action=na.omit, model=T) # We do an ANOVA on each linear model: a.40k.9<-anova(lm.40k.9) a.50k.9<-anova(lm.50k.9) a.60k.9<-anova(lm.60k.9) # We find the difference in teh SS Pop and the SS Resid: d.40k.9<-a.40k.9[1,3] - a.40k.9[2,3] d.50k.9<-a.50k.9[1,3] - a.50k.9[2,3] d.60k.9<-a.60k.9[1,3] - a.60k.9[2,3] # We divide the difference by the number of replicates: d.40k.9<-d.40k.9/length(data.40k$Fitness) d.50k.9<-d.50k.9/length(data.50k$Fitness) d.60k.9<-d.60k.9/length(data.60k$Fitness) # Get the signed square root of these differences: differences.9<-c(d.40k.9, d.50k.9, d.60k.9) Signed.Sq.Roots.9<-vector("list", length(differences.9)) for (i in 1:length(differences.9)) { value<-differences.9[i] if (value > 0 ) { signed.sq.root<-value^0.5 } if (value < 0 ) { abs.value<-abs(value) abs.sq.root<-abs.value^0.5 signed.sq.root<-(-1 * abs.sq.root) } Signed.Sq.Roots.9[i]<-signed.sq.root } Signed.Sq.Roots.9<-as.numeric(as.character(Signed.Sq.Roots.9)) # Now we repeat these procedures, excluding the hypermutator populations: data.40k.6<-subset(data.40k, Evolved.Population=="Ara - 5" | Evolved.Population=="Ara - 6" | Evolved.Population=="Ara + 1" | Evolved.Population=="Ara + 2" | Evolved.Population=="Ara + 4" | Evolved.Population=="Ara + 5") data.50k.6<-subset(data.50k, Evolved.Population=="Ara - 5" | Evolved.Population=="Ara - 6" | Evolved.Population=="Ara + 1" | Evolved.Population=="Ara + 2" | Evolved.Population=="Ara + 4" | Evolved.Population=="Ara + 5") data.60k.6<-subset(data.60k, Evolved.Population=="Ara - 5" | Evolved.Population=="Ara - 6" | Evolved.Population=="Ara + 1" | Evolved.Population=="Ara + 2" | Evolved.Population=="Ara + 4" | Evolved.Population=="Ara + 5") lm.40k.6<-lm(Fitness ~ Evolved.Population, data=data.40k.6, na.action=na.omit, model=T) lm.50k.6<-lm(Fitness ~ Evolved.Population, data=data.50k.6, na.action=na.omit, model=T) lm.60k.6<-lm(Fitness ~ Evolved.Population, data=data.60k.6, na.action=na.omit, model=T) # We do an ANOVA on each linear model: a.40k.6<-anova(lm.40k.6) a.50k.6<-anova(lm.50k.6) a.60k.6<-anova(lm.60k.6) # We find the difference in teh SS Pop and the SS Resid: d.40k.6<-a.40k.6[1,3] - a.40k.6[2,3] d.50k.6<-a.50k.6[1,3] - a.50k.6[2,3] d.60k.6<-a.60k.6[1,3] - a.60k.6[2,3] # We divide the difference by the number of replicates: d.40k.6<-d.40k.6/length(data.40k.6$Fitness) d.50k.6<-d.50k.6/length(data.50k.6$Fitness) d.60k.6<-d.60k.6/length(data.60k.6$Fitness) # Get the signed square root of these differences: differences.6<-c(d.40k.6, d.50k.6, d.60k.6) Signed.Sq.Roots.6<-vector("list", length(differences.6)) for (i in 1:length(differences.6)) { value<-differences.6[i] if (value > 0 ) { signed.sq.root<-value^0.5 } if (value < 0 ) { abs.value<-abs(value) abs.sq.root<-abs.value^0.5 signed.sq.root<-(-1 * abs.sq.root) } Signed.Sq.Roots.6[i]<-signed.sq.root } Signed.Sq.Roots.6<-as.numeric(as.character(Signed.Sq.Roots.6)) # And now again, excluding population Ara + 1: data.40k.5<-subset(data.40k, Evolved.Population=="Ara - 5" | Evolved.Population=="Ara - 6" | Evolved.Population=="Ara + 2" | Evolved.Population=="Ara + 4" | Evolved.Population=="Ara + 5") data.50k.5<-subset(data.50k, Evolved.Population=="Ara - 5" | Evolved.Population=="Ara - 6" | Evolved.Population=="Ara + 2" | Evolved.Population=="Ara + 4" | Evolved.Population=="Ara + 5") data.60k.5<-subset(data.60k, Evolved.Population=="Ara - 5" | Evolved.Population=="Ara - 6" | Evolved.Population=="Ara + 2" | Evolved.Population=="Ara + 4" | Evolved.Population=="Ara + 5") lm.40k.5<-lm(Fitness ~ Evolved.Population, data=data.40k.5, na.action=na.omit, model=T) lm.50k.5<-lm(Fitness ~ Evolved.Population, data=data.50k.5, na.action=na.omit, model=T) lm.60k.5<-lm(Fitness ~ Evolved.Population, data=data.60k.5, na.action=na.omit, model=T) # We do an ANOVA on each linear model: a.40k.5<-anova(lm.40k.5) a.50k.5<-anova(lm.50k.5) a.60k.5<-anova(lm.60k.5) # We find the difference in teh SS Pop and the SS Resid: d.40k.5<-a.40k.5[1,3] - a.40k.5[2,3] d.50k.5<-a.50k.5[1,3] - a.50k.5[2,3] d.60k.5<-a.60k.5[1,3] - a.60k.5[2,3] # We divide the difference by the number of replicates: d.40k.5<-d.40k.5/length(data.40k.5$Fitness) d.50k.5<-d.50k.5/length(data.50k.5$Fitness) d.60k.5<-d.60k.5/length(data.60k.5$Fitness) # Get the signed square root of these differences: differences.5<-c(d.40k.5, d.50k.5, d.60k.5) Signed.Sq.Roots.5<-vector("list", length(differences.5)) for (i in 1:length(differences.5)) { value<-differences.5[i] if (value > 0 ) { signed.sq.root<-value^0.5 } if (value < 0 ) { abs.value<-abs(value) abs.sq.root<-abs.value^0.5 signed.sq.root<-(-1 * abs.sq.root) } Signed.Sq.Roots.5[i]<-signed.sq.root } Signed.Sq.Roots.5<-as.numeric(as.character(Signed.Sq.Roots.5)) # Now we organize this information into what we want as Table 7: Table.7<-matrix(nrow=14, ncol=3) Table.7[1, 1]<-"At 40,000 generations" Table.7[2, 1]<-"All nine populations" Table.7[3, 1]<-"Excluding three hypermutators" Table.7[4, 1]<-"Also excluding population Ara+1" Table.7[6, 1]<-"At 40,000 generations" Table.7[7, 1]<-"All nine populations" Table.7[8, 1]<-"Excluding three hypermutators" Table.7[9, 1]<-"Also excluding population Ara+1" Table.7[11, 1]<-"At 40,000 generations" Table.7[12, 1]<-"All nine populations" Table.7[13, 1]<-"Excluding three hypermutators" Table.7[14, 1]<-"Also excluding population Ara+1" Table.7[1, 2]<-"Vpop" Table.7[2, 2]<-round(((as.numeric(as.character(anova.table.40.full[3,4]))-as.numeric(as.character(anova.table.40.full[5,4])))/(as.numeric(as.character(anova.table.40.full[4,3]))+1)), digits=6) Table.7[3, 2]<-round(((as.numeric(as.character(anova.table.40.6[3,4]))-as.numeric(as.character(anova.table.40.6[5,4])))/(as.numeric(as.character(anova.table.40.6[4,3]))+1)), digits=6) Table.7[4, 2]<-round(((as.numeric(as.character(anova.table.40.5[3,4]))-as.numeric(as.character(anova.table.40.5[5,4])))/(as.numeric(as.character(anova.table.40.5[4,3]))+1)), digits=6) Table.7[6, 2]<-"Vpop" Table.7[7, 2]<-round(((as.numeric(as.character(anova.table.50.full[3,4]))-as.numeric(as.character(anova.table.50.full[5,4])))/(as.numeric(as.character(anova.table.50.full[4,3]))+1)), digits=6) Table.7[8, 2]<-round(((as.numeric(as.character(anova.table.50.6[3,4]))-as.numeric(as.character(anova.table.50.6[5,4])))/(as.numeric(as.character(anova.table.50.6[4,3]))+1)), digits=6) Table.7[9, 2]<-round(((as.numeric(as.character(anova.table.50.5[3,4]))-as.numeric(as.character(anova.table.50.5[5,4])))/(as.numeric(as.character(anova.table.50.5[4,3]))+1)), digits=6) Table.7[11, 2]<-"Vpop" Table.7[12, 2]<-round(((as.numeric(as.character(anova.table.60.full[3,4]))-as.numeric(as.character(anova.table.60.full[5,4])))/(as.numeric(as.character(anova.table.60.full[4,3]))+1)), digits=6) Table.7[13, 2]<-round(((as.numeric(as.character(anova.table.60.6[3,4]))-as.numeric(as.character(anova.table.60.6[5,4])))/(as.numeric(as.character(anova.table.60.6[4,3]))+1)), digits=6) Table.7[14, 2]<-round(((as.numeric(as.character(anova.table.60.5[3,4]))-as.numeric(as.character(anova.table.60.5[5,4])))/(as.numeric(as.character(anova.table.60.5[4,3]))+1)), digits=6) Table.7[1, 3]<-"sigmapop" Table.7[2, 3]<-round(((as.numeric(as.character(anova.table.40.full[3,4]))-as.numeric(as.character(anova.table.40.full[5,4])))/(as.numeric(as.character(anova.table.40.full[4,3]))+1))^.5, digits=4) Table.7[3, 3]<-round(((as.numeric(as.character(anova.table.40.6[3,4]))-as.numeric(as.character(anova.table.40.6[5,4])))/(as.numeric(as.character(anova.table.40.6[4,3]))+1))^.5, digits=4) Table.7[4, 3]<-round(((as.numeric(as.character(anova.table.40.5[3,4]))-as.numeric(as.character(anova.table.40.5[5,4])))/(as.numeric(as.character(anova.table.40.5[4,3]))+1))^.5, digits=4) Table.7[6, 3]<-"sigmapop" Table.7[7, 3]<-round(((as.numeric(as.character(anova.table.50.full[3,4]))-as.numeric(as.character(anova.table.50.full[5,4])))/(as.numeric(as.character(anova.table.50.full[4,3]))+1))^.5, digits=4) Table.7[8, 3]<-round(((as.numeric(as.character(anova.table.50.6[3,4]))-as.numeric(as.character(anova.table.50.6[5,4])))/(as.numeric(as.character(anova.table.50.6[4,3]))+1))^.5, digits=4) Table.7[9, 3]<-round(((as.numeric(as.character(anova.table.50.5[3,4]))-as.numeric(as.character(anova.table.50.5[5,4])))/(as.numeric(as.character(anova.table.50.5[4,3]))+1))^.5, digits=4) Table.7[11, 3]<-"sigmapop" Table.7[12, 3]<-round(((as.numeric(as.character(anova.table.60.full[3,4]))-as.numeric(as.character(anova.table.60.full[5,4])))/(as.numeric(as.character(anova.table.60.full[4,3]))+1))^.5, digits=4) Table.7[13, 3]<-round(((as.numeric(as.character(anova.table.60.6[3,4]))-as.numeric(as.character(anova.table.60.6[5,4])))/(as.numeric(as.character(anova.table.60.6[4,3]))+1))^.5, digits=4) Table.7[14, 3]<-round(((as.numeric(as.character(anova.table.60.5[3,4]))-as.numeric(as.character(anova.table.60.5[5,4])))/(as.numeric(as.character(anova.table.60.5[4,3]))+1))^.5, digits=4) # Tables were renumbered on submission, as many were moved to the Supplement. What used to be Table 7 is now Table 3. write.csv(file="Table.3.csv", Table.7) ######## # Now, for Tables 8 and 9, we want to find just the blocks which are complete. # These will be the cases where the particular variable for the Dummy.Replicate # column appears 27 times. data.complete<-{} possible.replicates<-unique(data.accepted$Dummy.Replicate) for (i in 1:length(possible.replicates) ) { current.rep<-possible.replicates[i] current.data<-subset(data.accepted, Dummy.Replicate==current.rep) rep.count<-length(current.data$Fitness) if(rep.count == 27) { data.complete<-rbind(data.complete, current.data) } } # This gives us 29 blocks of complete data. # Within this complete data, we need to figure out how much fitness changed within # each replicate, from 40,000 to 50,000 generations; from 50,000 to 60,000 generations; # and from 40,000 to 60,000 generations. This will involve restructing the data file. complete.replicates<-unique(data.complete$Dummy.Replicate) evolved.list<-unique(data.complete$Evolved.Population) #numeric.block<-seq(from=1, to=length(complete.replicates, by=1) data.by.change<-matrix(nrow=length(complete.replicates)*length(evolved.list), ncol=5) for (i in 1:length(complete.replicates)) { current.replicate<-complete.replicates[i] for (j in 1:length(evolved.list)) { current.pop<-evolved.list[j] current.data<-subset(data.complete, Dummy.Replicate==current.replicate) current.data<-subset(current.data, Evolved.Population==current.pop) current.40<-subset(current.data, Generation==40000) current.50<-subset(current.data, Generation==50000) current.60<-subset(current.data, Generation==60000) value.40<-current.40$Fitness value.50<-current.50$Fitness value.60<-current.60$Fitness diff.full<-value.60-value.40 diff.first<-value.50-value.40 diff.second<-value.60-value.50 data.by.change[((i-1)*length(evolved.list))+j, 1]<-as.character(current.pop) data.by.change[((i-1)*length(evolved.list))+j, 2]<-as.character(current.replicate) data.by.change[((i-1)*length(evolved.list))+j, 3]<-diff.full data.by.change[((i-1)*length(evolved.list))+j, 4]<-diff.first data.by.change[((i-1)*length(evolved.list))+j, 5]<-diff.second } } colnames(data.by.change)<-c("Evolved.Population", "Dummy.Replicate", "Full.Change", "Early.Change", "Late.Change") data.by.change<-as.data.frame(data.by.change) data.by.change$Full.Change<-as.numeric(as.character(data.by.change$Full.Change)) data.by.change$Early.Change<-as.numeric(as.character(data.by.change$Early.Change)) data.by.change$Late.Change<-as.numeric(as.character(data.by.change$Late.Change)) data.by.change$Deceleration<-data.by.change$Early.Change-data.by.change$Late.Change data.by.change.9<-data.by.change data.by.change.6<-subset(data.by.change, Evolved.Population=="Ara + 1" | Evolved.Population=="Ara + 2" | Evolved.Population=="Ara + 4" | Evolved.Population=="Ara + 5" | Evolved.Population=="Ara - 5" | Evolved.Population=="Ara - 6") # Now we can perform the necessary ANOVAs: anova.change.full.9<-aov(Full.Change ~ Error(Evolved.Population*Dummy.Replicate), data=data.by.change.9) anova.change.full.6<-aov(Full.Change ~ Error(Evolved.Population*Dummy.Replicate), data=data.by.change.6) ###################### # Now I need to rework the Friedman p function to handle a different type of data. Friedman.p.2<-function(data.set) { Rep.list<-unique(data.set$Dummy.Replicate) Pop.list<-unique(data.set$Evolved.Population) rank.matrix<-matrix(nrow=length(Pop.list), ncol=length(Rep.list)) for (i in 1:length(Rep.list)) { current.rep<-Rep.list[i] current.data<-subset(data.set, Dummy.Replicate==current.rep) Evolved.list<-current.data$Evolved.Population Change.list<-current.data$Full.Change current.matrix<-matrix(ncol=3, nrow=length(Pop.list)) current.matrix[,1]<-(as.character(Evolved.list)) current.matrix[,2]<-as.numeric(as.character(Change.list)) current.matrix[,2]<-as.numeric(as.character(current.matrix[,2]))+1000 current.matrix<-current.matrix[order(as.numeric(as.character(current.matrix[,2]))),] current.matrix[,3]<-seq(1:length(Pop.list)) current.matrix<-current.matrix[order(current.matrix[,1]),] current.matrix[,2]<-as.numeric(as.character(current.matrix[,2]))-1000 rank.matrix[,i]<-as.numeric(as.character(current.matrix[,3])) } rank.sum<-matrix(nrow=length(Pop.list), ncol=1) for (i in 1:length(Pop.list)) { current.data<-as.numeric(as.character(rank.matrix[i,])) current.sum<-sum(as.numeric(as.character(current.data))) rank.sum[i, 1]<-current.sum } a<-length(Pop.list) b<-length(Rep.list) squared.sum.ranks<-0 for (i in 1:length(Pop.list)) { current.rank.sum<-as.numeric(as.character(rank.sum[i,1])) value<-as.numeric(as.character(current.rank.sum))^2 squared.sum.ranks<-as.numeric(as.character(squared.sum.ranks))+as.numeric(as.character(value)) } computed.value<-((12/(a*b*(a+1)))*(as.numeric(as.character(squared.sum.ranks)))) - (3*b*(a+1)) calc.p<-1-pchisq(computed.value, (a-1)) return(calc.p) } ### F.p.change.9<-Friedman.p.2(data.by.change.9) F.p.change.6<-Friedman.p.2(data.by.change.6) # Now we pull out the right information from these ANOVAs, just as we did for Tables 4-6. # First for all the populations: anova.table.full.9<-matrix(nrow=6, ncol=7) anova.table.full.9[1,1]<-"All nine populations" anova.table.full.9[2,1]<-"Source" anova.table.full.9[2,2]<-"SS" anova.table.full.9[2,3]<-"df" anova.table.full.9[2,4]<-"MS" anova.table.full.9[2,5]<-"F" anova.table.full.9[2,6]<-"p" anova.table.full.9[2,7]<-"Friedman's p" anova.table.full.9[3,1]<-"Population" anova.table.full.9[3,2]<-round(as.numeric(as.character(summary(anova.change.full.9)[[1]][[1]][[2]])), digits=5) anova.table.full.9[3,3]<-as.numeric(as.character(summary(anova.change.full.9)[[1]][[1]][[1]])) anova.table.full.9[3,4]<-round(as.numeric(as.character(summary(anova.change.full.9)[[1]][[1]][[3]])), digits=5) anova.table.full.9[4,1]<-"Block" anova.table.full.9[4,2]<-round(as.numeric(as.character(summary(anova.change.full.9)[[2]][[1]][[2]])), digits=5) anova.table.full.9[4,3]<-as.numeric(as.character(summary(anova.change.full.9)[[2]][[1]][[1]])) anova.table.full.9[4,4]<-round(as.numeric(as.character(summary(anova.change.full.9)[[2]][[1]][[3]])), digits=5) anova.table.full.9[5,1]<-"Error" anova.table.full.9[5,2]<-round(as.numeric(as.character(summary(anova.change.full.9)[[3]][[1]][[2]])), digits=5) anova.table.full.9[5,3]<-as.numeric(as.character(summary(anova.change.full.9)[[3]][[1]][[1]])) anova.table.full.9[5,4]<-round(as.numeric(as.character(summary(anova.change.full.9)[[3]][[1]][[3]])), digits=5) anova.table.full.9[6,1]<-"Total" anova.table.full.9[3,5]<-round(as.numeric(as.character(summary(anova.change.full.9)[[1]][[1]][[3]]))/as.numeric(as.character(summary(anova.change.full.9)[[3]][[1]][[3]])), digits=5) anova.table.full.9[3,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.full.9[3,5])), as.numeric(as.character(summary(anova.change.full.9)[[1]][[1]][[1]])), as.numeric(as.character(summary(anova.change.full.9)[[3]][[1]][[1]])) ))), digits=4) anova.table.full.9[3,7]<-round(F.p.change.9, digits=5) anova.table.full.9[4,5]<-round(as.numeric(as.character(summary(anova.change.full.9)[[2]][[1]][[3]]))/as.numeric(as.character(summary(anova.change.full.9)[[3]][[1]][[3]])), digits=5) anova.table.full.9[4,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.full.9[4,5])), as.numeric(as.character(summary(anova.change.full.9)[[2]][[1]][[1]])), as.numeric(as.character(summary(anova.change.full.9)[[3]][[1]][[1]])) ))), digits=4) if(as.numeric(as.character(anova.table.full.9[3,6]))<0.0001) { anova.table.full.9[3,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.full.9[4,6]))<0.0001) { anova.table.full.9[4,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.full.9[3,7]))<0.0001) { anova.table.full.9[3,7]<-"< 0.0001"} anova.table.full.9[6,2]<-as.numeric(as.character(anova.table.full.9[3,2]))+as.numeric(as.character(anova.table.full.9[4,2]))+as.numeric(as.character(anova.table.full.9[5,2])) anova.table.full.9[6,3]<-as.numeric(as.character(anova.table.full.9[3,3]))+as.numeric(as.character(anova.table.full.9[4,3]))+as.numeric(as.character(anova.table.full.9[5,3])) anova.table.full.9[6,4]<-round(as.numeric(as.character(anova.table.full.9[6,2]))/as.numeric(as.character(anova.table.full.9[6,3])), digits=5) # Then for just the 6 non hypermutators: anova.table.full.6<-matrix(nrow=6, ncol=7) anova.table.full.6[1,1]<-"Excluding three hypermutator populations" anova.table.full.6[2,1]<-"Source" anova.table.full.6[2,2]<-"SS" anova.table.full.6[2,3]<-"df" anova.table.full.6[2,4]<-"MS" anova.table.full.6[2,5]<-"F" anova.table.full.6[2,6]<-"p" anova.table.full.6[2,7]<-"Friedman's p" anova.table.full.6[3,1]<-"Population" anova.table.full.6[3,2]<-round(as.numeric(as.character(summary(anova.change.full.6)[[1]][[1]][[2]])), digits=5) anova.table.full.6[3,3]<-as.numeric(as.character(summary(anova.change.full.6)[[1]][[1]][[1]])) anova.table.full.6[3,4]<-round(as.numeric(as.character(summary(anova.change.full.6)[[1]][[1]][[3]])), digits=5) anova.table.full.6[4,1]<-"Block" anova.table.full.6[4,2]<-round(as.numeric(as.character(summary(anova.change.full.6)[[2]][[1]][[2]])), digits=5) anova.table.full.6[4,3]<-as.numeric(as.character(summary(anova.change.full.6)[[2]][[1]][[1]])) anova.table.full.6[4,4]<-round(as.numeric(as.character(summary(anova.change.full.6)[[2]][[1]][[3]])), digits=5) anova.table.full.6[5,1]<-"Error" anova.table.full.6[5,2]<-round(as.numeric(as.character(summary(anova.change.full.6)[[3]][[1]][[2]])), digits=5) anova.table.full.6[5,3]<-as.numeric(as.character(summary(anova.change.full.6)[[3]][[1]][[1]])) anova.table.full.6[5,4]<-round(as.numeric(as.character(summary(anova.change.full.6)[[3]][[1]][[3]])), digits=5) anova.table.full.6[6,1]<-"Total" anova.table.full.6[3,5]<-round(as.numeric(as.character(summary(anova.change.full.6)[[1]][[1]][[3]]))/as.numeric(as.character(summary(anova.change.full.6)[[3]][[1]][[3]])), digits=5) anova.table.full.6[3,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.full.6[3,5])), as.numeric(as.character(summary(anova.change.full.6)[[1]][[1]][[1]])), as.numeric(as.character(summary(anova.change.full.6)[[3]][[1]][[1]])) ))), digits=4) anova.table.full.6[3,7]<-round(F.p.change.6, digits=5) anova.table.full.6[4,5]<-round(as.numeric(as.character(summary(anova.change.full.6)[[2]][[1]][[3]]))/as.numeric(as.character(summary(anova.change.full.6)[[3]][[1]][[3]])), digits=5) anova.table.full.6[4,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.full.6[4,5])), as.numeric(as.character(summary(anova.change.full.6)[[2]][[1]][[1]])), as.numeric(as.character(summary(anova.change.full.6)[[3]][[1]][[1]])) ))), digits=4) if(as.numeric(as.character(anova.table.full.6[3,6]))<0.0001) { anova.table.full.6[3,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.full.6[4,6]))<0.0001) { anova.table.full.6[4,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.full.6[3,7]))<0.0001) { anova.table.full.6[3,7]<-"< 0.0001"} anova.table.full.6[6,2]<-as.numeric(as.character(anova.table.full.6[3,2]))+as.numeric(as.character(anova.table.full.6[4,2]))+as.numeric(as.character(anova.table.full.6[5,2])) anova.table.full.6[6,3]<-as.numeric(as.character(anova.table.full.6[3,3]))+as.numeric(as.character(anova.table.full.6[4,3]))+as.numeric(as.character(anova.table.full.6[5,3])) anova.table.full.6[6,4]<-round(as.numeric(as.character(anova.table.full.6[6,2]))/as.numeric(as.character(anova.table.full.6[6,3])), digits=5) Table.8<-rbind(anova.table.full.9, anova.table.full.6) Table.8<-as.data.frame(Table.8) # Tables were renumbered on submission, as many were moved to the Supplement. What used to be Table 8 is now Table S6. write.csv(file="Table.S6.csv", Table.8) # Now for Table 9. It's a lot like Table 8, but based on the Deceleration column, not the Full.Change # That does require a slight modification to the Friedman p calculation function. Friedman.p.3<-function(data.set) { Rep.list<-unique(data.set$Dummy.Replicate) Pop.list<-unique(data.set$Evolved.Population) rank.matrix<-matrix(nrow=length(Pop.list), ncol=length(Rep.list)) for (i in 1:length(Rep.list)) { current.rep<-Rep.list[i] current.data<-subset(data.set, Dummy.Replicate==current.rep) Evolved.list<-current.data$Evolved.Population Deceleration.list<-current.data$Deceleration current.matrix<-matrix(ncol=3, nrow=length(Pop.list)) current.matrix[,1]<-(as.character(Evolved.list)) current.matrix[,2]<-as.numeric(as.character(Deceleration.list)) current.matrix[,2]<-as.numeric(as.character(current.matrix[,2]))+1000 current.matrix<-current.matrix[order(as.numeric(as.character(current.matrix[,2]))),] current.matrix[,3]<-seq(1:length(Pop.list)) current.matrix<-current.matrix[order(current.matrix[,1]),] # current.matrix[,2]<-as.numeric(as.character(current.matrix[,2]))-1000 rank.matrix[,i]<-as.numeric(as.character(current.matrix[,3])) } rank.sum<-matrix(nrow=length(Pop.list), ncol=1) for (i in 1:length(Pop.list)) { current.data<-as.numeric(as.character(rank.matrix[i,])) current.sum<-sum(as.numeric(as.character(current.data))) rank.sum[i, 1]<-current.sum } a<-length(Pop.list) b<-length(Rep.list) squared.sum.ranks<-0 for (i in 1:length(Pop.list)) { current.rank.sum<-as.numeric(as.character(rank.sum[i,1])) value<-as.numeric(as.character(current.rank.sum))^2 squared.sum.ranks<-as.numeric(as.character(squared.sum.ranks))+as.numeric(as.character(value)) } computed.value<-((12/(a*b*(a+1)))*(as.numeric(as.character(squared.sum.ranks)))) - (3*b*(a+1)) calc.p<-1-pchisq(computed.value, (a-1)) return(calc.p) } F.p.deceleration.9<-Friedman.p.3(data.by.change.9) F.p.deceleration.6<-Friedman.p.3(data.by.change.6) anova.change.deceleration.9<-aov(Deceleration ~ Error(Evolved.Population*Dummy.Replicate), data=data.by.change.9) anova.change.deceleration.6<-aov(Deceleration ~ Error(Evolved.Population*Dummy.Replicate), data=data.by.change.6) # Now we pull out the right information from these ANOVAs, just as we did for Tables 4-6. # First for all the populations: anova.table.deceleration.9<-matrix(nrow=6, ncol=7) anova.table.deceleration.9[1,1]<-"All nine populations" anova.table.deceleration.9[2,1]<-"Source" anova.table.deceleration.9[2,2]<-"SS" anova.table.deceleration.9[2,3]<-"df" anova.table.deceleration.9[2,4]<-"MS" anova.table.deceleration.9[2,5]<-"F" anova.table.deceleration.9[2,6]<-"p" anova.table.deceleration.9[2,7]<-"Friedman's p" anova.table.deceleration.9[3,1]<-"Population" anova.table.deceleration.9[3,2]<-round(as.numeric(as.character(summary(anova.change.deceleration.9)[[1]][[1]][[2]])), digits=5) anova.table.deceleration.9[3,3]<-as.numeric(as.character(summary(anova.change.deceleration.9)[[1]][[1]][[1]])) anova.table.deceleration.9[3,4]<-round(as.numeric(as.character(summary(anova.change.deceleration.9)[[1]][[1]][[3]])), digits=5) anova.table.deceleration.9[4,1]<-"Block" anova.table.deceleration.9[4,2]<-round(as.numeric(as.character(summary(anova.change.deceleration.9)[[2]][[1]][[2]])), digits=5) anova.table.deceleration.9[4,3]<-as.numeric(as.character(summary(anova.change.deceleration.9)[[2]][[1]][[1]])) anova.table.deceleration.9[4,4]<-round(as.numeric(as.character(summary(anova.change.deceleration.9)[[2]][[1]][[3]])), digits=5) anova.table.deceleration.9[5,1]<-"Error" anova.table.deceleration.9[5,2]<-round(as.numeric(as.character(summary(anova.change.deceleration.9)[[3]][[1]][[2]])), digits=5) anova.table.deceleration.9[5,3]<-as.numeric(as.character(summary(anova.change.deceleration.9)[[3]][[1]][[1]])) anova.table.deceleration.9[5,4]<-round(as.numeric(as.character(summary(anova.change.deceleration.9)[[3]][[1]][[3]])), digits=5) anova.table.deceleration.9[6,1]<-"Total" anova.table.deceleration.9[3,5]<-round(as.numeric(as.character(summary(anova.change.deceleration.9)[[1]][[1]][[3]]))/as.numeric(as.character(summary(anova.change.deceleration.9)[[3]][[1]][[3]])), digits=5) anova.table.deceleration.9[3,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.deceleration.9[3,5])), as.numeric(as.character(summary(anova.change.deceleration.9)[[1]][[1]][[1]])), as.numeric(as.character(summary(anova.change.deceleration.9)[[3]][[1]][[1]])) ))), digits=4) anova.table.deceleration.9[3,7]<-round(F.p.deceleration.9, digits=5) anova.table.deceleration.9[4,5]<-round(as.numeric(as.character(summary(anova.change.deceleration.9)[[2]][[1]][[3]]))/as.numeric(as.character(summary(anova.change.deceleration.9)[[3]][[1]][[3]])), digits=5) anova.table.deceleration.9[4,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.deceleration.9[4,5])), as.numeric(as.character(summary(anova.change.deceleration.9)[[2]][[1]][[1]])), as.numeric(as.character(summary(anova.change.deceleration.9)[[3]][[1]][[1]])) ))), digits=4) if(as.numeric(as.character(anova.table.deceleration.9[3,6]))<0.0001) { anova.table.deceleration.9[3,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.deceleration.9[4,6]))<0.0001) { anova.table.deceleration.9[4,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.deceleration.9[3,7]))<0.0001) { anova.table.deceleration.9[3,7]<-"< 0.0001"} anova.table.deceleration.9[6,2]<-as.numeric(as.character(anova.table.deceleration.9[3,2]))+as.numeric(as.character(anova.table.deceleration.9[4,2]))+as.numeric(as.character(anova.table.deceleration.9[5,2])) anova.table.deceleration.9[6,3]<-as.numeric(as.character(anova.table.deceleration.9[3,3]))+as.numeric(as.character(anova.table.deceleration.9[4,3]))+as.numeric(as.character(anova.table.deceleration.9[5,3])) anova.table.deceleration.9[6,4]<-round(as.numeric(as.character(anova.table.deceleration.9[6,2]))/as.numeric(as.character(anova.table.deceleration.9[6,3])), digits=5) # Then for just the 6 non hypermutators: anova.table.deceleration.6<-matrix(nrow=6, ncol=7) anova.table.deceleration.6[1,1]<-"Excluding three hypermutator populations" anova.table.deceleration.6[2,1]<-"Source" anova.table.deceleration.6[2,2]<-"SS" anova.table.deceleration.6[2,3]<-"df" anova.table.deceleration.6[2,4]<-"MS" anova.table.deceleration.6[2,5]<-"F" anova.table.deceleration.6[2,6]<-"p" anova.table.deceleration.6[2,7]<-"Friedman's p" anova.table.deceleration.6[3,1]<-"Population" anova.table.deceleration.6[3,2]<-round(as.numeric(as.character(summary(anova.change.deceleration.6)[[1]][[1]][[2]])), digits=5) anova.table.deceleration.6[3,3]<-as.numeric(as.character(summary(anova.change.deceleration.6)[[1]][[1]][[1]])) anova.table.deceleration.6[3,4]<-round(as.numeric(as.character(summary(anova.change.deceleration.6)[[1]][[1]][[3]])), digits=5) anova.table.deceleration.6[4,1]<-"Block" anova.table.deceleration.6[4,2]<-round(as.numeric(as.character(summary(anova.change.deceleration.6)[[2]][[1]][[2]])), digits=5) anova.table.deceleration.6[4,3]<-as.numeric(as.character(summary(anova.change.deceleration.6)[[2]][[1]][[1]])) anova.table.deceleration.6[4,4]<-round(as.numeric(as.character(summary(anova.change.deceleration.6)[[2]][[1]][[3]])), digits=5) anova.table.deceleration.6[5,1]<-"Error" anova.table.deceleration.6[5,2]<-round(as.numeric(as.character(summary(anova.change.deceleration.6)[[3]][[1]][[2]])), digits=5) anova.table.deceleration.6[5,3]<-as.numeric(as.character(summary(anova.change.deceleration.6)[[3]][[1]][[1]])) anova.table.deceleration.6[5,4]<-round(as.numeric(as.character(summary(anova.change.deceleration.6)[[3]][[1]][[3]])), digits=5) anova.table.deceleration.6[6,1]<-"Total" anova.table.deceleration.6[3,5]<-round(as.numeric(as.character(summary(anova.change.deceleration.6)[[1]][[1]][[3]]))/as.numeric(as.character(summary(anova.change.deceleration.6)[[3]][[1]][[3]])), digits=5) anova.table.deceleration.6[3,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.deceleration.6[3,5])), as.numeric(as.character(summary(anova.change.deceleration.6)[[1]][[1]][[1]])), as.numeric(as.character(summary(anova.change.deceleration.6)[[3]][[1]][[1]])) ))), digits=4) anova.table.deceleration.6[3,7]<-round(F.p.deceleration.6, digits=5) anova.table.deceleration.6[4,5]<-round(as.numeric(as.character(summary(anova.change.deceleration.6)[[2]][[1]][[3]]))/as.numeric(as.character(summary(anova.change.deceleration.6)[[3]][[1]][[3]])), digits=5) anova.table.deceleration.6[4,6]<-round(as.numeric(as.character(1-pf(as.numeric(as.character(anova.table.deceleration.6[4,5])), as.numeric(as.character(summary(anova.change.deceleration.6)[[2]][[1]][[1]])), as.numeric(as.character(summary(anova.change.deceleration.6)[[3]][[1]][[1]])) ))), digits=4) if(as.numeric(as.character(anova.table.deceleration.6[3,6]))<0.0001) { anova.table.deceleration.6[3,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.deceleration.6[4,6]))<0.0001) { anova.table.deceleration.6[4,6]<-"< 0.0001"} if(as.numeric(as.character(anova.table.deceleration.6[3,7]))<0.0001) { anova.table.deceleration.6[3,7]<-"< 0.0001"} anova.table.deceleration.6[6,2]<-as.numeric(as.character(anova.table.deceleration.6[3,2]))+as.numeric(as.character(anova.table.deceleration.6[4,2]))+as.numeric(as.character(anova.table.deceleration.6[5,2])) anova.table.deceleration.6[6,3]<-as.numeric(as.character(anova.table.deceleration.6[3,3]))+as.numeric(as.character(anova.table.deceleration.6[4,3]))+as.numeric(as.character(anova.table.deceleration.6[5,3])) anova.table.deceleration.6[6,4]<-round(as.numeric(as.character(anova.table.deceleration.6[6,2]))/as.numeric(as.character(anova.table.deceleration.6[6,3])), digits=5) Table.9<-rbind(anova.table.deceleration.9, anova.table.deceleration.6) Table.9<-as.data.frame(Table.9) # Tables were renumbered on submission, as many were moved to the Supplement. What used to be Table 9 is now Table S7. write.csv(file="Table.S7.csv", Table.9) ####### # Now we want to look at the in-text calculations of correlation between fitness gain # over an interval, and initial fitness at the start of that interval. initial.fitness.1<-Table.1[1:9, 3] initial.fitness.1<-as.numeric(as.character(initial.fitness.1)) change.1<-Table.2[1:9, 3] change.1<-as.numeric(as.character(change.1)) cor(initial.fitness.1, change.1) # Now we exclude Ara+1: initial.fitness.2<-initial.fitness.1[2:9] change.2<-change.1[2:9] cor(initial.fitness.2, change.2) cor.test(initial.fitness.2, change.2, alternative=c("two.sided")) ##################### # Finally, we want to create a graphical version of Table 1 Figure.1.Gens<-c(40000, 50000, 60000) Figure.1.A.Plus.1<-c(as.numeric(as.character(Table.1[1, 3])), as.numeric(as.character(Table.1[1, 6])), as.numeric(as.character(Table.1[1, 9]))) Figure.1.A.Plus.2<-c(as.numeric(as.character(Table.1[2, 3])), as.numeric(as.character(Table.1[2, 6])), as.numeric(as.character(Table.1[2, 9]))) Figure.1.A.Plus.3<-c(as.numeric(as.character(Table.1[3, 3])), as.numeric(as.character(Table.1[3, 6])), as.numeric(as.character(Table.1[3, 9]))) Figure.1.A.Plus.4<-c(as.numeric(as.character(Table.1[4, 3])), as.numeric(as.character(Table.1[4, 6])), as.numeric(as.character(Table.1[4, 9]))) Figure.1.A.Plus.5<-c(as.numeric(as.character(Table.1[5, 3])), as.numeric(as.character(Table.1[5, 6])), as.numeric(as.character(Table.1[5, 9]))) Figure.1.A.Minus.1<-c(as.numeric(as.character(Table.1[6, 3])), as.numeric(as.character(Table.1[6, 6])), as.numeric(as.character(Table.1[6, 9]))) Figure.1.A.Minus.4<-c(as.numeric(as.character(Table.1[7, 3])), as.numeric(as.character(Table.1[7, 6])), as.numeric(as.character(Table.1[7, 9]))) Figure.1.A.Minus.5<-c(as.numeric(as.character(Table.1[8, 3])), as.numeric(as.character(Table.1[8, 6])), as.numeric(as.character(Table.1[8, 9]))) Figure.1.A.Minus.6<-c(as.numeric(as.character(Table.1[9, 3])), as.numeric(as.character(Table.1[9, 6])), as.numeric(as.character(Table.1[9, 9]))) Figure.1.Means<-c(as.numeric(as.character(Table.1[10, 3])), as.numeric(as.character(Table.1[10, 6])), as.numeric(as.character(Table.1[10,9]))) Figure.1.A.Plus.1<-as.numeric(as.character(Figure.1.A.Plus.1)) Figure.1.A.Plus.2<-as.numeric(as.character(Figure.1.A.Plus.2)) Figure.1.A.Plus.3<-as.numeric(as.character(Figure.1.A.Plus.3)) Figure.1.A.Plus.4<-as.numeric(as.character(Figure.1.A.Plus.4)) Figure.1.A.Plus.5<-as.numeric(as.character(Figure.1.A.Plus.5)) Figure.1.A.Minus.1<-as.numeric(as.character(Figure.1.A.Minus.1)) Figure.1.A.Minus.4<-as.numeric(as.character(Figure.1.A.Minus.4)) Figure.1.A.Minus.5<-as.numeric(as.character(Figure.1.A.Minus.5)) Figure.1.A.Minus.6<-as.numeric(as.character(Figure.1.A.Minus.6)) Figure.1.Means<-as.numeric(as.character(Figure.1.Means)) graphics.off() pdf(file="Figure.1.pdf", width=6.9, height=6.9) par(mfrow=c(1,1), oma=c(3,3, 1.5, 1.5), bty="n", xpd=NA, mar=c(3, 3, 1, 1)) plot(x=Figure.1.Gens, y=Figure.1.Means, xlim=c(40000, 60000), ylim=c(0.90, 1.15), axes=FALSE, xlab="", ylab="", col="white") lines(x=Figure.1.Gens, y=Figure.1.A.Plus.1, lwd=2, col="darkred") lines(x=Figure.1.Gens, y=Figure.1.A.Plus.2, lwd=2, col="darkorange") lines(x=Figure.1.Gens, y=Figure.1.A.Plus.3, lwd=2, col="mediumblue") lines(x=Figure.1.Gens, y=Figure.1.A.Plus.4, lwd=2, col="forestgreen") lines(x=Figure.1.Gens, y=Figure.1.A.Plus.5, lwd=2, col="navyblue") lines(x=Figure.1.Gens, y=Figure.1.A.Minus.1, lwd=2, col="purple4") lines(x=Figure.1.Gens, y=Figure.1.A.Minus.4, lwd=2, col="darkgoldenrod") lines(x=Figure.1.Gens, y=Figure.1.A.Minus.5, lwd=2, col="chocolate4") lines(x=Figure.1.Gens, y=Figure.1.A.Minus.6, lwd=2, col="violetred4") lines(x=Figure.1.Gens, y=Figure.1.Means, lwd=5, col="black", lty=5) #points(x=Figure.1.Gens, y=Figure.1.Means, pch=16, col="black") axis(2, at=c(0.90, 1.15), labels=c("", ""), lwd.ticks=0, pos=0) axis(1, at=c(40000, 50000, 60000), labels=c("40,000", "50,000", "60,000"), pos=0.90) axis(2, at=c(0.90, 0.95, 1.00, 1.05, 1.10, 1.15), labels=c("0.90", "0.95", "1.00", "1.05", "1.10", "1.15"), las=2, pos=40000) mtext("Time (generations) ", side=1, outer=TRUE, cex=1.2) mtext("Relative fitness ", side=2, outer=TRUE, cex=1.2) dev.off() ######## # Now we have been asked to do an additional analysis. # For each of the populations in this data: # 1) Use Wiser et al data for that lineage to fit power law. # 2) Find predicted value of each curve at 60k generations. # 3) Divide all values in 2 by predicted value at 40k for Ara-5. # 4) Find correlation between 60k in new data set and values in 3. # Step 1: Get all predicted curves Published.data <- read.table("Concatenated.LTEE.data.all.csv", sep=",", header=T) power.extractor <- function(Pop) { just.pop<-subset(Published.data, Population==Pop) just.pop<-na.omit(just.pop) Gen.list<-unique(just.pop$Generation) x<-seq(from=0, to=max(Gen.list), by=50) pop.5<-subset(just.pop, Generation<5001) pop.10<-subset(just.pop, Generation<10001) pop.20<-subset(just.pop, Generation<20001) pop.30<-subset(just.pop, Generation<30001) pop.40<-subset(just.pop, Generation<40001) pop.50<-subset(just.pop, Generation<50001) model.2.50<-nls(Fitness ~ (b*Generation + 1)^a, data=pop.50, start=list(a=0.1, b=0.01), algorithm="port", trace=F, lower=c(.01, .0005), upper=c(.4, 0.1), na.action=na.omit, model=T, control=nls.control(maxiter=1000, warnOnly=F)) model.2.50.a<-coef(model.2.50)[1] model.2.50.b<-coef(model.2.50)[2] fit.value.40 <- (model.2.50.b*40000 + 1)^model.2.50.a fit.value.50 <- (model.2.50.b*50000 + 1)^model.2.50.a fit.value.60 <- (model.2.50.b*60000 + 1)^model.2.50.a output.matrix <- matrix(nrow=1, ncol=6) output.matrix[1, 1] <- Pop output.matrix[1, 2] <- model.2.50.a output.matrix[1, 3] <- model.2.50.b output.matrix[1, 4] <- fit.value.40 output.matrix[1, 5] <- fit.value.50 output.matrix[1, 6] <- fit.value.60 return(output.matrix) } values.power.Plus.1 <- power.extractor("Ara + 1") values.power.Plus.2 <- power.extractor("Ara + 2") values.power.Plus.3 <- power.extractor("Ara + 3") values.power.Plus.4 <- power.extractor("Ara + 4") values.power.Plus.5 <- power.extractor("Ara + 5") values.power.Minus.1 <- power.extractor("Ara - 1") values.power.Minus.4 <- power.extractor("Ara - 4") values.power.Minus.5 <- power.extractor("Ara - 5") values.power.Minus.6 <- power.extractor("Ara - 6") # Step 2: Find fitted values: fitted.power <- rbind(values.power.Plus.1, values.power.Plus.2, values.power.Plus.3, values.power.Plus.4, values.power.Plus.5, values.power.Minus.1, values.power.Minus.4, values.power.Minus.5, values.power.Minus.6) colnames(fitted.power) <- c("Population", "Parameter.a", "Parameter.b", "Fitted.value.40", "Fitted.value.50", "Fitted.value.60") fitted.frame <- as.data.frame(fitted.power) fitted.frame$Fitted.value.40<-as.numeric(as.character(fitted.frame$Fitted.value.40)) fitted.frame$Fitted.value.50<-as.numeric(as.character(fitted.frame$Fitted.value.50)) fitted.frame$Fitted.value.60<-as.numeric(as.character(fitted.frame$Fitted.value.60)) # Step 3: Divide by value in population Ara-5 at 40k. divisor <- as.numeric(as.character(values.power.Minus.5[1, 4])) scaled.power <- cbind(fitted.power[,1], fitted.power[,2], fitted.power[,3], fitted.frame$Fitted.value.40/divisor, fitted.frame$Fitted.value.50/divisor, fitted.frame$Fitted.value.60/divisor) colnames(scaled.power) <- c("Population", "Parameter.a", "Parameter.b", "Scaled.40", "Scaled.50", "Scaled.60") scaled.power.frame <- as.data.frame(scaled.power) scaled.power.frame$Parameter.a <- as.numeric(as.character(scaled.power.frame$Parameter.a)) scaled.power.frame$Parameter.b <- as.numeric(as.character(scaled.power.frame$Parameter.b)) scaled.power.frame$Scaled.40 <- as.numeric(as.character(scaled.power.frame$Scaled.40)) scaled.power.frame$Scaled.50 <- as.numeric(as.character(scaled.power.frame$Scaled.50)) scaled.power.frame$Scaled.60 <- as.numeric(as.character(scaled.power.frame$Scaled.60)) comparison.table <- cbind(as.character(Table.1[1:9,1]), as.numeric(as.character(Table.1[1:9,9])), as.character(scaled.power.frame$Population), as.numeric(as.character(scaled.power.frame$Scaled.60))) colnames(comparison.table) <- c("Population.Measured", "Measured.Fitness", "Population.Predicted", "Predicted.Fitness") comparison.frame <- as.data.frame(comparison.table) comparison.frame$Measured.Fitness <- as.numeric(as.character(comparison.frame$Measured.Fitness)) comparison.frame$Predicted.Fitness <- as.numeric(as.character(comparison.frame$Predicted.Fitness)) # Now we find the correlation between the measured and predicted values: cor.test(comparison.frame$Measured.Fitness, comparison.frame$Predicted Fitness, alternative=c("two.sided")