PICs
Goals
We are going to test for an evolutionary correlation between two traits in our dataset. We will use a variety of methods.
Read in tree and data and make them match
require(ape)## Loading required package: aperequire(geiger)## Loading required package: geigerbigTree <- read.tree("mammal_bird_final.txt")
fullData <- read.table("LS_ti_mammal_censored.txt", header=T)
# try this
name.check(bigTree, fullData)## $tree_not_data
## [1] "Accipiter_cooperii" "Accipiter_nisus"
## [3] "Accipiter_striatus" "Acomys_spinosissimus"
## [5] "Aethomys_chrysophilus" "Akodon_azarae"
## [7] "Alectoris_graeca" "Anas_platyrhynchos"
## [9] "Anhinga_anhinga" "Aotus_trivirgatus"
## [11] "Apodemus_agrarius" "Aptenodytes_patagonicus"
## [13] "Apteryx_australis" "Aramus_guarauna"
## [15] "Arctocebus_calabarensis" "Arctocephalus_forsteri"
## [17] "Arctocephalus_pusillus" "Ardea_herodias"
## [19] "Artibeus_jamaicensis" "Artibeus_lituratus"
## [21] "Arvicola_amphibius" "Asellia_tridens"
## [23] "Athene_cunicularia" "Atilax_paludinosus"
## [25] "Aythya_collaris" "Babyrousa_babyrussa"
## [27] "Balaenoptera_acutorostrata" "Balaenoptera_borealis"
## [29] "Barbastella_barbastellus" "Bathyergus_suillus"
## [31] "Bombycilla_garrulus" "Botaurus_lentiginosus"
## [33] "Brachyphylla_cavernarum" "Bubo_virginianus"
## [35] "Bubulcus_ibis" "Budorcas_taxicolor"
## [37] "Buteo_buteo" "Buteo_jamaicensis"
## [39] "Buteo_lineatus" "Butorides_striata"
## [41] "Cacatua_galerita" "Calonectris_diomedea"
## [43] "Camelus_dromedarius" "Canis_aureus"
## [45] "Canis_lupus" "Carduelis_carduelis"
## [47] "Carduelis_chloris" "Carollia_perspicillata"
## [49] "Castor_canadensis" "Casuarius_casuarius"
## [51] "Cathartes_aura" "Ceratotherium_simum"
## [53] "Cercopithecus_cephus" "Charadrius_vociferus"
## [55] "Cheirogaleus_major" "Cheiromeles_torquatus"
## [57] "Chiroderma_villosum" "Chrotopterus_auritus"
## [59] "Coccothraustes_coccothraustes" "Colobus_polykomos"
## [61] "Columba_livia" "Corvus_frugilegus"
## [63] "Coturnix_japonica" "Cricetomys_gambianus"
## [65] "Cricetulus_barabensis" "Cryptomys_hottentotus"
## [67] "Ctenomys_talarum" "Cynopterus_brachyotis"
## [69] "Dendrohyrax_arboreus" "Dendroica_coronata"
## [71] "Dendroica_pinus" "Desmodillus_auricularis"
## [73] "Diceros_bicornis" "Dromaius_novaehollandiae"
## [75] "Dugong_dugon" "Egretta_caerulea"
## [77] "Egretta_rufescens" "Egretta_thula"
## [79] "Eidolon_helvum" "Elephantulus_brachyrhynchus"
## [81] "Elephantulus_edwardii" "Elephantulus_intufi"
## [83] "Elephantulus_myurus" "Elephantulus_rupestris"
## [85] "Elephas_maximus" "Emballonura_semicaudata"
## [87] "Enhydra_lutris" "Eonycteris_spelaea"
## [89] "Epomops_franqueti" "Eptesicus_serotinus"
## [91] "Equus_burchellii" "Equus_caballus"
## [93] "Erethizon_dorsatum" "Erithacus_rubecula"
## [95] "Eudyptes_chrysocome" "Eudyptes_chrysolophus"
## [97] "Eudyptula_minor" "Eumetopias_jubatus"
## [99] "Euoticus_elegantulus" "Falco_mexicanus"
## [101] "Falco_sparverius" "Falco_subbuteo"
## [103] "Felis_silvestris" "Ficedula_hypoleuca"
## [105] "Galago_alleni" "Garrulus_glandarius"
## [107] "Gavia_immer" "Gazella_dorcas"
## [109] "Genetta_servalina" "Gerbilliscus_brantsii"
## [111] "Gerbilliscus_leucogaster" "Glossophaga_soricina"
## [113] "Glyphorynchus_spirurus" "Graphiurus_murinus"
## [115] "Hemitragus_jemlahicus" "Hexaprotodon_liberiensis"
## [117] "Hippopotamus_amphibius" "Hipposideros_diadema"
## [119] "Hydrochoerus_hydrochaeris" "Hydrurga_leptonyx"
## [121] "Hylobates_pileatus" "Hypsignathus_monstrosus"
## [123] "Ixobrychus_exilis" "Lagenorhynchus_obliquidens"
## [125] "Lagopus_lagopus" "Lagopus_muta"
## [127] "Lagothrix_lagotricha" "Lanius_ludovicianus"
## [129] "Larus_argentatus" "Lasiopodomys_brandtii"
## [131] "Lasiurus_borealis" "Laterallus_jamaicensis"
## [133] "Lavia_frons" "Lemniscomys_griselda"
## [135] "Lepilemur_leucopus" "Leptonychotes_weddellii"
## [137] "Lobodon_carcinophaga" "Loxia_pytyopsittacus"
## [139] "Loxodonta_africana" "Macaca_fascicularis"
## [141] "Macropus_eugenii" "Macroscelides_proboscideus"
## [143] "Marmota_marmota" "Megaceryle_alcyon"
## [145] "Megascops_asio" "Melanerpes_carolinus"
## [147] "Melogale_moschata" "Melopsittacus_undulatus"
## [149] "Melospiza_melodia" "Meriones_unguiculatus"
## [151] "Mesoplodon_bowdoini" "Mesoplodon_hectori"
## [153] "Mesoplodon_layardii" "Micaelamys_namaquensis"
## [155] "Microtus_agrestis" "Microtus_longicaudus"
## [157] "Microtus_ochrogaster" "Microtus_pennsylvanicus"
## [159] "Microtus_pinetorum" "Mimus_polyglottos"
## [161] "Miniopterus_schreibersii" "Monachus_schauinslandi"
## [163] "Morus_bassanus" "Morus_capensis"
## [165] "Motacilla_alba" "Mus_minutoides"
## [167] "Mustela_nivalis" "Mycteria_americana"
## [169] "Myodes_gapperi" "Myodes_glareolus"
## [171] "Myotis_daubentonii" "Myotis_myotis"
## [173] "Myotis_mystacinus" "Myotis_nattereri"
## [175] "Myotomys_sloggetti" "Myotomys_unisulcatus"
## [177] "Mystromys_albicaudatus" "Nandinia_binotata"
## [179] "Nasalis_larvatus" "Neophoca_cinerea"
## [181] "Neophocaena_phocaenoides" "Noctilio_leporinus"
## [183] "Notomys_alexis" "Nyctalus_noctula"
## [185] "Nyctanassa_violacea" "Nyctereutes_procyonoides"
## [187] "Nycteris_grandis" "Nycteris_thebaica"
## [189] "Nycticorax_nycticorax" "Nymphicus_hollandicus"
## [191] "Oceanites_oceanicus" "Okapia_johnstoni"
## [193] "Ommatophoca_rossii" "Opisthocomus_hoazin"
## [195] "Orcinus_orca" "Ornithorhynchus_anatinus"
## [197] "Oryctolagus_cuniculus" "Otomys_angoniensis"
## [199] "Otomys_irroratus" "Pagophilus_groenlandicus"
## [201] "Pan_troglodytes" "Pandion_haliaetus"
## [203] "Panthera_leo" "Panthera_onca"
## [205] "Panthera_tigris" "Parotomys_brantsii"
## [207] "Parotomys_littledalei" "Parus_caeruleus"
## [209] "Parus_major" "Passer_domesticus"
## [211] "Passer_montanus" "Pecari_tajacu"
## [213] "Pelecanus_occidentalis" "Perdix_perdix"
## [215] "Pernis_apivorus" "Peromyscus_leucopus"
## [217] "Peromyscus_maniculatus" "Petrodromus_tetradactylus"
## [219] "Phalacrocorax_auritus" "Phenacomys_intermedius"
## [221] "Phoca_vitulina" "Phocoena_phocoena"
## [223] "Phodopus_roborovskii" "Phoebastria_immutabilis"
## [225] "Phoebetria_fusca" "Phoenicurus_ochruros"
## [227] "Phoenicurus_phoenicurus" "Phylloscopus_collybita"
## [229] "Phylloscopus_trochilus" "Phyllostomus_hastatus"
## [231] "Physeter_catodon" "Phytotoma_rara"
## [233] "Pinicola_enucleator" "Pipistrellus_kuhlii"
## [235] "Platyrrhinus_vittatus" "Plecotus_auritus"
## [237] "Pluvialis_squatarola" "Poiana_richardsonii"
## [239] "Pontoporia_blainvillei" "Porphyrio_hochstetteri"
## [241] "Potamogale_velox" "Prunella_modularis"
## [243] "Pseudomys_hermannsburgensis" "Pseudorca_crassidens"
## [245] "Pteronotus_parnellii" "Pteropus_alecto"
## [247] "Pteropus_giganteus" "Pteropus_poliocephalus"
## [249] "Pteropus_vampyrus" "Puffinus_gravis"
## [251] "Puffinus_lherminieri" "Pusa_hispida"
## [253] "Pygoscelis_adeliae" "Pygoscelis_papua"
## [255] "Pyrrhula_pyrrhula" "Rangifer_tarandus"
## [257] "Rattus_norvegicus" "Rattus_rattus"
## [259] "Regulus_regulus" "Rhabdomys_pumilio"
## [261] "Rhea_americana" "Rhinoceros_unicornis"
## [263] "Rhinolophus_clivosus" "Rhinolophus_ferrumequinum"
## [265] "Rhinolophus_hipposideros" "Rhynchonycteris_naso"
## [267] "Rousettus_aegyptiacus" "Saccostomus_campestris"
## [269] "Sciurus_carolinensis" "Sciurus_niger"
## [271] "Scotophilus_kuhlii" "Sitta_carolinensis"
## [273] "Sotalia_fluviatilis" "Sousa_chinensis"
## [275] "Spermophilus_dauricus" "Spheniscus_demersus"
## [277] "Spheniscus_magellanicus" "Steatomys_pratensis"
## [279] "Stenella_coeruleoalba" "Struthio_camelus"
## [281] "Sturnira_lilium" "Sturnus_vulgaris"
## [283] "Sus_scrofa" "Syconycteris_australis"
## [285] "Sylvia_atricapilla" "Sylvia_borin"
## [287] "Symphalangus_syndactylus" "Tachyglossus_aculeatus"
## [289] "Tadarida_aegyptiaca" "Taphozous_melanopogon"
## [291] "Tapirus_indicus" "Tayassu_pecari"
## [293] "Tetrao_urogallus" "Thallomys_paedulcus"
## [295] "Thomomys_bottae" "Trichechus_manatus"
## [297] "Trichoglossus_haematodus" "Troglodytes_aedon"
## [299] "Troglodytes_troglodytes" "Tscherskia_triton"
## [301] "Turdus_grayi" "Turdus_merula"
## [303] "Turdus_migratorius" "Turdus_philomelos"
## [305] "Tursiops_truncatus" "Tyrannus_melancholicus"
## [307] "Tyto_alba" "Uromys_caudimaculatus"
## [309] "Ursus_arctos" "Vampyressa_pusilla"
## [311] "Vespertilio_murinus" "Vulpes_lagopus"
## [313] "Vulpes_vulpes" "Zalophus_californianus"
## [315] "Zenaida_macroura"
##
## $data_not_tree
## [1] "1" "10" "100" "101" "102" "103" "104" "105" "106" "107" "108"
## [12] "109" "11" "110" "111" "112" "113" "114" "115" "116" "117" "118"
## [23] "119" "12" "120" "121" "122" "123" "124" "125" "126" "127" "128"
## [34] "129" "13" "130" "131" "132" "133" "134" "135" "136" "137" "138"
## [45] "139" "14" "140" "141" "142" "143" "144" "145" "146" "147" "148"
## [56] "149" "15" "150" "151" "152" "153" "154" "155" "156" "157" "158"
## [67] "159" "16" "160" "161" "162" "163" "164" "165" "166" "167" "168"
## [78] "169" "17" "170" "171" "172" "173" "174" "175" "176" "177" "178"
## [89] "179" "18" "180" "181" "182" "183" "184" "185" "186" "187" "188"
## [100] "189" "19" "190" "191" "192" "193" "194" "195" "196" "197" "198"
## [111] "2" "20" "21" "22" "23" "24" "25" "26" "27" "28" "29"
## [122] "3" "30" "31" "32" "33" "34" "35" "36" "37" "38" "39"
## [133] "4" "40" "41" "42" "43" "44" "45" "46" "47" "48" "49"
## [144] "5" "50" "51" "52" "53" "54" "55" "56" "57" "58" "59"
## [155] "6" "60" "61" "62" "63" "64" "65" "66" "67" "68" "69"
## [166] "7" "70" "71" "72" "73" "74" "75" "76" "77" "78" "79"
## [177] "8" "80" "81" "82" "83" "84" "85" "86" "87" "88" "89"
## [188] "9" "90" "91" "92" "93" "94" "95" "96" "97" "98" "99"# problem is species names in the data file
# geiger and other packages assume that names are in rownames
# make species names the rownames of the matrix
mammalData <- fullData[,-1]
rownames(mammalData) <- fullData[,1]
# is it better now?
name.check(bigTree, mammalData)## $tree_not_data
## [1] "Accipiter_cooperii" "Accipiter_nisus"
## [3] "Accipiter_striatus" "Alectoris_graeca"
## [5] "Anas_platyrhynchos" "Anhinga_anhinga"
## [7] "Aptenodytes_patagonicus" "Apteryx_australis"
## [9] "Aramus_guarauna" "Ardea_herodias"
## [11] "Athene_cunicularia" "Aythya_collaris"
## [13] "Bombycilla_garrulus" "Botaurus_lentiginosus"
## [15] "Bubo_virginianus" "Bubulcus_ibis"
## [17] "Buteo_buteo" "Buteo_jamaicensis"
## [19] "Buteo_lineatus" "Butorides_striata"
## [21] "Cacatua_galerita" "Calonectris_diomedea"
## [23] "Carduelis_carduelis" "Carduelis_chloris"
## [25] "Casuarius_casuarius" "Cathartes_aura"
## [27] "Charadrius_vociferus" "Coccothraustes_coccothraustes"
## [29] "Columba_livia" "Corvus_frugilegus"
## [31] "Coturnix_japonica" "Dendroica_coronata"
## [33] "Dendroica_pinus" "Dromaius_novaehollandiae"
## [35] "Egretta_caerulea" "Egretta_rufescens"
## [37] "Egretta_thula" "Erithacus_rubecula"
## [39] "Eudyptes_chrysocome" "Eudyptes_chrysolophus"
## [41] "Eudyptula_minor" "Falco_mexicanus"
## [43] "Falco_sparverius" "Falco_subbuteo"
## [45] "Ficedula_hypoleuca" "Garrulus_glandarius"
## [47] "Gavia_immer" "Glyphorynchus_spirurus"
## [49] "Ixobrychus_exilis" "Lagopus_lagopus"
## [51] "Lagopus_muta" "Lanius_ludovicianus"
## [53] "Larus_argentatus" "Laterallus_jamaicensis"
## [55] "Lavia_frons" "Leptonychotes_weddellii"
## [57] "Loxia_pytyopsittacus" "Megaceryle_alcyon"
## [59] "Megascops_asio" "Melanerpes_carolinus"
## [61] "Melopsittacus_undulatus" "Melospiza_melodia"
## [63] "Mimus_polyglottos" "Morus_bassanus"
## [65] "Morus_capensis" "Motacilla_alba"
## [67] "Mycteria_americana" "Nyctanassa_violacea"
## [69] "Nycticorax_nycticorax" "Nymphicus_hollandicus"
## [71] "Oceanites_oceanicus" "Opisthocomus_hoazin"
## [73] "Pandion_haliaetus" "Parus_caeruleus"
## [75] "Parus_major" "Passer_domesticus"
## [77] "Passer_montanus" "Pelecanus_occidentalis"
## [79] "Perdix_perdix" "Pernis_apivorus"
## [81] "Phalacrocorax_auritus" "Phoebastria_immutabilis"
## [83] "Phoebetria_fusca" "Phoenicurus_ochruros"
## [85] "Phoenicurus_phoenicurus" "Phylloscopus_collybita"
## [87] "Phylloscopus_trochilus" "Phytotoma_rara"
## [89] "Pinicola_enucleator" "Pluvialis_squatarola"
## [91] "Porphyrio_hochstetteri" "Prunella_modularis"
## [93] "Puffinus_gravis" "Puffinus_lherminieri"
## [95] "Pygoscelis_adeliae" "Pygoscelis_papua"
## [97] "Pyrrhula_pyrrhula" "Regulus_regulus"
## [99] "Rhea_americana" "Sitta_carolinensis"
## [101] "Spheniscus_demersus" "Spheniscus_magellanicus"
## [103] "Struthio_camelus" "Sturnus_vulgaris"
## [105] "Sylvia_atricapilla" "Sylvia_borin"
## [107] "Tetrao_urogallus" "Trichoglossus_haematodus"
## [109] "Troglodytes_aedon" "Troglodytes_troglodytes"
## [111] "Turdus_grayi" "Turdus_merula"
## [113] "Turdus_migratorius" "Turdus_philomelos"
## [115] "Tyrannus_melancholicus" "Tyto_alba"
## [117] "Zenaida_macroura"
##
## $data_not_tree
## character(0)# easy pruner - remove species from the tree that are not in the matrix
prune <- treedata(bigTree, mammalData)## Warning in treedata(bigTree, mammalData): The following tips were not found in 'data' and were dropped from 'phy':
## Accipiter_cooperii
## Accipiter_nisus
## Accipiter_striatus
## Alectoris_graeca
## Anas_platyrhynchos
## Anhinga_anhinga
## Aptenodytes_patagonicus
## Apteryx_australis
## Aramus_guarauna
## Ardea_herodias
## Athene_cunicularia
## Aythya_collaris
## Bombycilla_garrulus
## Botaurus_lentiginosus
## Bubo_virginianus
## Bubulcus_ibis
## Buteo_buteo
## Buteo_jamaicensis
## Buteo_lineatus
## Butorides_striata
## Cacatua_galerita
## Calonectris_diomedea
## Carduelis_carduelis
## Carduelis_chloris
## Casuarius_casuarius
## Cathartes_aura
## Charadrius_vociferus
## Coccothraustes_coccothraustes
## Columba_livia
## Corvus_frugilegus
## Coturnix_japonica
## Dendroica_coronata
## Dendroica_pinus
## Dromaius_novaehollandiae
## Egretta_caerulea
## Egretta_rufescens
## Egretta_thula
## Erithacus_rubecula
## Eudyptes_chrysocome
## Eudyptes_chrysolophus
## Eudyptula_minor
## Falco_mexicanus
## Falco_sparverius
## Falco_subbuteo
## Ficedula_hypoleuca
## Garrulus_glandarius
## Gavia_immer
## Glyphorynchus_spirurus
## Ixobrychus_exilis
## Lagopus_lagopus
## Lagopus_muta
## Lanius_ludovicianus
## Larus_argentatus
## Laterallus_jamaicensis
## Lavia_frons
## Leptonychotes_weddellii
## Loxia_pytyopsittacus
## Megaceryle_alcyon
## Megascops_asio
## Melanerpes_carolinus
## Melopsittacus_undulatus
## Melospiza_melodia
## Mimus_polyglottos
## Morus_bassanus
## Morus_capensis
## Motacilla_alba
## Mycteria_americana
## Nyctanassa_violacea
## Nycticorax_nycticorax
## Nymphicus_hollandicus
## Oceanites_oceanicus
## Opisthocomus_hoazin
## Pandion_haliaetus
## Parus_caeruleus
## Parus_major
## Passer_domesticus
## Passer_montanus
## Pelecanus_occidentalis
## Perdix_perdix
## Pernis_apivorus
## Phalacrocorax_auritus
## Phoebastria_immutabilis
## Phoebetria_fusca
## Phoenicurus_ochruros
## Phoenicurus_phoenicurus
## Phylloscopus_collybita
## Phylloscopus_trochilus
## Phytotoma_rara
## Pinicola_enucleator
## Pluvialis_squatarola
## Porphyrio_hochstetteri
## Prunella_modularis
## Puffinus_gravis
## Puffinus_lherminieri
## Pygoscelis_adeliae
## Pygoscelis_papua
## Pyrrhula_pyrrhula
## Regulus_regulus
## Rhea_americana
## Sitta_carolinensis
## Spheniscus_demersus
## Spheniscus_magellanicus
## Struthio_camelus
## Sturnus_vulgaris
## Sylvia_atricapilla
## Sylvia_borin
## Tetrao_urogallus
## Trichoglossus_haematodus
## Troglodytes_aedon
## Troglodytes_troglodytes
## Turdus_grayi
## Turdus_merula
## Turdus_migratorius
## Turdus_philomelos
## Tyrannus_melancholicus
## Tyto_alba
## Zenaida_macroura# returns list with two items: prune$phy is the tree and prune$dat is the matrix and both match
prunedTree <- prune$phy
name.check(prunedTree, mammalData)## [1] "OK"Standard linear regression NO PHYLO
Let’s test for a relationship between trait1 and trait2 using standard linear regression.
plot(trait2 ~ trait1, data=mammalData)
nonPhyloModel <- lm(trait2 ~ trait1, data=mammalData)
summary(nonPhyloModel)##
## Call:
## lm(formula = trait2 ~ trait1, data = mammalData)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0868 -0.6848 -0.1542 0.6158 3.1560
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.59786 0.22209 -47.72 <2e-16 ***
## trait1 2.05654 0.04056 50.70 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.046 on 196 degrees of freedom
## Multiple R-squared: 0.9292, Adjusted R-squared: 0.9288
## F-statistic: 2571 on 1 and 196 DF, p-value: < 2.2e-16anova(nonPhyloModel)## Analysis of Variance Table
##
## Response: trait2
## Df Sum Sq Mean Sq F value Pr(>F)
## trait1 1 2812.62 2812.62 2570.7 < 2.2e-16 ***
## Residuals 196 214.44 1.09
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(trait2 ~ trait1, data=mammalData)
abline(nonPhyloModel, col="red")
plot(nonPhyloModel)
t1 <- mammalData[,"trait1"]
t2 <- mammalData[,"trait2"]
t1_squared <- t1^2
austinModel <- lm(t2 ~ t1 + t1_squared)
anova(austinModel)## Analysis of Variance Table
##
## Response: t2
## Df Sum Sq Mean Sq F value Pr(>F)
## t1 1 2812.62 2812.62 2661.2608 < 2.2e-16 ***
## t1_squared 1 8.35 8.35 7.9044 0.005435 **
## Residuals 195 206.09 1.06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1predictedValues<-predict(austinModel)
plot(t2~t1)
points(t1, predictedValues, pch=19, cex=0.5, col="red")
require(car)## Loading required package: carAIC(nonPhyloModel)## [1] 583.697AIC(austinModel)## [1] 577.8294# conclusion = austin winsPIC analysis
Now we will use independent contrasts to look for an evolutionary correlation
require(phytools)## Loading required package: phytools## Loading required package: mapsphylomorphospace(prunedTree, mammalData[,c("trait1","trait2")], label="off", node.size=c(0.5, 0.7))
abline(nonPhyloModel, col="red")
t1 <- mammalData[,"trait1"]
t2 <- mammalData[,"trait2"]
# this is super bad news
t1_pic<-pic(t1, prunedTree)
# t1 does not have names = usually get wrong answer
names(t1) <- rownames(mammalData)
names(t2) <- rownames(mammalData)
# this is correct now
t1_pic<-pic(t1, prunedTree)
t2_pic<-pic(t2, prunedTree)
plot(t2_pic~t1_pic)
picModel<- lm(t2_pic ~ t1_pic + 0)
summary(picModel)##
## Call:
## lm(formula = t2_pic ~ t1_pic + 0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.70207 -0.09079 0.00345 0.08953 0.63270
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## t1_pic 1.67632 0.08169 20.52 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1573 on 196 degrees of freedom
## Multiple R-squared: 0.6824, Adjusted R-squared: 0.6808
## F-statistic: 421.1 on 1 and 196 DF, p-value: < 2.2e-16plot(t2_pic~t1_pic)
abline(picModel, col="blue")
rr<-residuals(picModel)
hist(rr)
plot(picModel)
plot(prunedTree, cex=0.2)
nodelabels(cex=0.5)
PGLS version
require(nlme)## Loading required package: nlmepglsModel <- gls(trait2~trait1, correlation=corBrownian(phy=prunedTree), data=mammalData)
summary(pglsModel)## Generalized least squares fit by REML
## Model: trait2 ~ trait1
## Data: mammalData
## AIC BIC logLik
## 520.5517 530.386 -257.2758
##
## Correlation Structure: corBrownian
## Formula: ~1
## Parameter estimate(s):
## numeric(0)
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -8.420445 1.110084 -7.585409 0
## trait1 1.676319 0.081690 20.520486 0
##
## Correlation:
## (Intr)
## trait1 -0.412
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.2285526 -0.5613046 -0.1962622 0.1739419 1.9416041
##
## Residual standard error: 2.02735
## Degrees of freedom: 198 total; 196 residualplot(trait2~trait1, data=mammalData, pch=19, cex=0.5)
abline(nonPhyloModel, col="red")
abline(pglsModel, col="blue")
coefficients(picModel)## t1_pic
## 1.676319coefficients(pglsModel)## (Intercept) trait1
## -8.420445 1.676319# to facilitate model comparison in PGLS force ML instead of REML
pglsModel_ML <- gls(trait2~trait1, correlation=corBrownian(phy=prunedTree), data=mammalData, method="ML")
summary(pglsModel_ML)## Generalized least squares fit by maximum likelihood
## Model: trait2 ~ trait1
## Data: mammalData
## AIC BIC logLik
## 519.23 529.0948 -256.615
##
## Correlation Structure: corBrownian
## Formula: ~1
## Parameter estimate(s):
## numeric(0)
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -8.420445 1.110084 -7.585409 0
## trait1 1.676319 0.081690 20.520486 0
##
## Correlation:
## (Intr)
## trait1 -0.412
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.2348048 -0.5641611 -0.1972610 0.1748272 1.9514851
##
## Residual standard error: 2.017085
## Degrees of freedom: 198 total; 196 residualcoefficients(pglsModel)## (Intercept) trait1
## -8.420445 1.676319coefficients(pglsModel_ML)## (Intercept) trait1
## -8.420445 1.676319pglsLambda <- gls(trait2~trait1, correlation=corPagel(0.5, phy=prunedTree, fixed=FALSE), data=mammalData)
summary(pglsLambda)## Generalized least squares fit by REML
## Model: trait2 ~ trait1
## Data: mammalData
## AIC BIC logLik
## 484.9451 498.0576 -238.4726
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.9118228
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -9.699775 0.8114129 -11.95418 0
## trait1 1.905912 0.0740676 25.73208 0
##
## Correlation:
## (Intr)
## trait1 -0.511
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.9025900 -0.5730860 -0.1616109 0.3034165 2.3799271
##
## Residual standard error: 1.441062
## Degrees of freedom: 198 total; 196 residualpglsModel_2 <- gls(trait2~trait1*Lifestyle, correlation=corBrownian(phy=prunedTree), data=mammalData)
summary(pglsModel_2)## Generalized least squares fit by REML
## Model: trait2 ~ trait1 * Lifestyle
## Data: mammalData
## AIC BIC logLik
## 524.6794 547.4819 -255.3397
##
## Correlation Structure: corBrownian
## Formula: ~1
## Parameter estimate(s):
## numeric(0)
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) -6.311608 1.5748667 -4.007709 0.0001
## trait1 1.418196 0.1709287 8.297006 0.0000
## LifestyleFlying -3.330016 1.6920532 -1.968032 0.0505
## LifestyleTerrestrial -2.211888 1.4406688 -1.535321 0.1264
## trait1:LifestyleFlying 0.329312 0.2749398 1.197761 0.2325
## trait1:LifestyleTerrestrial 0.248514 0.2015801 1.232830 0.2191
##
## Correlation:
## (Intr) trait1 LfstyF LfstyT tr1:LF
## trait1 -0.741
## LifestyleFlying -0.599 0.737
## LifestyleTerrestrial -0.674 0.830 0.657
## trait1:LifestyleFlying 0.461 -0.622 -0.765 -0.516
## trait1:LifestyleTerrestrial 0.607 -0.816 -0.588 -0.952 0.507
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.15884968 -0.41946029 0.03558324 0.37993492 2.03146888
##
## Residual standard error: 2.025188
## Degrees of freedom: 198 total; 192 residualanova(pglsModel_2)## Denom. DF: 192
## numDF F-value p-value
## (Intercept) 1 0.9335 0.3352
## trait1 1 421.9900 <.0001
## Lifestyle 2 1.2288 0.2949
## trait1:Lifestyle 2 0.9806 0.3770