#this program builds two datasets - the spring dataset and muscle dataset #then it combines the two into one dataset #then it looks at scaling of both parameters within species #then it moves to phylogeny-based analyses looking at cross-species scaling and correlates of appendage type #open libraries relevant to analyses library(ape) library(geiger) library(picante) library(caper) library(apTreeshape) library(ggplot2) library(splines) #clear all variables from workspace rm(list=ls()) #set the directory and open the files setwd("C:/Users/snp2/Documents/stomatopods/muscles/2013 muscles/FINAL data to datadryad") ##springdataset alldatatoparse=read.table(file="springdataset.txt",header=T,sep="\t") ##muscledataset allmuscle=read.table(file="muscledataset.txt",header=T,sep="\t") ################################################################################ #### Making SPRING dataset alldatatoparse=read.table(file="compdataall112310v1.txt",header=T,sep="\t") #only look at 100% pulls (no partial ones) alldata = subset(alldatatoparse, alldatatoparse$pull==0) #all following have 100% pulls (partials have been removed) #make proport max slope into percent alldata[,15]=alldata$proportmaxslope*100 #take a look at the data names(alldata) #[1] "Family" "Genus" "Species" "ID" "Appendage" "Trial" # [7] "Saddle_Status" "Date.of.measurement" "Sex" "database.filename" "Remarks.re.photo.msmts" "pull" #[13] "raw.filename" "imaxslope" "proportmaxslope" "maxslope" "maxslopeintercept" "maxsloper2d" #[19] "slope60_90" "intercept60_90" "r.2.adj.loading.60.95." "slope5_40" "intercept5_40" "r.2.adj.loading.5.40." #[25] "slope_all" "intercept_all" "r.2.adj.loading.all.curve" "slopeunload5_15" "interceptunload5_15" "r.2.adj.unloading.5.15." #[31] "slopeunloadall" "interceptunloadall" "r.2.adj.unloading.all.curve" "maxexten" "maxforce" "workload" #[37] "workunload" "resilience" "logmaxexten" "logmaxforce" "logworkload" "logworkunload" #[43] "bodymass" "cpl" "bodylen" "Lmeruswidth" "Ldacpropwidth" "Rmeruswidth" #[49] "Rdactpropwidth" "merus.length" "V.length" "logmeruslen" "logvlen" #to make life easier, just include columns with data we're currently interested in allcolumns = c(1:7, 9, 14:17, 34:53) reldata = alldata[,allcolumns] # names(reldata) # [1] "Family" "Genus" "Species" "ID" "Appendage" "Trial" "Saddle_Status" "Sex" # [9] "imaxslope" "proportmaxslope" "maxslope" "maxslopeintercept" "maxexten" "maxforce" "workload" "workunload" #[17] "resilience" "logmaxexten" "logmaxforce" "logworkload" "logworkunload" "bodymass" "cpl" "bodylen" #[25] "Lmeruswidth" "Ldacpropwidth" "Rmeruswidth" "Rdactpropwidth" "merus.length" "V.length" "logmeruslen" "logvlen" #separate out intact saddle data, and males and females subsets (not using severed saddles in this study) intact = subset(reldata, reldata$Saddle_Status=="intact") #intact males intactmales = subset(intact, intact$Sex=="m") #intact females intactfemales = subset(intact, intact$Sex=="f") ##SADDLES INTACT ALL (males + females_ #start calculating descriptive statistics for each species with intact saddles, lumping trials within individual Species <- factor(intact$Species) ID <- factor(intact$ID) Trial <- factor(intact$Trial) Appendage <- factor(intact$Appendage) #gives means for all individuals sz = ncol(intact) intactindivid = aggregate(intact[,9:sz],list(Species, ID), mean) #need to get rid of chiragra 1, because it is missing most body size data intactindivid = intactindivid[-1,] #now get means for species with intact saddles Species2 <- factor(intactindivid$Group.1) sz = ncol(intactindivid) intspeciesmeans = aggregate(intactindivid[,3:sz],list(Species2), mean) intspeciesds = aggregate(intactindivid[,3:sz],list(Species2), FUN = "sd") intspeciesmin = aggregate(intactindivid[,3:sz],list(Species2), FUN = "min") intspeciesmax = aggregate(intactindivid[,3:sz],list(Species2), FUN = "max") #here are the individual counts for intact males and females -- I checked an these all add up correctly with the males and females separated out below # bredini : 1 # californiensis : 3 # chiragra : 9 # empusa: 2 # falcatus : 18 # festae : 1 # maculata : 5 # oerstedii : 6 # smithii : 11 # sulcata : 3 # wennerae : 2 #remove festae because it is not in the phylogeny spdataMF = intspeciesmeans[-6,] spdataMF <- data.frame(spdataMF) ##SADDLES INTACT FEMALES #start calculating descriptive statistics for each species with intact saddles, lumping trials within individual SpeciesIF <- factor(intactfemales$Species) IDIF <- factor(intactfemales$ID) TrialIF <- factor(intactfemales$Trial) AppendageIF <- factor(intactfemales$Appendage) #gives individual means sz = ncol(intactfemales) intactindividIF = aggregate(intactfemales[,9:sz],list(SpeciesIF, IDIF), mean) #need to get rid of chiragra 1, because it is missing most body size data intactindividIF = intactindividIF[-1,] #now get means for species with intact saddles Species2IF <- factor(intactindividIF$Group.1) sz = ncol(intactindividIF) intspeciesmeansIF = aggregate(intactindividIF[,3:sz],list(Species2IF), mean) intspeciesdsIF = aggregate(intactindividIF[,3:sz],list(Species2IF), FUN = "sd") intspeciesminIF = aggregate(intactindividIF[,3:sz],list(Species2IF), FUN = "min") intspeciesmaxIF = aggregate(intactindividIF[,3:sz],list(Species2IF), FUN = "max") #here are the individual counts for intact females # chiragra : 4 # falcatus : 12 # maculata : 1 # oerstedii : 5 # smithii : 6 # sulcata : 1 # wennerae : 1 spdataIF = intspeciesmeansIF spdataIF <- data.frame(spdataIF) #SADDLES INTACT MALES #start calculating descriptive statistics for each species with intact saddles, lumping trials within individual SpeciesIM <- factor(intactmales$Species) IDIM <- factor(intactmales$ID) TrialIM <- factor(intactmales$Trial) AppendageIM <- factor(intactmales$Appendage) #gives individual means sz = ncol(intactmales) intactindividIM = aggregate(intactmales[,9:sz],list(SpeciesIM, IDIM), mean) #now get means for species with intact saddles Species2IM <- factor(intactindividIM$Group.1) sz = ncol(intactindividIM) intspeciesmeansIM = aggregate(intactindividIM[,3:sz],list(Species2IM), mean) intspeciesdsIM = aggregate(intactindividIM[,3:sz],list(Species2IM), FUN = "sd") intspeciesminIM = aggregate(intactindividIM[,3:sz],list(Species2IM), FUN = "min") intspeciesmaxIM = aggregate(intactindividIM[,3:sz],list(Species2IM), FUN = "max") #here are the individual counts for intact males: # bredini:1 # californiensis:3 # chiragra:5 # empusa:2 # falcatus:6 # festae: 1 # maculata:4 # oerstedii:1 # smithii:5 # sulcata:2 # wennerae:1 #remove festae because it is not in the phylogeny spdataIM = intspeciesmeansIM[-6,] spdataIM <- data.frame(spdataIM) ############################################################################### ##making MUSCLE datasets malemuscles = subset(allmuscle, allmuscle$Sex=="m") malemusclesDF = data.frame(malemuscles) femalemuscles = subset(allmuscle, allmuscle$Sex=="f") femalemusclesDF = data.frame(femalemuscles) #sort matrix by species name allmuscleDF = data.frame(allmuscle) speciesort = allmuscleDF[order(allmuscleDF$Species), ] #these are just at species level sz = ncol(allmuscle) musclespeciesALL = aggregate(allmuscleDF[,8:sz], list(allmuscleDF$Species),mean, na.rm = TRUE) musclespeciesMALES = aggregate(malemusclesDF[,8:sz], list(malemusclesDF$Species),mean, na.rm = TRUE) musclespeciesFEMALES = aggregate(femalemusclesDF[,8:sz], list(femalemusclesDF$Species),mean, na.rm = TRUE) #parse out individual species data # # austrinus = 4, 1 m, 3 f # bahiahondensis = 1 f # bredini = 13, 12f, 1m # californiensis = 13, 12m, 1f # chiragra = 7, 4m 3f #iffy female, male comparison # ciliata =18, 9m, 7f #definite female male comparison # empusa = 2m # falcatus =14, 5m, 9f # good f m comparison # # [9] festae = 12, 1f, 11 m # hieroglyphica = 1m # latirostris 4, 3f, 1m # maculata=5, 4m ,1f # oerstedii =7 2f 5m # ornata = 2, 1f, 1m # oxyrhyncha =2. 1f, 1m # pygmaea = 3f # # [17] smithii = 7, 3m 4f #iffy female male comparison # spinocarinatus = 1m # sulcata = 1m # tweediei = 5, 1m , 4f # vicina =3, 2f, 1m # wennerae =5, 1m, 4f #combo males + females #commented out taxa with too few individs for scaling austrinusALL = subset(allmuscle, allmuscle$Species=="austrinus") #bahiahondensisALL = subset(allmuscle, allmuscle$Species=="bahiahondensis") brediniALL = subset(allmuscle, allmuscle$Species=="bredini") californiensisALL = subset(allmuscle, allmuscle$Species=="californiensis") chiragraALL = subset(allmuscle, allmuscle$Species=="chiragra") ciliataALL = subset(allmuscle, allmuscle$Species=="ciliata") #empusaALL = subset(allmuscle, allmuscle$Species=="empusa") falcatusALL = subset(allmuscle, allmuscle$Species=="falcatus") festaeALL = subset(allmuscle, allmuscle$Species=="festae") #hieroglyphicaALL = subset(allmuscle, allmuscle$Species=="hieroglyphica") latirostrisALL = subset(allmuscle, allmuscle$Species=="latirostris") maculataALL = subset(allmuscle, allmuscle$Species=="maculata") oerstediiALL = subset(allmuscle, allmuscle$Species=="oerstedii") #ornataALL = subset(allmuscle, allmuscle$Species=="ornata") #oxyrhynchaALL = subset(allmuscle, allmuscle$Species=="oxyrhyncha") #pygmaeaALL = subset(allmuscle, allmuscle$Species=="pygmaea") smithiiALL = subset(allmuscle, allmuscle$Species=="smithii") #spinocarinatusALL = subset(allmuscle, allmuscle$Species=="spinocarinatus") #sulcataALL = subset(allmuscle, allmuscle$Species=="sulcata") tweedieiALL = subset(allmuscle, allmuscle$Species=="tweediei") #vicinaALL = subset(allmuscle, allmuscle$Species=="vicina") wenneraeALL = subset(allmuscle, allmuscle$Species=="wennerae") ############################################################################## #make a table of body sizes for paper #now get means for species with intact saddles sz = ncol(allmuscle) musclespeciessd = aggregate(allmuscleDF[,8:sz], list(allmuscleDF$Species),sd, na.rm = TRUE) write.table(musclespeciessd, "clipboard", sep="\t") musclespeciesmean = aggregate(allmuscleDF[,8:sz], list(allmuscleDF$Species),mean, na.rm = TRUE) write.table(musclespeciesmean, "clipboard", sep="\t") musclespeciesmin = aggregate(allmuscleDF[,8:sz], list(allmuscleDF$Species),min, na.rm = TRUE) write.table(musclespeciesmin, "clipboard", sep="\t") musclespeciesmax = aggregate(allmuscleDF[,8:sz], list(allmuscleDF$Species),max, na.rm = TRUE) write.table(musclespeciesmax, "clipboard", sep="\t") #set decimal places to reflect accuracy SLmean = format(round(musclespeciesmean$MeanmeanSL,2)) SLDF = format(round(musclespeciessd$MeanmeanSL,2)) SLmin = format(round(musclespeciesmin$MeanmeanSL,2)) SLmax = format(round(musclespeciesmax$MeanmeanSL,2)) tabledata = paste(SLmean, "(", SLmin, "-", SLmax,")", "+", SLDF) write.table(tabledata, "clipboard", sep="\t") damean = format(round(musclespeciesmean$MeanPCSAisAsin2amm2,2)) daDF = format(round(musclespeciessd$MeanPCSAisAsin2amm2,2)) damin = format(round(musclespeciesmin$MeanPCSAisAsin2amm2,2)) damax = format(round(musclespeciesmax$MeanPCSAisAsin2amm2,2)) tabledata = paste(damean, "(", damin, "-", damax,")", "+", daDF) write.table(tabledata, "clipboard", sep="\t") ############################################################################### #MERGE muscle and spring datasets #first merge species without regard for ID colnames(musclespeciesALL)[1] <- "Species" colnames(spdataMF)[1] <- "Species" mergeSPECIES = merge(musclespeciesALL, spdataMF, by=c("Species")) #this variable has all of the overlapping species data for muscles and springs #gives 10 species #now try to link individual data together #goal is to use muscle dataset as limiting factor and to include only spring data that match individuals from muscle data #this dataset would be used for looking at scaling between muscle and spring parameters colnames(intactindivid)[1] <- "Species" colnames(intactindivid)[2] <- "ID" intactindividDF = data.frame(intactindivid) intactindividDF = intactindividDF[order(intactindivid$Species), ] allmuscleDF = data.frame(allmuscleDF) allmuscleDF = allmuscleDF[order(allmuscleDF$Family, allmuscleDF$Genus, allmuscleDF$Species), ] mergeINDIVID = merge(allmuscleDF, intactindividDF,by=c("Species", "ID")) #this variable gives overlapping individ data - results in 5 species, 2-6 individ/species ############################################################################### ## Look at the data #let's do this in an efficient way using ggplot LOVE IT! #add a column to make continuous appendage type into categorical #or can just say factor(variable) and will make it categorical musclespeciesALL["AppenType"] <- NA musclespeciesALL$AppenType[musclespeciesALL$undif0spearer1smash2==0]<-"U" musclespeciesALL$AppenType[musclespeciesALL$undif0spearer1smash2==1]<-"P" musclespeciesALL$AppenType[musclespeciesALL$undif0spearer1smash2==2]<-"M" #this is unrefined in qplot #pdf(file = "unloggedspeciesgraphs.pdf", width = 6, height = 6) # # If inside a script, you will need to explicitly print() plots # qplot(mpg, wt, data = mtcars) # qplot(wt, mpg, data = mtcars) j = 10#y variable of interest for(i in 3:6){ xname = colnames(musclespeciesALL)[i] yname =colnames(musclespeciesALL)[j] xdata = musclespeciesALL[[i]] ydata = musclespeciesALL[[j]] c <-qplot(xdata, ydata, data=musclespeciesALL, colour= AppenType, shape= AppenType, geom = c("point", "smooth"), method = "lm", ) #put the data on the graph muscleplots <- c + labs(x=xname, y=yname) +theme(axis.title.x=element_text(size=20, face = "bold",family = "sans"), axis.title.y=element_text(size=20, face = "bold",family = "sans"), axis.text.x = element_text(angle=90, vjust=0.5, size=16), axis.text.y = element_text(angle=90, vjust=0.5, size=16)) #adjust the appearance print (muscleplots) #make it show up on the screen #plot(musclespeciesALL[[i]], musclespeciesALL[[j]],xlab = xname, ylab=yname ) #Sys.sleep(2) #pauses for 1 second } #dev.off() #this is using ggplot #pdf(file = "loggedspeciesgraphsv4.pdf", width = 6, height = 6) #now do the same thing logged j = 18 #y variable of interest (7:10 unlogged variables, 15:18 logged) for(i in 16:18){ xname = colnames(musclespeciesALL)[i] yname =colnames(musclespeciesALL)[j] xdata = musclespeciesALL[[i]] ydata = musclespeciesALL[[j]] p <- ggplot(musclespeciesALL, aes(xdata, ydata,colour= AppenType, shape= AppenType))+geom_point(aes(size=2))+ geom_smooth(method="lm") muscleplots <- p + labs(x=xname, y=yname) +theme(axis.title.x=element_text(size=20, face = "bold",family = "sans"), axis.title.y=element_text(size=20, face = "bold",family = "sans"), axis.text.x = element_text(angle=90, vjust=0.5, size=16), axis.text.y = element_text(angle=90, vjust=0.5, size=16) ) #adjust the appearance print(muscleplots) } #dev.off() #merus scaling #pdf(file = "loggedmerusscaling.pdf", width = 6, height = 6) j = 14 #y variable of interest (7:10 unlogged variables, 15:18 logged) for(i in 11:13){ xname = colnames(musclespeciesALL)[i] yname =colnames(musclespeciesALL)[j] xdata = musclespeciesALL[[i]] ydata = musclespeciesALL[[j]] p <- ggplot(musclespeciesALL, aes(xdata, ydata,colour= AppenType, shape= AppenType))+geom_point(aes(size=2))+ geom_smooth(method="lm") muscleplots <- p + labs(x=xname, y=yname) +theme(axis.title.x=element_text(size=20, face = "bold",family = "sans"), axis.title.y=element_text(size=20, face = "bold",family = "sans"), axis.text.x = element_text(angle=90, vjust=0.5, size=16), axis.text.y = element_text(angle=90, vjust=0.5, size=16) ) print(muscleplots) } #dev.off() #plot muscle and springdata #mergeSPECIES provides merged averaged species data #mergeINDIVID provides data for which individuals are matched #pdf(file = "musclespringgraphs2nd.pdf", width = 6, height = 6) AppenType = factor(mergeSPECIES$undif0spearer1smash2) j = 18 #y variable of interest (17:18 logged PCSA and SL) for(i in 27:31){ xname = colnames(mergeSPECIES)[i] yname =colnames(mergeSPECIES)[j] xdata = mergeSPECIES[[i]] ydata = mergeSPECIES[[j]] p <- ggplot(mergeSPECIES, aes(xdata, ydata,colour= AppenType, shape= AppenType))+geom_point(aes(size=2))+ geom_smooth(method="lm") muscleplots <- p + labs(x=xname, y=yname) +theme(axis.title.x=element_text(size=20, face = "bold",family = "sans"), axis.title.y=element_text(size=20, face = "bold",family = "sans"), axis.text.x = element_text(angle=90, vjust=0.5, size=16), axis.text.y = element_text(angle=90, vjust=0.5, size=16) ) print(muscleplots) } #dev.off() ############################################################################### #time to do some stats #look at distributions for phylogeny testing - test for normality, etc. #make a loop to get all the p-values for a Shapiro Wilk test for each species specname = mergeSPECIES #did this for musclespeciesFEMALES and musclespeciesALL #find the size of the dataset to be scanned i=3 sizem = ncol(specname) shap=vector(length=sizem) #loop through size of dataset and enter result into the appropriate column for(i in 3:sizem) { res = shapiro.test(specname[[i]]) shap[i] = res$p.value } shap #make it easier to copy and paste in Excel write.table(shap, "clipboard", sep="\t") specname = smithiiALL #cycled through the four species with enough data #find the size of the dataset to be scanned i=9 sizem = ncol(specname) shap=vector(length=sizem) #loop through size of dataset and enter result into the appropriate column for(i in 9:sizem) { res = shapiro.test(specname[,i]) shap[i] = res$p.value } shap #do m-f t-test on chiragra, ciliata, falcatus, smithii #do t-test comparing males and females within species males = subset(smithiiALL,smithiiALL$Sex == 'm') females = subset(smithiiALL, smithiiALL$Sex == 'f') #tresult$parameter = the df of test tresult = t.test(males$MeanAvgmeruslengthmm, females$MeanAvgmeruslengthmm) mlresult = cbind(tresult$data.name, tresult$statistic, tresult$parameter, tresult$p.value) tresult = t.test(males$Meanavgcarapacelengthmm,females$Meanavgcarapacelengthmm) cplresult = cbind(tresult$data.name, tresult$statistic, tresult$parameter, tresult$p.value) tresult= t.test(males$Meanavgpinnationangledeg,females$Meanavgpinnationangledeg) pinangleresult = cbind(tresult$data.name, tresult$statistic, tresult$parameter, tresult$p.value) tresult = t.test(males$Meanavgapodemeareamm2,females$Meanavgapodemeareamm2) apoarearesult = cbind(tresult$data.name, tresult$statistic, tresult$parameter, tresult$p.value) tresult = t.test(males$MeanPCSAisAsin2amm2,females$MeanPCSAisAsin2amm2) pcsaresult = cbind(tresult$data.name, tresult$statistic, tresult$parameter, tresult$p.value) tresult = t.test(males$MeanmeanSL,females$MeanmeanSL) SLresult = cbind(tresult$data.name, tresult$statistic, tresult$parameter, tresult$p.value) tresult = rbind(mlresult, cplresult, pinangleresult, apoarearesult, pcsaresult, SLresult) tresult #add a couple more analyses males = subset(falcatusALL,falcatusALL$Sex == 'm') females = subset(falcatusALL,falcatusALL$Sex == 'f') tresult = t.test(males$Meanavgwetweightg, females$Meanavgwetweightg) massresult = cbind(tresult$data.name, tresult$statistic, tresult$parameter, tresult$p.value) tresult = t.test(males$Meanavgbodylengthmm, females$Meanavgbodylengthmm) blresult = cbind(tresult$data.name, tresult$statistic, tresult$parameter, tresult$p.value) tresult = rbind(massresult, blresult) tresult #make it easier to copy and paste in Excel write.table(tresult, "clipboard", sep="\t") #linear regressions with species data #each data point is a species dataset = musclespeciesALL overallresults=matrix(nrow=4, ncol=8) j = 14 #y variable of interest (7:10 unlogged variables, 15:18 logged) beg = 11 end= 13 for(i in beg:end){ xname = colnames(dataset)[i] yname =colnames(dataset)[j] xdata = dataset[[i]] ydata = dataset[[j]] categ = factor(dataset[[2]]) model<-lm(ydata~xdata) results <-summary(model) intercept = summary(model)$coefficients[1] slope = summary(model)$coefficients[2] fstat = summary(model)$fstatistic[1] df1 = summary(model)$fstatistic[2] df2 = summary(model)$fstatistic[3] pval = summary(model)$coefficients[8] rsqd = summary(model)$r.squared rsqdadj = summary(model)$adj.r.squared lmresults = cbind(intercept, slope, fstat, df1, df2,pval, rsqd, rsqdadj) overallresults [i-beg+1,] = lmresults plot(ydata~xdata, sub = pval) abline(model) text(xname, yname, categ, cex=0.5, pos=3, col="blue") } overallresults #make it easier to copy and paste in Excel write.table(overallresults, "clipboard", sep="\t") ######### # #linear regressions within species: # falcatusALL # chiragraALL # ciliataALL # smithiiALL #use unlogged data here, because the unlogged data are normally distributed, whereas logged are not dataset = ciliataALL overallresults=matrix(nrow=4, ncol=10) j = 16 #y variable of interest (15 is PCSA, 16 is SL) beg = 10 end= 12 for(i in beg:end){ xname = colnames(dataset)[i] yname =colnames(dataset)[j] xdata = dataset[[i]] ydata = dataset[[j]] categ = factor(dataset[[2]]) model<-lm(ydata~xdata) results <-summary(model) intercept = summary(model)$coefficients[1] slope = summary(model)$coefficients[2] fstat = summary(model)$fstatistic[1] df1 = summary(model)$fstatistic[2] df2 = summary(model)$fstatistic[3] pval = summary(model)$coefficients[8] rsqd = summary(model)$r.squared rsqdadj = summary(model)$adj.r.squared lmresults = cbind(yname, xname, intercept, slope, fstat, df1, df2,pval, rsqd, rsqdadj) overallresults [i-beg+1,] = lmresults plot(ydata~xdata, sub = pval) abline(model) text(xname, yname, categ, cex=0.5, pos=3, col="blue") } overallresults #make it easier to copy and paste in Excel write.table(overallresults, "clipboard", sep="\t") #compare regression slopes and distribution for ciliata and falcatus to see if scaling differs #http://r-eco-evo.blogspot.com/2011/08/comparing-two-regression-slopes-by.html #combine the two datasets cilfalcompare = rbind(ciliataALL, falcatusALL) cilfalcompare = data.frame(cilfalcompare) mod1 <- aov(MeanPCSAisAsin2amm2~MeanAvgmeruslengthmm*Species, data=cilfalcompare) summary(mod1) mod2 <- aov(MeanPCSAisAsin2amm2~MeanAvgmeruslengthmm+Species, data=cilfalcompare) summary(mod2) anova(mod1,mod2) #look at spring scaling and muscle scaling within and across species smithiimerge = subset(mergeINDIVID, mergeINDIVID$Species=="smithii") chiragramerge = subset(mergeINDIVID, mergeINDIVID$Species=="chiragra") maculatamerge = subset(mergeINDIVID, mergeINDIVID$Species=="maculata") #want to look at workload (32),workunload (33), maxslope (28), maxforce (31), resilience (34), maxexten(30) #in context of PCSA (15) and SL (17); how to incorporate body size? what is x v. y variable? dataset = chiragramerge overallresults=matrix(nrow=7, ncol=10) j = 15 #y variable of interest (13:16) beg = 27 end= 33 for(i in beg:end){ xname = colnames(dataset)[i] yname =colnames(dataset)[j] xdata = dataset[[i]] ydata = dataset[[j]] categ = factor(dataset[[2]]) model<-lm(ydata~xdata) results <-summary(model) intercept = summary(model)$coefficients[1] slope = summary(model)$coefficients[2] fstat = summary(model)$fstatistic[1] df1 = summary(model)$fstatistic[2] df2 = summary(model)$fstatistic[3] pval = summary(model)$coefficients[8] rsqd = summary(model)$r.squared rsqdadj = summary(model)$adj.r.squared lmresults = cbind(yname, xname, intercept, slope, fstat, df1, df2,pval, rsqd, rsqdadj) overallresults [i-beg+1,] = lmresults plot(ydata~xdata, sub = pval) abline(model) text(xname, yname, categ, cex=0.5, pos=3, col="blue") } overallresults #make it easier to copy and paste in Excel write.table(overallresults, "clipboard", sep="\t") dataset = mergeSPECIES #logPCSA (17), logSL (18) overallresults=matrix(nrow=11, ncol=10) j = 18 #y variable of interest (13:16) beg = 21 end= 31 for(i in beg:end){ xname = colnames(dataset)[i] yname =colnames(dataset)[j] xdata = dataset[[i]] ydata = dataset[[j]] categ = factor(dataset[[2]]) model<-lm(ydata~xdata) results <-summary(model) intercept = summary(model)$coefficients[1] slope = summary(model)$coefficients[2] fstat = summary(model)$fstatistic[1] df1 = summary(model)$fstatistic[2] df2 = summary(model)$fstatistic[3] pval = summary(model)$coefficients[8] rsqd = summary(model)$r.squared rsqdadj = summary(model)$adj.r.squared lmresults = cbind(yname, xname, intercept, slope, fstat, df1, df2,pval, rsqd, rsqdadj) overallresults [i-beg+1,] = lmresults plot(ydata~xdata, sub = pval) abline(model) text(xname, yname, categ, cex=0.5, pos=3, col="blue") } overallresults #make it easier to copy and paste in Excel write.table(overallresults, "clipboard", sep="\t") ############################################################################### ## set up phylogeny using Porter et al 2009 JEB, added a few taxa see lab notebookJan 22, 2013 musclespeciesALL <- data.frame(musclespeciesALL) read.nexus("FINALTREEFORANALYSES.nex") -> tr #check to make sure that tree and data names match # matching <-name.check(tr, musclespeciesALL, data.names=musclespeciesALL$Species) # matching$Tree.not.data # matching$Data.not.tree # row.names(musclespeciesALL)<-musclespeciesALL$Species #need to root tree for analyses tr3 <-root(tr,"Haustrali") plot(tr3) tr3 <- multi2di(tr3) #makes tree dichotomousotherwise there is a polytomy at the root #combine phylog and data into one object pods <- comparative.data(phy=tr3, data = musclespeciesALL, names.col = Species, vcv=TRUE, na.omit = FALSE, warn.dropped = TRUE) #show any dropped species pods$dropped$tips #should be some dropped here pods$dropped$unmatched.rows #should be nothing dropped here ############################################################################### #########Run phylogenetic tests #look at correlation between muscles, sizes dataset = musclespeciesALL #logPCSA (17), logSL (18) overallresults=matrix(nrow=2, ncol=10) j = 14 #y variable of interest (14:18) beg = 11 end= 13 for(i in beg:end){ xname = colnames(dataset)[i] yname =colnames(dataset)[j] xdata = dataset[[i]] ydata = dataset[[j]] categ = factor(dataset[[2]]) model.pgls <- pgls(ydata ~ xdata, data=pods, lambda = 'ML') summary(model.pgls) par(ask=T) #put stats results together into a matrix summary(model.pgls)$param.C$lambda$opt results <- cbind(xname, yname, summary(model.pgls)$coef, summary(model.pgls)$df, summary(model.pgls)$param.C$lambda$opt, summary(model.pgls)$param.C$lambda$bounds.val, summary(model.pgls)$param.C$lambda$bounds.p) results overallresults<- rbind(overallresults, results) } #delete the two first rows of NA's from table overallresults <- overallresults[-1,] overallresults <- overallresults[-1,] overallresults #make it easier to copy and paste in Excel write.table(overallresults, "clipboard", sep="\t") ############################################# overallresults=matrix(nrow=3, ncol=10) #look at correlation between appendage type and the same as above dataset = musclespeciesALL #code with non-spearer v. spearer spnosp <- musclespeciesALL$undif0spearer1smash2 #change hemisquilla to non spearer (i.e., change from 0 to 2) spnosp[spnosp == 0] <- 2 #code with non-hammer v. hammer hamnoham <- musclespeciesALL$undif0spearer1smash2 #change hemisquilla to non hammer (i.e., change from 0 to 1) hamnoham[hamnoham == 0] <- 1 #logPCSA (17), logSL (18) j = 18 #y variable of interest (14:18) beg = 17 end= 17 for(i in beg:end){ xname = colnames(dataset)[i] yname =colnames(dataset)[j] xdata = dataset[[i]] ydata = dataset[[j]] categ = factor(dataset[[2]]) model.pgls <- pgls(ydata ~ xdata*spnosp, data=pods, lambda = 'ML') summary(model.pgls) par(ask=T) #put stats results together into a matrix summary(model.pgls)$param.C$lambda$opt results <- cbind(xname, yname, summary(model.pgls)$coef, summary(model.pgls)$df, summary(model.pgls)$param.C$lambda$opt, summary(model.pgls)$param.C$lambda$bounds.val, summary(model.pgls)$param.C$lambda$bounds.p) results overallresults<- rbind(overallresults, results) } #delete the two first rows of NA's from table overallresults #make it easier to copy and paste in Excel write.table(overallresults, "clipboard", sep="\t") ####################################### #now to pgls on just appendage type and quant variable (not size) overallresults=matrix(nrow=3, ncol=9) #look at correlation between appendage type and the same as above dataset = musclespeciesALL #code with non-spearer v. spearer spnosp <- musclespeciesALL$undif0spearer1smash2 #change hemisquilla to non spearer (i.e., change from 0 to 2) spnosp[spnosp == 0] <- 2 #code with non-hammer v. hammer hamnoham <- musclespeciesALL$undif0spearer1smash2 #change hemisquilla to non hammer (i.e., change from 0 to 1) hamnoham[hamnoham == 0] <- 1 beg = 11 end= 18 for(j in beg:end){ xname =colnames(dataset)[j] xdata = dataset[[j]] model.pgls <- pgls(hamnoham ~ xdata, data=pods, lambda = 'ML') summary(model.pgls) par(ask=T) #put stats results together into a matrix summary(model.pgls)$param.C$lambda$opt results <- cbind(xname, summary(model.pgls)$coef, summary(model.pgls)$df, summary(model.pgls)$param.C$lambda$opt, summary(model.pgls)$param.C$lambda$bounds.val, summary(model.pgls)$param.C$lambda$bounds.p) results overallresults<- rbind(overallresults, results) } #delete the two first rows of NA's from table overallresults #make it easier to copy and paste in Excel write.table(overallresults, "clipboard", sep="\t") ######## #Now do pgls on data subsetted by appendage type to see if phylogenetic signal disappears #set up datasets to look at subsets of appendage types spearers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==1) smashers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==2) nonspeardata <- musclespeciesALL nonspeardata$undif0spearer1smash2[nonspeardata$undif0spearer1smash2 == 0] <- 2 #change californiensis to smasher code nonspearers<- subset(nonspeardata, nonspeardata$undif0spearer1smash2==2) nonhammerdata <- musclespeciesALL nonhammerdata$undif0spearer1smash2[nonhammerdata$undif0spearer1smash2 == 0] <- 1 #change californiensis to spearer code nonhammers <- subset(nonhammerdata, nonhammerdata$undif0spearer1smash2==1) #need to check for normality #look at distributions for phylogeny testing - test for normality, etc. #make a loop to get all the p-values for a Shapiro Wilk test for each species specname = nonhammers sizem = ncol(specname) shap=vector(length=sizem) i=3 #loop through size of dataset and enter result into the appropriate column for (i in 3:(sizem-1)) { res = shapiro.test(specname[[i]]) shap[i] = res$p.value } shap write.table(shap, "clipboard", sep="\t") #now run pgls on these subset datasets dataset = data.frame(nonhammers) dataset =droplevels(dataset) ) Species = factor(dataset$Species) #match with tree #combine phylog and data into one object categpods <- comparative.data(phy=tr3, data = dataset, names.col = Species, vcv=TRUE, na.omit = FALSE, warn.dropped = TRUE) #show any dropped species pods$dropped$tips #should be some dropped here pods$dropped$unmatched.rows #should be nothing dropped here overallresults=matrix(nrow=3, ncol=10) j = 18 #y variable of interest (14:18) (14 is log merus length) beg = 11 #xdata: 11:14 end= 14 for(i in beg:end){ xname = colnames(dataset)[i] yname =colnames(dataset)[j] xdata = dataset[[i]] ydata = dataset[[j]] model.pgls <- pgls(ydata ~ xdata, data=categpods, lambda = 'ML') summary(model.pgls) par(ask=T) #put stats results together into a matrix summary(model.pgls)$param.C$lambda$opt results <- cbind(xname, yname, summary(model.pgls)$coef, summary(model.pgls)$df, summary(model.pgls)$param.C$lambda$opt, summary(model.pgls)$param.C$lambda$bounds.val, summary(model.pgls)$param.C$lambda$bounds.p) results overallresults<- rbind(overallresults, results) } overallresults #make it easier to copy and paste in Excel write.table(overallresults, "clipboard", sep="\t") ######## #correlation between spring parameters and muscle parameters using pgls #use mergeSPECIES dataset dataset = mergeSPECIES #match with tree read.nexus("FINALTREEFORANALYSES.nex") -> tr #check to make sure that tree and data names match # matching <-name.check(tr, musclespeciesALL, data.names=musclespeciesALL$Species) # matching$Tree.not.data # matching$Data.not.tree # row.names(musclespeciesALL)<-musclespeciesALL$Species #need to root tree for analyses tr3 <-root(tr,"Haustrali") plot(tr3) tr3 <- multi2di(tr3) #makes tree dichotomousotherwise there is a polytomy at the root #combine phylog and data into one object springpods <- comparative.data(phy=tr3, data = dataset, names.col = Species, vcv=TRUE, na.omit = FALSE, warn.dropped = TRUE) #show any dropped species springpods$dropped$tips #should be some dropped here springpods$dropped$unmatched.rows #should be nothing dropped here overallresults=matrix(nrow=3, ncol=10) j = 18 #x variable of interest 17:18 beg = 27 #ydata: (27:31) end= 29 for(i in beg:end){ xname = colnames(dataset)[j] yname =colnames(dataset)[i] xdata = dataset[[j]] ydata = dataset[[i]] model.pgls <- pgls(ydata ~ xdata, data=springpods, lambda = 'ML') summary(model.pgls) par(ask=T) #put stats results together into a matrix summary(model.pgls)$param.C$lambda$opt results <- cbind(xname, yname, summary(model.pgls)$coef, summary(model.pgls)$df, summary(model.pgls)$param.C$lambda$opt, summary(model.pgls)$param.C$lambda$bounds.val, summary(model.pgls)$param.C$lambda$bounds.p) results overallresults<- rbind(overallresults, results) } overallresults #make it easier to copy and paste in Excel write.table(overallresults, "clipboard", sep="\t") #################################### #now include appendage type as factor in muscle-spring analysis spnosp <- mergeSPECIES$undif0spearer1smash2 #change hemisquilla to non spearer (i.e., change from 0 to 2) spnosp[spnosp == 0] <- 2 #code with non-hammer v. hammer hamnoham <- mergeSPECIES$undif0spearer1smash2 #change hemisquilla to non hammer (i.e., change from 0 to 1) hamnoham[hamnoham == 0] <- 1 overallresults=matrix(nrow=3, ncol=10) j = 18 #x variable of interest 17:18 beg = 27 #ydata: (27:31) end= 29 for(i in beg:end){ xname = colnames(dataset)[j] yname =colnames(dataset)[i] xdata = dataset[[j]] ydata = dataset[[i]] model.pgls <- pgls(ydata ~ xdata*spnosp, data=springpods, lambda = 'ML') summary(model.pgls) par(ask=T) #put stats results together into a matrix summary(model.pgls)$param.C$lambda$opt results <- cbind(xname, yname, summary(model.pgls)$coef, summary(model.pgls)$df, summary(model.pgls)$param.C$lambda$opt, summary(model.pgls)$param.C$lambda$bounds.val, summary(model.pgls)$param.C$lambda$bounds.p) results overallresults<- rbind(overallresults, results) } overallresults #make it easier to copy and paste in Excel write.table(overallresults, "clipboard", sep="\t") #need to look at phylogenetically corrected data with slopes overlaid #can't just plot raw species data this time, because of phylogenetic signal #plot independent contrasts with regression line overlaid c <-qplot(xdata, ydata, data=mergeSPECIES, colour= AppenType, shape= AppenType, geom = c("point", "smooth"), method = "lm", ) #put the data on the graph plots <- c + labs(x=xname, y=yname) +theme(axis.title.x=element_text(size=20, face = "bold",family = "sans"), axis.title.y=element_text(size=20, face = "bold",family = "sans"), abline(model.pgls) axis.text.x = element_text(angle=90, vjust=0.5, size=16), axis.text.y = element_text(angle=90, vjust=0.5, size=16)) #adjust the appearance print (plots) #make it show up on the screen #********************************************************** #MAKE FIGURES FOR PAPER #within species scaling ######### # #linear regressions within species: # falcatusALL # chiragraALL # ciliataALL # smithiiALL #use unlogged data here, because the unlogged data are normally distributed, whereas logged are not dataset = data.frame(ciliataALL) overallresults=matrix(nrow=4, ncol=10) j = 15 #y variable of interest (15 is PCSA, 16 is SL) ylab = "PCSA (mm^2) i= 12#10 is cpl, 12 is merus len, xlab = "Carapace Length (mm)", xlab = "Merus Length (mm)" xdata = dataset[[i]] ydata = dataset[[j]] model = lm(ydata~xdata) dev.off() par(ask = FALSE) plot(ydata~xdata, xlab = "Merus Length (mm)", ylab = "PCSA (mm^2)",pch=19, cex=1.5, cex.lab=1.5, cex.axis = 1.5, xlim =c(6, 14), ylim=c(5, 25)) abline(model) dataset = data.frame(falcatusALL) overallresults=matrix(nrow=4, ncol=10) j = 15 #y variable of interest (15 is PCSA, 16 is SL) ylab = "PCSA (mm^2) i= 12#10 is cpl, 12 is merus len, xlab = "Carapace Length (mm)", xlab = "Merus Length (mm)" xdata = dataset[[i]] ydata = dataset[[j]] model = lm(ydata~xdata) dev.off() par(ask = FALSE) plot(ydata~xdata, xlab = "Merus Length (mm)", ylab = "PCSA (mm^2)",pch=19, cex=1.5, cex.lab=1.5, cex.axis = 1.5, xlim =c(6, 14), ylim=c(5, 25)) abline(model) dataset = data.frame(ciliataALL) overallresults=matrix(nrow=4, ncol=10) j = 12 #y variable of interest (15 is PCSA, 16 is SL) ylab = "PCSA (mm^2) i= 10#10 is cpl, 12 is merus len, xlab = "Carapace Length (mm)", xlab = "Merus Length (mm)" xdata = dataset[[i]] ydata = dataset[[j]] model = lm(ydata~xdata) dev.off() par(ask = FALSE) plot(ydata~xdata, xlab = "Carapace Length (mm)", ylab = "Merus Length (mm)",pch=19, cex=1.5, cex.lab=1.5, cex.axis = 1.5, ylim =c(6, 14), xlim =c(7, 17)) abline(model) dataset = data.frame(falcatusALL) overallresults=matrix(nrow=4, ncol=10) j = 12 #y variable of interest (15 is PCSA, 16 is SL) ylab = "PCSA (mm^2) i= 10#10 is cpl, 12 is merus len, xlab = "Carapace Length (mm)", xlab = "Merus Length (mm)" xdata = dataset[[i]] ydata = dataset[[j]] model = lm(ydata~xdata) dev.off() par(ask = FALSE) plot(ydata~xdata, xlab = "Carapace Length (mm)", ylab = "Merus Length (mm)",pch=19, cex=1.5, cex.lab=1.5, cex.axis = 1.5, ylim =c(6, 14), xlim =c(7, 17)) abline(model) #now cross species stuff #want: PCSA v. ML, SL v. ML, ML v. CPL, SPR force v. PCSA #let's do this in an efficient way using ggplot LOVE IT! #add a column to make continuous appendage type into categorical #or can just say factor(variable) and will make it categorical musclespeciesALL["AppenType"] <- NA musclespeciesALL$AppenType[musclespeciesALL$undif0spearer1smash2==0]<-"U" musclespeciesALL$AppenType[musclespeciesALL$undif0spearer1smash2==1]<-"P" musclespeciesALL$AppenType[musclespeciesALL$undif0spearer1smash2==2]<-"M" j = 10 #y""Sarcomere Length (micron)" i=6 #x 6 "Merus Length (mm)" xvar = musclespeciesALL[[i]] yvar = musclespeciesALL[[j]] par(mar=c(4.5, 5, .5, 0.5)) #set margins around graph plot(yvar~xvar, col="white", xlab = "Merus Length (mm)", ylab = expression(paste("Sarcomere Length (",mu, " m)")), cex.lab=1.5, cex.axis = 1.5) spearers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==1) smashers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==2) undif <- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==0) points(spearers$MeanmeanSL~spearers$MeanAvgmeruslengthmm, col = 4, lwd=3,pch=15, cex=2) points(smashers$MeanmeanSL~smashers$MeanAvgmeruslengthmm, col = 25, lwd=3,pch=16, cex=2) points(undif$MeanmeanSL~undif$MeanAvgmeruslengthmm, col = 6, lwd=3,pch=17, cex=2) j = 9 #y" 9 PCSA i=6 #x 6 "Merus Length (mm)" xvar = musclespeciesALL[[i]] yvar = musclespeciesALL[[j]] par(mar=c(4.5, 5, .5, 0.5)) #set margins around graph plot(yvar~xvar, col="white", xlab = "Merus Length (mm)", ylab = parse(text = "PCSA (mm^2)"), cex.lab=1.5, cex.axis = 1.5) spearers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==1) smashers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==2) undif <- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==0) points(spearers$MeanPCSAisAsin2amm2~spearers$MeanAvgmeruslengthmm, col = 4, lwd=3,pch=15, cex=2) points(smashers$MeanPCSAisAsin2amm2~smashers$MeanAvgmeruslengthmm, col = 25, lwd=3,pch=16, cex=2) points(undif$MeanPCSAisAsin2amm2~undif$MeanAvgmeruslengthmm, col = 6, lwd=3,pch=17, cex=2) j = 6 #6 "Merus Length (mm)" i=4 #x 4 "Carapace Length (mm)" xvar = musclespeciesALL[[i]] yvar = musclespeciesALL[[j]] par(mar=c(4.5, 5, .5, 0.5)) #set margins around graph plot(yvar~xvar, col="white", xlab = "Merus Length (mm)", ylab = "Carapace Length (mm)", cex.lab=1.5, cex.axis = 1.5) spearers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==1) smashers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==2) undif <- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==0) points(spearers$MeanAvgmeruslengthmm~spearers$Meanavgcarapacelengthmm, col = 4, lwd=3,pch=15, cex=2) points(smashers$MeanAvgmeruslengthmm~smashers$Meanavgcarapacelengthmm, col = 25, lwd=3,pch=16, cex=2) points(undif$MeanAvgmeruslengthmm~undif$Meanavgcarapacelengthmm, col = 6, lwd=3,pch=17, cex=2) j = 24 #maxforce i=9 #x MeanPCSAisAsin2amm2 xvar = mergeSPECIES[[i]] yvar = mergeSPECIES[[j]] par(mar=c(4.5, 5, .5, 0.5)) #set margins around graph plot(yvar~xvar, col="white", xlab = parse(text = "PCSA (mm^2)"), ylab = "Maximum Spring Force (N)", cex.lab=1.5, cex.axis = 1.5) spearers<- subset(mergeSPECIES, mergeSPECIES$undif0spearer1smash2==1) smashers<- subset(mergeSPECIES, mergeSPECIES$undif0spearer1smash2==2) undif <- subset(mergeSPECIES, mergeSPECIES$undif0spearer1smash2==0) points(spearers$maxforce~spearers$MeanPCSAisAsin2amm2, col = 4, lwd=3,pch=15, cex=2) points(smashers$maxforce~smashers$MeanPCSAisAsin2amm2, col = 25, lwd=3,pch=16, cex=2) points(undif$maxforce~undif$MeanPCSAisAsin2amm2, col = 6, lwd=3,pch=17, cex=2) #### LOGGED j = 18 #yLOGSarcomere Length (micron)" i=14 #x 6 LOG "Merus Length (mm)" xvar = musclespeciesALL[[i]] yvar = musclespeciesALL[[j]] par(mar=c(4.5, 5, .5, 0.5)) #set margins around graph plot(yvar~xvar, col="white", xlab = "Log Merus Length (mm)", ylab = expression(paste("Log Sarcomere Length (",mu, " m)")), cex.lab=1.5, cex.axis = 1.5) spearers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==1) smashers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==2) undif <- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==0) points(spearers$Meanlogsl~spearers$Meanlogmeruslengthmm, col = 4, lwd=3,pch=15, cex=2) points(smashers$Meanlogsl~smashers$Meanlogmeruslengthmm, col = 25, lwd=3,pch=16, cex=2) points(undif$Meanlogsl~undif$Meanlogmeruslengthmm, col = 6, lwd=3,pch=17, cex=2) j = 10 #y Sarcomere Length (micron)" i=6 #x 6 "Merus Length (mm)" xvar = musclespeciesALL[[i]] yvar = musclespeciesALL[[j]] par(mar=c(4.5, 5, .5, 0.5)) #set margins around graph plot(yvar~xvar, col="white", xlab = "Merus Length (mm)", ylab = expression(paste("Sarcomere Length (",mu, " m)")), cex.lab=1.5, cex.axis = 1.5) spearers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==1) smashers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==2) undif <- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==0) points(spearers$MeanmeanSL~spearers$MeanAvgmeruslengthmm, col = 4, lwd=3,pch=15, cex=2) points(smashers$MeanmeanSL~smashers$MeanAvgmeruslengthmm, col = 25, lwd=3,pch=16, cex=2) points(undif$MeanmeanSL~undif$MeanAvgmeruslengthmm, col = 6, lwd=3,pch=17, cex=2) j = 17 #y" LOGPCSA i=14 #x LOG Merus Length (mm)" xvar = musclespeciesALL[[i]] yvar = musclespeciesALL[[j]] par(mar=c(4.5, 5, .8, 0.5)) #set margins around graph plot(yvar~xvar, col="white", xlab = "Log Merus Length (mm)", ylab = "Log PCSA (mm^2)", cex.lab=1.5, cex.axis = 1.5) spearers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==1) smashers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==2) undif <- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==0) points(spearers$Meanlogpcsa~spearers$Meanlogmeruslengthmm, col = 4, lwd=3,pch=15, cex=2) points(smashers$Meanlogpcsa~smashers$Meanlogmeruslengthmm, col = 25, lwd=3,pch=16, cex=2) points(undif$Meanlogpcsa~undif$Meanlogmeruslengthmm, col = 6, lwd=3,pch=17, cex=2) pods <- comparative.data(phy=tr3, data = musclespeciesALL, names.col = Species, vcv=TRUE, na.omit = FALSE, warn.dropped = TRUE) #show any dropped species pods$dropped$tips #should be some dropped here pods$dropped$unmatched.rows #should be nothing dropped here model.pgls <- pgls(yvar ~ xvar, data=pods, lambda = 'ML') abline(model.pgls) j = 14 #LOG Merus Length (mm)" i=12 #x "LOG Carapace Length (mm)" xvar = musclespeciesALL[[i]] yvar = musclespeciesALL[[j]] par(mar=c(4.5, 5, .5, 0.5)) #set margins around graph plot(yvar~xvar, col="white", xlab = "Log Carapace Length (mm)", ylab = "Log Merus Length (mm)", cex.lab=1.5, cex.axis = 1.5) spearers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==1) smashers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==2) undif <- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==0) points(spearers$Meanlogmeruslengthmm~spearers$Meanlogavgcplmm, col = 4, lwd=3,pch=15, cex=2) points(smashers$Meanlogmeruslengthmm~smashers$Meanlogavgcplmm, col = 25, lwd=3,pch=16, cex=2) points(undif$Meanlogmeruslengthmm~undif$Meanlogavgcplmm, col = 6, lwd=3,pch=17, cex=2) pods <- comparative.data(phy=tr3, data = musclespeciesALL, names.col = Species, vcv=TRUE, na.omit = FALSE, warn.dropped = TRUE) #show any dropped species pods$dropped$tips #should be some dropped here pods$dropped$unmatched.rows #should be nothing dropped here model.pgls <- pgls(yvar ~ xvar, data=pods, lambda = 'ML') abline(model.pgls) j = 29 #logmaxforce i=17 #x Meanlogpcsa xvar = mergeSPECIES[[i]] yvar = mergeSPECIES[[j]] par(mar=c(4.5, 5, .5, 0.5)) #set margins around graph plot(yvar~xvar, col="white", xlab = "Log PCSA (mm^2)", ylab = "Log Maximum Spring Force (N)", cex.lab=1.5, cex.axis = 1.5) spearers<- subset(mergeSPECIES, mergeSPECIES$undif0spearer1smash2==1) smashers<- subset(mergeSPECIES, mergeSPECIES$undif0spearer1smash2==2) undif <- subset(mergeSPECIES, mergeSPECIES$undif0spearer1smash2==0) points(spearers$logmaxforce~spearers$Meanlogpcsa, col = 4, lwd=3,pch=15, cex=2) points(smashers$logmaxforce~smashers$Meanlogpcsa, col = 25, lwd=3,pch=16, cex=2) points(undif$logmaxforce~undif$Meanlogpcsa, col = 6, lwd=3,pch=17, cex=2) read.nexus("FINALTREEFORANALYSES.nex") -> tr tr3 <-root(tr,"Haustrali") tr3 <- multi2di(tr3) #makes tree dichotomousotherwise there is a polytomy at the root #combine phylog and data into one object springpods <- comparative.data(phy=tr3, data = mergeSPECIES, names.col = Species, vcv=TRUE, na.omit = FALSE, warn.dropped = TRUE) #show any dropped species springpods$dropped$tips #should be some dropped here springpods$dropped$unmatched.rows #should be nothing dropped here model.pgls <- pgls(yvar ~ xvar, data=springpods, lambda = 'ML') abline(model.pgls) #### data not logged, but plotted on log scale ## Figure 5 on Evolution muscle paper j = 10 #[10] "MeanmeanSL" i=6 #"MeanAvgmeruslengthmm" xvar = musclespeciesALL[[i]] yvar = musclespeciesALL[[j]] par(mar=c(4.5, 5, .5, 0.5)) #set margins around graph plot(yvar~xvar, xaxt="n", yaxt = "n", log = "xy",col="white", xlab = "Merus Length (mm)", ylab = expression(paste("Sarcomere Length (",mu, " m)")), cex.lab=1.5) #run after the plot function closes axis(1, cex.axis = 1.5,at=c(-10, 5, 10, 20, 30, 50), labels=c("0", "5", "10", "20", "30", "40")) axis(2, cex.axis = 1.5) spearers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==1) smashers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==2) undif <- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==0) points(spearers$MeanmeanSL~spearers$MeanAvgmeruslengthmm, col = 4, lwd=3,pch=15, cex=2) points(smashers$MeanmeanSL~smashers$MeanAvgmeruslengthmm, col = 25, lwd=3,pch=16, cex=2) points(undif$MeanmeanSL~undif$MeanAvgmeruslengthmm, col = 6, lwd=3,pch=17, cex=2) j = 9 #y" MeanPCSAisAsin2amm2 i=6 #xMeanAvgmeruslengthmm xvar = musclespeciesALL[[i]] yvar = musclespeciesALL[[j]] par(mar=c(4.5, 5, .8, 0.5)) #set margins around graph plot(yvar~xvar, xaxt="n", yaxt = "n", log = "xy",col="white", xlab = "Merus Length (mm)", ylab = expression(paste("PCSA (mm^2)")), cex.lab=1.5) #run after the plot function closes axis(1, cex.axis = 1.5,at=c(-10, 5, 10, 20, 30, 50), labels=c("-10", "5", "10", "20", "30", "50")) axis(2, cex.axis = 1.5,at=c(-1, 5, 10, 20, 30, 70), labels=c("-1", "5", "10", "20", "30", "70")) spearers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==1) smashers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==2) undif <- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==0) points(spearers$MeanPCSAisAsin2amm2~spearers$MeanAvgmeruslengthmm, col = 4, lwd=3,pch=15, cex=2) points(smashers$MeanPCSAisAsin2amm2~smashers$MeanAvgmeruslengthmm, col = 25, lwd=3,pch=16, cex=2) points(undif$MeanPCSAisAsin2amm2~undif$MeanAvgmeruslengthmm, col = 6, lwd=3,pch=17, cex=2) pods <- comparative.data(phy=tr3, data = musclespeciesALL, names.col = Species, vcv=TRUE, na.omit = FALSE, warn.dropped = TRUE) #show any dropped species pods$dropped$tips #should be some dropped here pods$dropped$unmatched.rows #should be nothing dropped here j = 17 #y" LOGPCSA i=14 #x LOG Merus Length (mm)" xvar = musclespeciesALL[[i]] yvar = musclespeciesALL[[j]] model.pgls <- pgls(yvar ~ xvar, data=pods, lambda = 'ML') abline(model.pgls) ############## j = 6 #xMeanAvgmeruslengthmm i=4 #x "Meanavgcarapacelengthmm" xvar = musclespeciesALL[[i]] yvar = musclespeciesALL[[j]] par(mar=c(4.5, 5, .5, 0.5)) #set margins around graph plot(yvar~xvar, col="white", xlab = "Carapace Length (mm)", log = "xy", yaxt = "n", ylab = "Merus Length (mm)", cex.lab=1.5, cex.axis = 1.5) axis(2, cex.axis = 1.5, at=c(0,5, 10, 20, 30), labels=c("0", "5", "10", "20", "30")) spearers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==1) smashers<- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==2) undif <- subset(musclespeciesALL, musclespeciesALL$undif0spearer1smash2==0) points(spearers$MeanAvgmeruslengthmm~spearers$Meanavgcarapacelengthmm, col = 4, lwd=3,pch=15, cex=2) points(smashers$MeanAvgmeruslengthmm~smashers$Meanavgcarapacelengthmm, col = 25, lwd=3,pch=16, cex=2) points(undif$MeanAvgmeruslengthmm~undif$Meanavgcarapacelengthmm, col = 6, lwd=3,pch=17, cex=2) j = 14 #LOG Merus Length (mm)" i=12 #x "LOG Carapace Length (mm)" xvar = musclespeciesALL[[i]] yvar = musclespeciesALL[[j]] pods <- comparative.data(phy=tr3, data = musclespeciesALL, names.col = Species, vcv=TRUE, na.omit = FALSE, warn.dropped = TRUE) #show any dropped species pods$dropped$tips #should be some dropped here pods$dropped$unmatched.rows #should be nothing dropped here model.pgls <- pgls(yvar ~ xvar, data=pods, lambda = 'ML') abline(model.pgls) ############ j = 24 #maxforce i=9 #x MeanPCSAisAsin2amm2 xvar = mergeSPECIES[[i]] yvar = mergeSPECIES[[j]] par(mar=c(4.5, 5, .5, 0.5)) #set margins around graph plot(yvar~xvar, col="white", xlab = "PCSA (mm^2)", log = "xy", ylab = "Maximum Spring Force (N)", cex.lab=1.5, cex.axis = 1.5) spearers<- subset(mergeSPECIES, mergeSPECIES$undif0spearer1smash2==1) smashers<- subset(mergeSPECIES, mergeSPECIES$undif0spearer1smash2==2) undif <- subset(mergeSPECIES, mergeSPECIES$undif0spearer1smash2==0) points(spearers$maxforce~spearers$MeanPCSAisAsin2amm2, col = 4, lwd=3,pch=15, cex=2) points(smashers$maxforce~smashers$MeanPCSAisAsin2amm2, col = 25, lwd=3,pch=16, cex=2) points(undif$maxforce~undif$MeanPCSAisAsin2amm2, col = 6, lwd=3,pch=17, cex=2) j = 29 #logmaxforce i=17 #x Meanlogpcsa xvar = mergeSPECIES[[i]] yvar = mergeSPECIES[[j]] read.nexus("FINALTREEFORANALYSES.nex") -> tr tr3 <-root(tr,"Haustrali") tr3 <- multi2di(tr3) #makes tree dichotomousotherwise there is a polytomy at the root #combine phylog and data into one object springpods <- comparative.data(phy=tr3, data = mergeSPECIES, names.col = Species, vcv=TRUE, na.omit = FALSE, warn.dropped = TRUE) #show any dropped species springpods$dropped$tips #should be some dropped here springpods$dropped$unmatched.rows #should be nothing dropped here model.pgls <- pgls(yvar ~ xvar, data=springpods, lambda = 'ML') abline(model.pgls) #set defaults for graphs theme_set(theme_bw(24)) j = 10 #y""Sarcomere Length (micron)" i=4 #x "Carapace Length (mm)" xdata = musclespeciesALL[[i]] ydata = musclespeciesALL[[j]] c <-qplot(xdata, ydata, data=musclespeciesALL, colour= AppenType, shape= AppenType, size = 12), method = "lm", ) #put the data on the graph muscleplots <- c + labs(x="Carapace Length (mm)", y="Sarcomere Length (micron)") +theme(axis.title.x=element_text(size=20, face = "bold",family = "sans"), axis.title.y=element_text(size=20, face = "bold",family = "sans"), axis.text.x = element_text(angle=90, vjust=0.5, size=16), axis.text.y = element_text(angle=90, vjust=0.5, size=16), panel.grid.minor = element_blank(), panel.grid.major = element_blank()) print (muscleplots) #make it show up on the screen #, geom = c("point", "smooth" #adjust the appearance #this is using ggplot #pdf(file = "loggedspeciesgraphsv4.pdf", width = 6, height = 6) #now do the same thing logged j = 18 #y variable of interest (7:10 unlogged variables, 15:18 logged) for(i in 16:18){ xname = colnames(musclespeciesALL)[i] yname =colnames(musclespeciesALL)[j] xdata = musclespeciesALL[[i]] ydata = musclespeciesALL[[j]] p <- ggplot(musclespeciesALL, aes(xdata, ydata,colour= AppenType, shape= AppenType))+geom_point(aes(size=3))+ geom_smooth(method="lm") muscleplots <- p + labs(x=xname, y=yname) +theme(axis.title.x=element_text(size=20, face = "bold",family = "sans"), axis.title.y=element_text(size=20, face = "bold",family = "sans"), axis.text.x = element_text(angle=90, vjust=0.5, size=16), axis.text.y = element_text(angle=90, vjust=0.5, size=16) ) #adjust the appearance print(muscleplots) } #dev.off() #merus scaling #pdf(file = "loggedmerusscaling.pdf", width = 6, height = 6) j = 14 #y variable of interest (7:10 unlogged variables, 15:18 logged) for(i in 11:13){ xname = colnames(musclespeciesALL)[i] yname =colnames(musclespeciesALL)[j] xdata = musclespeciesALL[[i]] ydata = musclespeciesALL[[j]] p <- ggplot(musclespeciesALL, aes(xdata, ydata,colour= AppenType, shape= AppenType))+geom_point(aes(size=2))+ geom_smooth(method="lm") muscleplots <- p + labs(x=xname, y=yname) +theme(axis.title.x=element_text(size=20, face = "bold",family = "sans"), axis.title.y=element_text(size=20, face = "bold",family = "sans"), axis.text.x = element_text(angle=90, vjust=0.5, size=16), axis.text.y = element_text(angle=90, vjust=0.5, size=16) ) print(muscleplots) } #dev.off() #plot muscle and springdata #mergeSPECIES provides merged averaged species data #mergeINDIVID provides data for which individuals are matched #pdf(file = "musclespringgraphs2nd.pdf", width = 6, height = 6) AppenType = factor(mergeSPECIES$undif0spearer1smash2) j = 18 #y variable of interest (17:18 logged PCSA and SL) for(i in 27:31){ xname = colnames(mergeSPECIES)[i] yname =colnames(mergeSPECIES)[j] xdata = mergeSPECIES[[i]] ydata = mergeSPECIES[[j]] p <- ggplot(mergeSPECIES, aes(xdata, ydata,colour= AppenType, shape= AppenType))+geom_point(aes(size=2))+ geom_smooth(method="lm") muscleplots <- p + labs(x=xname, y=yname) +theme(axis.title.x=element_text(size=20, face = "bold",family = "sans"), axis.title.y=element_text(size=20, face = "bold",family = "sans"), axis.text.x = element_text(angle=90, vjust=0.5, size=16), axis.text.y = element_text(angle=90, vjust=0.5, size=16) ) print(muscleplots) } #dev.off() log = "xy",