###Reading the database for working base<-read.csv(file="PhaetornisHeliconia2.csv", header=TRUE) ####Correlations Phaethornis - Heliconia________________________________########### base<-read.csv(file="PhaetornisHeliconia2.csv", header=TRUE) names(base) summary(base) boxplot(base$Bill length~base$OGU, col="blue") boxplot(base$Bill curvature~base$OGU, col="blue") boxplot(base$Bill width~base$OGU, col="blue") ####Working first on OGU 1 (Nayarit) OGU1<-subset(base,OGU==1) summary(OGU1) dev.off() #### First, Bill length boxplot(OGU1$Bill length~OGU1$Species, main="Length OGU 1", ylab="mm",cex.labs=2) #Testing for normality hist(OGU1$Bill length,col="red") #Expecting normality?? shapiro.test(OGU1$Bill length) #Result: #W = 0.95274, p-value = 0.004024 #Not normally distributed ##Testing for significant differences wi.Bill length.OGU1<-wilcox.test(Bill length~Species,data=OGU1) wi.Bill length.OGU1 #Result: #Wilcoxon rank sum test with continuity correction #data: Bill length by Species #W = 690, p-value = 0.377 #alternative hypothesis: true location shift is not equal to #No significant differences !!!!!! ####Now, Bill curvature boxplot(OGU1$Bill curvature~OGU1$Species, main="Curvature OGU1", ylab="mm",cex.labs=2) #Testing for normality hist(OGU1$Bill curvature,col="red") #Expecting normality?? shapiro.test(OGU1$Bill curvature) #Result is: #W = 0.9562, p-value = 0.01266 #Not normally distributed ##Testing for significant differences wi.curva.OGU1<-wilcox.test(Bill curvature~Species,data=OGU1) wi.curva.OGU1 #Result is: #Wilcoxon rank sum test with continuity correction #ata: Bill curvature by Species #W = 927, p-value = 2.992e-05 #alternative hypothesis: true location shift is not equal to 0 #Testing for significant differences!!!!!! ###Now Bill width boxplot(OGU1$Bill width~OGU1$Species, main="Width OGU1", ylab="mm",cex.labs=2) #Testing for normality hist(OGU1$Bill width,col="red") #Expecting normality?? shapiro.test(OGU1$Bill width) #Result is: #W = 0.96867, p-value = 0.0805 #They are normally distributed ##Testing for significant differences wi.width.OGU1<-wilcox.test(OGU1$Bill width~OGU1$Species) wi.width.OGU1 #Result is: #Wilcoxon rank sum test with continuity correction #data: OGU1$Bill width by OGU1$Species #W = 871, p-value = 0.0002897 #alternative hypothesis: true location shift is not equal to 0 #Testing for significant differences!!!!!! ####Building coorelations ###Bill length summary(OGU1)#hel 54, phae 29 phae<-subset(OGU1,Species=="Phaethornis") #Only Phaethornis from OGU1 heli<-subset(OGU1,Species=="Heliconia") #Only heliconia from OGU1 heli_Bill length<-sample(heli$Bill length,29,replace=F) #Pearson Parametric Correlation using ordered values cor.test(sort(phae$Bill length),sort(heli_Bill length),method="pearson") plot(sort(phae$Bill length),sort(heli_Bill length),main="Correlations Bill length OGU 1", xlab="Bill length Phaetornis",ylab="Bill length Heliconia",pch=20) #Result is: #data: sort(phae$Bill length) and sort(heli_Bill length) #t = 19.106, df = 27, p-value < 2.2e-16 #alternative hypothesis: true correlation is not equal to 0 #95 percent confidence interval: # 0.9258980 0.9835974 #sample estimates: # cor #0.9649507 ###Bill curvature summary(OGU1)#hel 51, phae 23 summary(OGU1$Bill curvature) phae<-subset(OGU1,Species=="Phaethornis") #Only Phaethornis from OGU1 heli<-subset(OGU1,Species=="Heliconia") #Only heliconia from OGU1 heli_curva<-sample(heli$Bill curvature,23,replace=F) heli_curva #Pearson Parametric Correlation using ordered values cor.test(sort(phae$Bill curvature),sort(heli_curva),method="pearson") plot(sort(phae$Bill curvature),sort(heli_curva),main="Correlation Bill curvature OGU 1", xlab="Bill curvature Phaethornis",ylab="Bill curvature Heliconia",pch=20) #Result is: #data: sort(phae$Bill curvature) and sort(heli_curva) #t = 8.6565, df = 21, p-value = 2.273e-08 #alternative hypothesis: true correlation is not equal to 0 #95 percent confidence interval: # 0.7418605 0.9499362 #sample estimates: # cor #0.8837998 ###Bill width summary(OGU1)#hel 41, phae 28 phae<-subset(OGU1,Species=="Phaethornis")#Only Phaethornis from OGU1 phae_Bill width<-sample(phae$Bill width,24,replace=F) phae_Bill width heli<-subset(OGU1,Species=="Heliconia") #Only heliconia from OGU1 heli_Bill width<-sample(heli$Bill width,24,replace=F) heli_Bill width<-c(4.4,3.7,3.3,4.4,3.3,2.1,5.2,5.5,1.9,4.9,4.9,2.9,4.1,2.2,3.6,2.9,4.0,4.28,3.7,4.4,4.10,2.1,4.5,3.5) heli_Bill width #Pearson Parametric Correlation using ordered values cor.test(sort(phae_Bill width),sort(heli_Bill width),method="pearson") plot(sort(phae_Bill width),sort(heli_Bill width),main="Correlation Bill width OGU 1", xlab="Bill width Phaetornis",ylab="Bill width Heliconia",pch=20) #Result is: #data: sort(phae_Bill width) and sort(heli_Bill width) #t = 10.607, df = 22, p-value = 4.089e-10 #alternative hypothesis: true correlation is not equal to 0 #95 percent confidence interval: # 0.8100147 0.9627679 #sample estimates: # cor #0.9145693 ####Working on OGU 2 (Jalisco) OGU2<-subset(base,OGU==2) summary(OGU2) #### First Bill length boxplot(OGU2$Bill length~OGU2$Species, main="Length OGU 2", ylab="mm",cex.labs=2) #Testing for normality hist(OGU2$Bill length,col="red") #Expecting normality?? shapiro.test(OGU2$Bill length) #Result is: #W = 0.92958, p-value = 0.0243 #Not normally distributed ##Testing for significant differences wi.Bill length.OGU2<-wilcox.test(Bill length~Species,data=OGU2) wi.Bill length.OGU2 #Resultado is: #Wilcoxon rank sum test with continuity correction #data: Bill length by Species #W = 261, p-value = 0.0002557 #alternative hypothesis: true location shift is not equal to 0 #There are significant differences!!!!!! ####Now working on Bill curvature boxplot(OGU2$Bill curvature~OGU2$Species, main="Curvature OGU 2", ylab="mm",cex.labs=2) #Testing for normality hist(OGU2$Bill curvature,col="red") #Expecting normality?? shapiro.test(OGU2$Bill curvature) #Resultado is: #W = 0.93632, p-value = 0.05314 #Not normally distributed ##Testing for significant differences wi.curva.OGU2<-wilcox.test(Bill curvature~Species,data=OGU2) wi.curva.OGU2 #Result is: #Wilcoxon rank sum test with continuity correction #data: Bill curvature by Species #W = 104, p-value = 0.6808 #alternative hypothesis: true location shift is not equal to 0 #No significant differences!!!!!! ###Now working on Bill width boxplot(OGU2$Bill width~OGU2$Species, main="Width OGU 2", ylab="mm",cex.labs=2) #Testing normality hist(OGU2$Bill width,col="red") #Expecting normality??? shapiro.test(OGU2$Bill width) #Result is: #W = 0.95445, p-value = 0.1551 #Normally distributed ##Testing for significant differences t.test.OGU2<-t.test(Bill width~Species,data=OGU2) t.test.OGU2 #Result is: # #data: Bill width by Species #t = 1.5377, df = 32.363, p-value = 0.1338 #alternative hypothesis: true difference in means is not equal to 0 #95 percent confidence interval: # -0.09877424 0.70831530 #sample estimates: # mean in group Heliconia mean in group Phaethornis #3.347826 3.043056 #No significant differences !!!!!! ####Building correlatons ###Bill length summary(OGU2)#hel 23, phae 13 phae<-subset(OGU2,Species=="Phaethornis") #Only Ohaethornis from OGU1 heli<-subset(OGU2,Species=="Heliconia") #Only heliconia from OGU1 heli_Bill length<-sample(heli$Bill length,13,replace=F) heli_Bill length #Correlacion paremtrica de earson usando los valores ordenados cor.test(sort(phae$Bill length),sort(heli_Bill length),method="pearson") plot(sort(phae$Bill length),sort(heli_Bill length), main="Correlación Bill length OGU 2", xlab="Bill length Phaethornis", ylab="Bill length Heliconia", pch=20) #Result is: #data: sort(phae$Bill length) and sort(heli_Bill length) #t = 19.106, df = 27, p-value < 2.2e-16 #alternative hypothesis: true correlation is not equal to 0 #95 percent confidence interval: # 0.9258980 0.9835974 #sample estimates: # cor #0.9649507 ###Bill curvature summary(OGU2)#hel 23, phae 10 summary(OGU2$Bill curvature) phae<-subset(OGU2,Species=="Phaethornis") #Only Ohaethornis from OGU1 heli<-subset(OGU2,Species=="Heliconia") #Only heliconia from OGU1 heli_curva<-sample(heli$Bill curvature,10,replace=F) heli_curva #Pearson Parametric Correlation using ordered values cor.test(sort(phae$Bill curvature),sort(heli_curva),method="pearson") plot(sort(phae$Bill curvature),sort(heli_curva),main="Correlation Bill curvature OGU 2",xlab="Bill curvature Phaethornis",ylab="Bill curvature Heliconia",pch=20) #Result is: #data: sort(phae$Bill curvature) and sort(heli_curva) #t = 8.675, df = 8, p-value = 2.427e-05 #alternative hypothesis: true correlation is not equal to 0 #95 percent confidence interval: # 0.8000125 0.9885876 #sample estimates: # cor #0.9507417 ###Bill width summary(OGU2)#hel 23, phae 12 summary(OGU2$Bill width) phae<-subset(OGU2,Species=="Phaethornis") #Only Ohaethornis from OGU1 heli<-subset(OGU2,Species=="Heliconia") #Only heliconia from OGU1 heli_an<-sample(heli$Bill width,12,replace=F) heli_an #Pearson Parametric Correlation using ordered values cor.test(sort(phae$Bill width),sort(heli_an),method="pearson") plot(sort(phae$Bill width),sort(heli_an),main="Correlation Bill width OGU 2",xlab="Bill width Phaethornis",ylab="Bill width Heliconia", pch=20) #Result is: #data: sort(phae$Bill width) and sort(heli_an) #t = 5.0674, df = 10, p-value = 0.0004867 #alternative hypothesis: true correlation is not equal to 0 #95 percent confidence interval: # 0.5348793 0.9565444 #sample estimates: # cor #0.8483611 ####Finally, working on OGU 3 (Oaxaca)####### OGU3<-subset(base,OGU==3) summary(OGU3) ####First, working on Bill length boxplot(OGU3$Bill length~OGU3$Species, main="Length OGU 3",ylab="mm",cex.labs=2) #Testing for normality hist(OGU3$Bill length,col="red") #Expecting normality??? shapiro.test(OGU3$Bill length) #Result is: #W = 0.76752, p-value = 5.915e-11 #Not normally distributed ##Testing for significant differences wi.Bill length.OGU3<-wilcox.test(Bill length~Species,data=OGU3) wi.Bill length.OGU3 #Result is: #Wilcoxon rank sum test with continuity correction #data: Bill length by Species #W = 1533, p-value = 1.225e-10 #alternative hypothesis: true location shift is not equal to 0 #There are significant differences!!!!!! ####Now, working on Bill curvature boxplot(OGU3$Bill curvature~OGU3$Species, main="Curvature OGU 3", ylab="mm",cex.labs=2) #Testing for normality hist(OGU3$Bill curvature,col="red") #Expecting normality?? shapiro.test(OGU3$Bill curvature) #Result is: #W = 0.88863, p-value = 3.694e-06 #Not normally distributed ##Testing for significant differences wi.curva.OGU3<-wilcox.test(Bill curvature~Species,data=OGU3) wi.curva.OGU3 #Result is: #Wilcoxon rank sum test with continuity correction #data: Bill curvature by Species #W = 173, p-value = 7.817e-06 #alternative hypothesis: true location shift is not equal to 0 #There are significant differences!!!!!! ###Now, working on Bill width boxplot(OGU3$Bill width~OGU3$Species,main="Width OGU 3", ylab="mm",cex.labs=2) #Testing for normality hist(OGU3$Bill width,col="red") #Expecting normality?? shapiro.test(OGU3$Bill width) #el resultado es: #W = 0.82653, p-value = 7.856e-09 #Not normally distributed ##Testing for significant differences wi.an.OGU3<-wilcox.test(Bill width~Species,data=OGU3) wi.an.OGU3 #Result is: # #data: Bill width by Species #data: Bill width by Species #W = 1474, p-value = 2.44e-12 #There are significant differences!!!!!! ####Building correlations ###Bill length summary(OGU3)#hel 22, phae 75 phae<-subset(OGU3,Species=="Phaethornis") #Only Ohaethornis from OGU3 heli<-subset(OGU3,Species=="Heliconia") #Only heliconia from OGU3 phae_Bill length<-sample(phae$Bill length,22,replace=F) #Obtenemos 22 valores de Bill length de heliconias phae_Bill length #Correlacion paremtrica de earson usando los valores ordenados cor.test(sort(heli$Bill length),sort(phae_Bill length),method="pearson") plot(sort(heli$Bill length),sort(phae_Bill length),main="Correlación Bill length OGU 3",xlab="Bill length Heliconia", ylab="Bill length Phaethornis",pch=20) #El resultado es: #data: sort(heli$Bill length) and sort(phae_Bill length) #t = 6.4541, df = 20, p-value = 2.708e-06 #alternative hypothesis: true correlation is not equal to 0 #95 percent confidence interval: # 0.6126667 0.9235246 #sample estimates: # cor #0.8219587 ###Bill curvature summary(OGU3)#hel 18, phae 75 summary(OGU3$Bill curvature) phae<-subset(OGU3,Species=="Phaethornis") #Only Ohaethornis from OGU3 heli<-subset(OGU3,Species=="Heliconia") #Only heliconia from OGU3 phae_curva<-sample(phae$Bill curvature,18,replace=F) #Obtenemos 18 valores de Bill length de heliconias phae_curva #Pearson Parametric Correlation using ordered values cor.test(sort(heli$Bill curvature),sort(phae_curva),method="pearson") plot(sort(heli$Bill curvature),sort(phae_curva),main="Correlación Bill curvature OGU 3", xlab="Bill curvature Heliconia",ylab="Bill curvature Phaethornis",pch=20) #Result is: #data: sort(heli$Bill curvature) and sort(phae_curva) #t = 7.1849, df = 16, p-value = 2.171e-06 #alternative hypothesis: true correlation is not equal to 0 #95 percent confidence interval: # 0.6871523 0.9521830 #sample estimates: # cor #0.8737227 ###Bill width summary(OGU3)#hel 22, phae 75 summary(OGU3$Bill width) phae<-subset(OGU3,Species=="Phaethornis") #Only Ohaethornis from OGU3 heli<-subset(OGU3,Species=="Heliconia") #Only heliconia from OGU3 phae_an<-sample(phae$Bill width,22,replace=F) #Obtenemos 22 valores de Bill length de heliconias phae_an #Pearson Parametric Correlation using ordered values cor.test(sort(heli$Bill width),sort(phae_an),method="pearson") plot(sort(heli$Bill width),sort(phae_an),main="Coorelación Bill width OGU 3",xlab="Bill width Heliconia",ylab="Bill width Phaethornis",pch=20) #Result is: # data: sort(heli$Bill width) and sort(phae_an) #t = 17.442, df = 20, p-value = 1.447e-13 #alternative hypothesis: true correlation is not equal to 0 #95 percent confidence interval: # 0.9247051 0.9871320 #sample estimates: # cor #0.968666 ######Building Principal Components Analysis (PCA) for each OGU!! #Using OGU1 (Nayarit), OGU2 (Jalisco) and OGU3 (Oaxaca) #Intalling package with the function for PCA analyses install.packages("car") #Load package library(car) #First OGU1 PCA<-princomp(~Bill length+Bill curvature+Bill width, data=OGU1, scores=T, cor=T) Species.PCA<-c() for(j in row.names(PCA$scores)){ Species.PCA<-c(Species.PCA,OGU1[j,"Species"],use.names=T) } #Factor loadings: PCA$loadings[,] #Plot (without confidence elipses) plot(PCA$scores[,1],PCA$scores[,2],col=Species.PCA, pch=19,main="OGU 1",ylab="PCA component 2",xlab="PCA component 1") #text(PCA$scores[,1],PCA$scores[,2],Species.PCA, col=Species.PCA) abline(h=0, lty=2) abline(v=0, lty=2) #Building confidence elipses on the PCA plot dataEllipse(PCA$scores[,1],PCA$scores[,2], #Coordinates "x" and "y" as.factor(Species.PCA), #Grouping factor col=c(1,2), # Colours, one per factor levels=c(0.95), #Percentage of points to be included within the confidence ellipse fill=T, #Filling out the ellipse? group.labels=c("Heliconia","Phaethornis"), #Tags for the ellipses xlab="CP1",ylab="CP2",main="", #Plot titles grid=F) #Drawing a grid? abline(h=0,lty=2) #x=0 and y=0 abline(v=0,lty=2) library(scatterplot3d) scatterplot3d(OGU1$Bill length,OGU1$Bill curvature,OGU1$Bill width, color=c(OGU1$Species)) #Now working on OGU2 PCA<-princomp(~Bill length+Bill curvature+Bill width, data=OGU2, scores=T, cor=T) Species.PCA<-c() for(j in row.names(PCA$scores)){ Species.PCA<-c(Species.PCA,OGU2[j,"Species"],use.names=T) } #Factor loadings: PCA$loadings[,] plot(PCA$scores[,1],PCA$scores[,2],col=Species.PCA, pch=19,main="OGU 2",ylab="PCA component 2",xlab="Component 1") #text(PCA$scores[,1],PCA$scores[,2],Species.PCA, col=Species.PCA) abline(h=0, lty=2) abline(v=0, lty=2) #Building confidence elipses on the PCA plot dataEllipse(PCA$scores[,1],PCA$scores[,2], #Coordinates "x" and "y" as.factor(Species.PCA), #Grouping factor col=c(1,2), # Colours, one per factor levels=c(0.95), #Percentage of points to be included within the confidence ellipse fill=T, #Filling out the ellipse? group.labels=c("Heliconia","Phaethornis"), #Tags for the ellipses xlab="CP1",ylab="CP2",main="", #Plot titles grid=F) #Drawing a grid? abline(h=0,lty=2) ##x=0 and y=0 abline(v=0,lty=2) scatterplot3d(OGU2$Bill length,OGU2$Bill curvature,OGU2$Bill width, color=c(OGU2$Species)) #Finally, working on OGU3 PCA<-princomp(~Bill length+Bill curvature+Bill width, data=OGU3, scores=T, cor=T) Species.PCA<-c() for(j in row.names(PCA$scores)){ Species.PCA<-c(Species.PCA,OGU3[j,"Species"],use.names=T) } #Factor loadings: PCA$loadings[,] plot(PCA$scores[,1],PCA$scores[,2],col=Species.PCA, pch=19,main="OGU 3",ylab="PCA Componente 2",xlab="Componente 1") #text(PCA$scores[,1],PCA$scores[,2],Species.PCA, col=Species.PCA) abline(h=0, lty=2) abline(v=0, lty=2) dataEllipse(PCA$scores[,1],PCA$scores[,2], #Coordinates "x" and "y" as.factor(Species.PCA), #Grouping factor col=c(1,2), # Colours, one per factor levels=c(0.95), #Percentage of points to be included within the confidence ellipse fill=T, #Filling out the ellipse? group.labels=c("Heliconia","Phaethornis"), #Tags for the ellipses xlab="CP1",ylab="CP2",main="", #Plot titles grid=F) #Drawing a grid? abline(h=0,lty=2) #x=0 and y=0 abline(v=0,lty=2) scatterplot3d(OGU3$Bill length,OGU3$Bill curvature,OGU3$Bill width, color=c(OGU3$Species)) ################Plot for all variables and tests for each OGU############################# base<-read.csv(file="PhaetornisHeliconia2.csv", header=TRUE) names(base) summary(base) boxplot(base$Bill length~base$OGU, col="blue") boxplot(base$Bill curvature~base$OGU, col="blue") boxplot(base$Bill width~base$OGU, col="blue") #### First, the Phaethornis hummingbirds phaethornis<-subset(base,Species=="Phaethornis") summary(phaethornis) dev.off() boxplot(phaethornis$Bill length~phaethornis$OGU, xlab="OGU", ylab="mm", main="bill length",cex.labs=2) boxplot(phaethornis$Bill curvature~phaethornis$OGU, xlab="OGU", ylab="mm", main="bill curvature",cex.labs=2) boxplot(phaethornis$Bill width~phaethornis$OGU, xlab="OGU", ylab="mm", main="bill width",cex.labs=2) ####Testing normality and between OGUs differeces### #Bill length phaBill length<-shapiro.test(phaethornis$Bill length) phaBill length #Result is #data: phaethornis$Bill length #W = 0.98984, p-value = 0.5521 #They are normally distributed #Using ANOVA (aov) phaBill lengthanova<-aov(Bill length~OGU,data=phaethornis) summary(phaBill lengthanova) #Resultado is: # Df Sum Sq Mean Sq F value Pr(>F) #OGU 1 107.9 107.90 24.59 2.52e-06 *** # Residuals 113 495.9 4.39 --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #2 observations deleted due to missingness #Bill width phaBill width<-shapiro.test(phaethornis$Bill width) phaBill width #Result is: #data: phaethornis$Bill width #W = 0.95781, p-value = 0.001838 #Not normally distributed #We then used a Kruskal-Wallis test phaBill widthkruskal<-kruskal.test(Bill width~OGU,data=phaethornis) phaBill widthkruskal #Result is: # Kruskal-Wallis chi-squared = 2.191, df = 2, p-value = 0.3344 #Bill curvature phacurva<-shapiro.test(phaethornis$Bill curvature) phacurva #Result is: #data: phaethornis$Bill width #data: phaethornis$Bill curvature #W = 0.74446, p-value = 1.246e-11 #Not normally distributed #We then used a Kruskal-Wallis test phacurvakruskal<-kruskal.test(Bill curvature~OGU,data=phaethornis) phacurvakruskal #Result is: # Kruskal-Wallis chi-squared = 0.9019, df = 2, p-value = 0.637 #### Then, the Heliconia flowers heliconias<-subset(base,Species=="Heliconia") summary(heliconias) dev.off() boxplot(heliconias$Bill length~heliconias$OGU, xlab="OGU", ylab="mm", main="corolla length",cex.labs=2) boxplot(heliconias$Bill curvature~heliconias$OGU, xlab="OGU", ylab="mm", main="corolla curvature",cex.labs=2) boxplot(heliconias$Bill width~heliconias$OGU, xlab="OGU", ylab="mm", main="corolla width",cex.labs=2) ####Testing normality and between OGUs differeces### #Bill length helBill length<-shapiro.test(heliconias$Bill length) helBill length #Result is: #data: heliconias$Bill length #W = 0.87962, p-value = 1.958e-07 #Not normally distributed #Using ANOVA (Kruskal) helBill lengthanova<-kruskal.test(Bill length~OGU,data=heliconias) helBill lengthanova #Result is: # Kruskal-Wallis chi-squared = 58.934, df = 2, p-value = 1.595e-13 #Bill width helBill width<-shapiro.test(heliconias$Bill width) helBill width #Result is: #data: heliconias$Bill width #W = 0.98033, p-value = 0.2142 #They are normally distributed #Using ANOVA helBill widthanova<-aov(Bill width~OGU,data=heliconias) summary(helBill widthanova) #Result is: # Df Sum Sq Mean Sq F value Pr(>F) #OGU 1 23.19 23.186 26.14 1.97e-06 *** # Residuals 84 74.50 0.887 --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #13 observations deleted due to missingness #Bill curvature helcurva<-shapiro.test(heliconias$Bill curvature) helcurva #Result is: #data: heliconias$Bill curvature #W = 0.97297, p-value = 0.05499 #They are normally distributed #Using ANOVA helcurvanova<-aov(Bill curvature~OGU,data=heliconias) summary(helcurvanova) #Result is: # Df Sum Sq Mean Sq F value Pr(>F) #OGU 1 0.003790 0.00379 38.01 2.02e-08 *** # Residuals 89 0.008877 0.00010 --- # Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #8 observations deleted due to missingness ### Canonic Correspondence Analysis (CCA) ### base<-read.csv(file="PhaetornisHeliconia2.csv", header=TRUE) names(base) summary(base) #install.packages ("CCA") library(CCA) ####Working first on OGU 1(Nayarit) OGU1<-subset(base,OGU==1) summary(OGU1) phae<-subset(OGU1,Species=="Phaethornis")[,14:16] heli<-subset(OGU1,Species=="Heliconia")[,14:16] heli<-heli[sample(1:54,29,replace=F),] #CCA cc_OGU1<-cc(phae,heli) cc_OGU1 barplot(cc_OGU1$cor, xlab = "Dimension",ylab = "Canonical correlations", names.arg =1:3, ylim = c(0,1)) plt.cc(cc_OGU1) ####Working on OGU 2 (Jalisco) OGU2<-subset(base,OGU==2) summary(OGU2) phae<-subset(OGU2,Species=="Phaethornis")[,14:16] heli<-subset(OGU2,Species=="Heliconia")[,14:16] heli<-heli[sample(1:23,13,replace=F),] #CCA cc_OGU2<-cc(phae,heli) cc_OGU2 barplot(cc_OGU2$cor, xlab = "Dimension",ylab = "Canonical correlations", names.arg =1:3, ylim = c(0,1)) plt.cc(cc_OGU2) ####Working on OGU 3 (Oaxaca) OGU3<-subset(base,OGU==3) summary(OGU3) phae<-subset(OGU3,Species=="Phaethornis")[,14:16] heli<-subset(OGU3,Species=="Heliconia")[,14:16] phae<-phae[sample(1:75,22,replace=F),] #CCA cc_OGU3<-cc(phae,heli) cc_OGU3 barplot(cc_OGU3$cor, xlab = "Dimension",ylab = "Canonical correlations", names.arg =1:3, ylim = c(0,1)) plt.cc(cc_OGU3) ####Variables correlations### base<- phae<-subset(base,Species=="Phaethornis") cor(na.omit(phae[,14:16])) heli<-subset(base,Species=="Heliconia") cor(na.omit(heli[,14:16]))