Thursday, June 6, 2013

There are 145 unique species of babblers to be used in my analysis.

> library(plyr)
> babblers.species.count<-ddply(.data=new.babblers, .(Genus, Species), species.count=length(Wing.Length), summarize)
> babblers.species.count

                   Genus         Species species.count
1             Actinodura        egertoni             6
2             Actinodura    morrisoniana             3
3             Actinodura      nipalensis             3
4             Actinodura         waldeni             5
5                Alcippe    brunneicauda             4
6                Alcippe          grotei            11
7                Alcippe      morrisonia             5
8                Alcippe      peracensis             2
9                Alcippe    poioicephala            11
10               Alcippe     rufogularis             5
11                 Babax      lanceolata             4
12             Chrysomma         sinense             6
13  Chrysomma (Moupinia)      poecilotis             4
14                 Cutia      nipalensis             7
15              Dumentia      hyperythra             6
16        Gampsorhynchus         rufulus             8
17              Garrulax         affinis            11
18              Garrulax     albogularis             5
19              Garrulax         canorus            14
20              Garrulax       chinensis            40
21              Garrulax      cineraceus             4
22              Garrulax erythrocephalus             1
23              Garrulax        formosus             4
24              Garrulax     leucolophus            41
25              Garrulax        lineatus             6
26              Garrulax           maesi             3
27              Garrulax       merulinus             1
28              Garrulax         milleti             1
29              Garrulax          milnei             4
30              Garrulax        mitratus             6
31              Garrulax       monileger             8
32              Garrulax        nuchalis             3
33              Garrulax       ocellatus             7
34              Garrulax       palliatus             4
35              Garrulax      pectoralis             1
36              Garrulax  perspicillatus            14
37              Garrulax poecilorhynchus             4
38              Garrulax          sannio            12
39              Garrulax         sinesis             9
40              Garrulax       squamatus             3
41              Garrulax        striatus             6
42              Garrulax     subunicolor             6
43              Garrulax      variegatum             3
44              Garrulax         vassali             7
45              Garrulax        virgatus             2
46         Heterophaesia       annectens             4
47         Heterophaesia     auricularis             3
48         Heterophaesia      capistrata             6
49         Heterophaesia      desgodinsi             1
50         Heterophaesia        gracilis             4
51         Heterophaesia    maelanoleuca             9
52         Heterophaesia       picaoides             7
53         Heterophaesia       pulchella             3
54          Heterophasia      desgodinsi             7
55           Illadopsis       albipectus             4
56           Illadopsis         cleaveri             4
57           Illadopsis       fulvescens             8
58           Illadopsis           puveli             2
59           Illadopsis      pyrrhoptera             2
60           Illadopsis        rufescens             3
61           Illadopsis       rufipennis             6
62               Kenopia         striata             2
63             Leiothrix     argentauris            15
64             Leiothrix           lutea             4
65            Leonardina           woodi             1
66             Liocichla       phoenicea             1
67             Liocichla         steerii             3
68              Macronus         gularis            25
69              Macronus        ptilosus             3
70              Macronus     striaticeps            12
71          Malacocincla         abbotti             4
72          Malacocincla     cinereiceps             2
73          Malacocincla     malaccensis             6
74          Malacocincla        sepiaria             3
75          Malacopteron          affine             3
76          Malacopteron      albogulare             4
77          Malacopteron        cinereum             1
78          Malacopteron     magnirostre             7
79          Malacopteron          magnum             5
80          Malacopteron    palawanenese             3
81                 Malia           grata             6
82         Micromacronus        sordidus             2
83                 Minla      ignotincta             6
84                 Minla        strigula             7
85             Napothera    brevicaudata             4
86             Napothera          crassa             2
87             Napothera     crispifrons             7
88             Napothera     epilepidota            10
89            Pellorneum      albiventre             9
90            Pellorneum    capristratum            15
91            Pellorneum      pyrrogenys             6
92            Pellorneum        ruficeps             5
93            Pellorneum        tickelli             1
94           Phyllanthus      atripennis             6
95              Pnoepyga         pusilla             1
96          Pomatorhinus    erythrogenys            18
97          Pomatorhinus    ferruginosus             4
98          Pomatorhinus      horsfieldi             6
99          Pomatorhinus      hypoleucos             4
100         Pomatorhinus        montanus             4
101         Pomatorhinus      ruficollis             9
102         Pomatorhinus        swinhoei             2
103          Ptilocichla         falcata             2
104          Ptilocichla   leucogrammica             3
105          Ptilocichla     mindanensis             3
106            Ptyrticus        turdinus             6
107         Rhopocichla         atriceps             4
108              Rimator    malacoptilus             3
109         Schoeniparus         brunnea            10
110         Schoeniparus     castaneceps             1
111         Schoeniparus         cinerea             3
112         Schoeniparus           dubia             6
113         Schoeniparus     rufogularis             4
114          Spelaeornis    chocolatinus             5
115          Spelaeornis        formosus             2
116          Spelaeornis   troglodytoide             3
117         Sphenocichla           humei             3
118            Stachyris         ambigua             4
119            Stachyris        chrysaea            11
120            Stachyris    erythroptera             8
121            Stachyris      grammiceps             2
122            Stachyris        leucotis             3
123            Stachyris        maculata             3
124            Stachyris    melanothorax             6
125            Stachyris       nigriceps             7
126            Stachyris     nigricollis             4
127            Stachyris           oglei             3
128            Stachyris    poliocephala             4
129            Stachyris        pyrrhops             3
130            Stachyris        ruficeps            10
131            Stachyris       rufifrons             3
132            Stachyris       striolata             5
133            Stachyris       thoracica             4
134              Timalia         pileata             8
135          Trichastoma         bicolor             4
136          Trichastoma       celebense            11
137          Trichastoma       rostratum             5
138          Trichastoma        tickelli            10
139          Trichastona        tickelli             2
140            Turdoides         bicolor             3
141            Turdoides         gularis             4
142            Turdoides       jardineii            10
143            Turdoides        plebejus            16
144            Turdoides      reinwardii             3
145        Xiphirhynchus   superciliaris             3

Had a minor panic attack because I realized my missing tarsus values were being imputed with mean of the bill length which of course would lead to extreme variance in the trait. Anyways, fixed it right up and still got somewhat similar results. Whew.
New PCA result:


Fixed script:
#impute missing values by taking mean of all other values of of other species ofthat character 

i<-1 #set the index
species.removed<-c()
#sift through all babbler species
while (i <= length(babblers$UniqueID_1)){
  #missing wing lengths; if the Wing Length is NA for that value
  if(is.na(babblers$Wing.Length[i]) == TRUE){
    #select the wing lengths of all birds which match this species with missing wing length
    wing.species.vector<-babblers$Wing.Length[which(babblers$Species[i] == babblers$Species)]
    #make a new vector without the missing value
    wing.species.filler<-wing.species.vector[which(is.na(wing.species.vector) == FALSE)]
    #calculate the mean
    wing.impute<-mean(wing.species.filler)
    #fill in missing value with the mean for that bird
    babblers$Wing.Length[i]<-wing.impute
  }
  
  #repeat the above method for all characters
  if(is.na(babblers$Tail.Length[i]) == TRUE){
    tail.species.vector<-babblers$Tail.Length[which(babblers$Species[i] == babblers$Species)]
    tail.species.filler<-tail.species.vector[which(is.na(tail.species.vector) == FALSE)]
    tail.impute<-mean(tail.species.filler)
    babblers$Tail.Length[i]<-tail.impute
  }
  
  
  
  if(is.na(babblers$Bill.Length[i]) == TRUE){
    BL.species.vector<-babblers$Bill.Length[which(babblers$Species[i] == babblers$Species)]
    BL.species.filler<-BL.species.vector[which(is.na(BL.species.vector) == FALSE)]
    BL.impute<-mean(BL.species.filler)
    babblers$Bill.Length[i]<-BL.impute
  }
  
  if(is.na(babblers$Bill.Width[i]) == TRUE){
    BW.species.vector<-babblers$Bill.Width[which(babblers$Species[i] == babblers$Species)]
    BW.species.filler<-BW.species.vector[which(is.na(BW.species.vector) == FALSE)]
    BW.impute<-mean(BW.species.filler)
    babblers$Bill.Width[i]<-BW.impute
  }
  
  if(is.na(babblers$Bill.Depth[i]) == TRUE){
    BD.species.vector<-babblers$Bill.Depth[which(babblers$Species[i] == babblers$Species)]
    BD.species.filler<-BD.species.vector[which(is.na(BD.species.vector) == FALSE)]
    BD.impute<-mean(BD.species.filler)
    babblers$Bill.Depth[i]<-BD.impute
  }
  
  if(is.na(babblers$Tarsus.Length[i]) == TRUE){
    TL.species.vector<-babblers$Tarsus.Length[which(babblers$Species[i] == babblers$Species)]
    TL.species.filler<-TL.species.vector[which(is.na(TL.species.vector) == FALSE)]
    TL.impute<-mean(TL.species.filler)
    babblers$Tarsus.Length[i]<-TL.impute
  }
  
  if(is.na(babblers$Hallux[i]) == TRUE){
    hallux.species.vector<-babblers$Hallux[which(babblers$Species[i] == babblers$Species)]
    hallux.species.filler<-hallux.species.vector[which(is.na(hallux.species.vector) == FALSE)]
    hallux.impute<-mean(hallux.species.filler)
    babblers$Hallux[i]<-hallux.impute
  }

  i<-i+1
}


PCAs and reweights of the values for each species is finally done and in a format that can be easily used to run the evol.rate.mcmc()


Script to impute data:
#impute missing values by taking mean of all other values of of other species ofthat character 

i<-1 #set the index
species.removed<-c()
#sift through all babbler species
while (i <= length(babblers$UniqueID_1)){
  #missing wing lengths; if the Wing Length is NA for that value
  if(is.na(babblers$Wing.Length[i]) == TRUE){
    #select the wing lengths of all birds which match this species with missing wing length
    wing.species.vector<-babblers$Wing.Length[which(babblers$Species[i] == babblers$Species)]
    #make a new vector without the missing value
    wing.species.filler<-wing.species.vector[which(is.na(wing.species.vector) == FALSE)]
    #calculate the mean
    wing.impute<-mean(wing.species.filler)
    #fill in missing value with the mean for that bird
    babblers$Wing.Length[i]<-wing.impute
  }
  
  #repeat the above method for all characters
  if(is.na(babblers$Tail.Length[i]) == TRUE){
    tail.species.vector<-babblers$Tail.Length[which(babblers$Species[i] == babblers$Species)]
    tail.species.filler<-tail.species.vector[which(is.na(tail.species.vector) == FALSE)]
    tail.impute<-mean(tail.species.filler)
    babblers$Tail.Length[i]<-tail.impute
  }
  
  
  
  if(is.na(babblers$Bill.Length[i]) == TRUE){
    BL.species.vector<-babblers$Bill.Length[which(babblers$Species[i] == babblers$Species)]
    BL.species.filler<-BL.species.vector[which(is.na(BL.species.vector) == FALSE)]
    BL.impute<-mean(BL.species.filler)
    babblers$Bill.Length[i]<-BL.impute
  }
  
  if(is.na(babblers$Bill.Width[i]) == TRUE){
    BW.species.vector<-babblers$Bill.Width[which(babblers$Species[i] == babblers$Species)]
    BW.species.filler<-BW.species.vector[which(is.na(BW.species.vector) == FALSE)]
    BW.impute<-mean(BW.species.filler)
    babblers$Bill.Width[i]<-BW.impute
  }
  
  if(is.na(babblers$Bill.Depth[i]) == TRUE){
    BD.species.vector<-babblers$Bill.Depth[which(babblers$Species[i] == babblers$Species)]
    BD.species.filler<-BD.species.vector[which(is.na(BD.species.vector) == FALSE)]
    BD.impute<-mean(BD.species.filler)
    babblers$Bill.Depth[i]<-BD.impute
  }
  
  if(is.na(babblers$Tarsus.Length[i]) == TRUE){
    TL.species.vector<-babblers$Tarsus.Length[which(babblers$Species[i] == babblers$Species)]
    TL.species.filler<-TL.species.vector[which(is.na(TL.species.vector) == FALSE)]
    TL.impute<-mean(BD.species.filler)
    babblers$Tarsus.Length[i]<-TL.impute
  }
  
  if(is.na(babblers$Hallux[i]) == TRUE){
    hallux.species.vector<-babblers$Hallux[which(babblers$Species[i] == babblers$Species)]
    hallux.species.filler<-hallux.species.vector[which(is.na(hallux.species.vector) == FALSE)]
    hallux.impute<-mean(hallux.species.filler)
    babblers$Hallux[i]<-hallux.impute
  }

  i<-i+1
}

Script to scale by geometric mean:
#make empty columns for geometric mean scale
babblers$gm.wing<-0
babblers$gm.tail<-0
babblers$gm.BL<-0
babblers$gm.BW<-0
babblers$gm.BD<-0
babblers$gm.tarsus<-0
babblers$gm.hallux<-0

#Size-adjust the raw data by dividing each individual by the geometric mean 

j<-1
could.not.scale<-c()
while(j<=length(babblers$UniqueID_1)){
  #grab all the data of that species and put in vector
  temp<-c(babblers$Wing.Length[j], babblers$Tail.Length[j], babblers$Bill.Length[j],
          babblers$Bill.Width[j], babblers$Bill.Depth[j], babblers$Tarsus.Length[j], 
          babblers$Hallux[j])
  temp<-na.omit(temp)
  #calculate the geometric mean of meausrements for that species
  geo.mean<-exp(mean(log(temp)))
  
  if(is.nan(geo.mean)==FALSE){
  #divide all values for that species by the geometric mean
  geo.wing<-(babblers$Wing.Length[j])/geo.mean
  babblers$gm.wing[j]<-geo.wing
  
  geo.tail<-babblers$Tail.Length[j]/geo.mean
  babblers$gm.tail[j]<-geo.tail
  
  geo.BL<-babblers$Bill.Length[j]/geo.mean
  babblers$gm.BL[j]<-geo.BL
  
  geo.BW<-babblers$Bill.Width[j]/geo.mean
  babblers$gm.BW[j]<-geo.BW
  
  geo.BD<-babblers$Bill.Depth[j]/geo.mean
  babblers$gm.BD[j]<-geo.BD
  
  geo.tarsus<-babblers$Tarsus.Length[j]/geo.mean
  babblers$gm.tarsus[j]<-geo.tarsus
  
  geo.hallux<-babblers$Hallux[j]/geo.mean
  babblers$gm.hallux[j]<-geo.hallux
 }
 else{
    append(could.not.scale, babblers$UniqueID_1[j])
  }
  #reset the temp
  temp<-c()
  
  #increment to next species
  j<-j+1
}

cat("The following species could not be scaled:", could.not.scale, "\n")

Script to do PCA with geometric mean values and put rescaled scores into new data frame:
#PCA for all measurements
pca.result <- prcomp(~babblers$gm.wing + babblers$gm.tail + babblers$gm.BL + babblers$gm.BW + babblers$gm.BD + 
       babblers$gm.tarsus + babblers$gm.hallux, data = babblers, retx = T)

#standard deviation of PCA
sd <- pca.result$sdev
#make vector of weights
loadings <- pca.result$rotation
#value of the rotated data (scores)
scores <- pca.result$x

#make a new data frame
new.babblers<-as.data.frame(babblers)

#make empty columns in the data frame
new.babblers$pca1<-0
new.babblers$pca2<-0
new.babblers$pca3<-0
new.babblers$pca4<-0
new.babblers$pca5<-0
new.babblers$pca6<-0
new.babblers$pca7<-0

#put PCA Scores into the data frame
new.babblers$pca1<-scores[,1]
new.babblers$pca2<-scores[,2]
new.babblers$pca3<-scores[,3]
new.babblers$pca4<-scores[,4]
new.babblers$pca5<-scores[,5]
new.babblers$pca6<-scores[,6]
new.babblers$pca7<-scores[,7]

new.babblers$pca# are the traits that will be used in the Bayesian analysis.
This site was really useful for understanding PCA:
http://strata.uga.edu/software/pdf/pcaTutorial.pdf
I am dropping Actinodura ramsayi.

                                     
LY_0134
Actinodura ramsayi

8.3       NA 10.44 5.34 5.86 25.57        NA
LY_0135 Actinodura ramsayi 7.9       NA 10.86 6.24 5.64          NA                 NA
I only have two specimens for Actinodura ramsayi but am missing the tail and hallux measurements for both so cannot impute the missing data. It really messes up the PCA and I will not be able to get trait data for this node if I do not fill these in. Should I take the mean of those samples from the entire Actinodura genus? Or drop the species?

Tuesday, June 4, 2013

Updated morphological data to include measurements from samples from data from Vietnam. Strangely, variance does not seem to be held just in one trait. 





When running PRGmatic v1.7, if you get an error:

"Can't exec "./cap3": No such file or directory at ./PRGmatic.pl line 116, <STDIN> line 8.
contigFile is Species_01_trimmed.fasta.cap.contigs
No such file or directory at ./PRGmatic.pl line 127, <STDIN> line 8."


Go to PRGmatic folder in home directory (or wherever it is), click >Dependent Software and you probably need to unzip "cap3.macosx.intel64.zip". Unzip this this folder and move all of the contents to the PRGmatic main directory.