library(ape)
library(ctv)
library(maps)
library(phytools)
library(ade4)
library(adephylo)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
BEAST on XDESE (2.1-2.4.8) ran on CSIPES server
500,000,000 generations
4 chains
116 terminal taxa
Trees sampled every 1000 generations
4 concatenated gene regions with a total of 4072 characters
tree <- read.tree("../Trich_1Milgen_Bayes_Contree")
Species <- tree$tip.label
write.csv(Species, file = "Trich_Species.csv")
Use Species column as row names.
Brooding_Data <- read.csv("Trich_Brooding_Data.csv", row.names = 1)
Transform the dataframe to a factor with 2 levels “no/yes” and species as attr.
## Megadromus_bullatus Megadromus_guerinii Megadromus_antarcticus
## Yes Yes Yes
## Megadromus_sp Megadromus_enyesi Megadromus_walkeri
## Yes Yes Yes
## Levels: No Yes
Final character states mapped on to terminal taxa.
plotTree(tree,fsize=0.4,type = "fan", ftype="i",lwd=1, offset = 1)
cols<-setNames(c("red","blue"),c("Does not Brood", "Broods"))
tiplabels(pie=to.matrix(brood.mode[tree$tip.label],
levels(brood.mode)),piecol=cols,cex=0.3)
add.simmap.legend(colors = cols,prompt=FALSE,x=0.9*par()$usr[1],
y=0.8*par()$usr[3],fsize=0.8)
Since this is a binary character, the Equal Rates model is used which has a single parameter, setting equal forward and reverse rates.
fit <- rerootingMethod(tree, brood.mode, model="ER")
fit
## Ancestral character estimates using re-rooting method
## of Yang et al. (1995):
## No Yes
## 117 0.997645 0.002355
## 118 0.999922 0.000078
## 119 0.999942 0.000058
## 120 0.999901 0.000099
## 121 0.999726 0.000274
## 122 0.979697 0.020303
## 123 0.016271 0.983729
## 124 0.000137 0.999863
## 125 0.000004 0.999996
## 126 0.000000 1.000000
## 127 0.000000 1.000000
## 128 0.998914 0.001086
## 129 0.999973 0.000027
## 130 0.999995 0.000005
## 131 0.999996 0.000004
## 132 0.999928 0.000072
## 133 0.992118 0.007882
## 134 0.995633 0.004367
## 135 0.989729 0.010271
## 136 0.999599 0.000401
## 137 0.999880 0.000120
## 138 0.999968 0.000032
## 139 0.999996 0.000004
## 140 1.000000 0.000000
## 141 0.999995 0.000005
## 142 0.999999 0.000001
## 143 0.999998 0.000002
## 144 0.031849 0.968151
## 145 0.000088 0.999912
## 146 0.000000 1.000000
## 147 0.000000 1.000000
## 148 0.000000 1.000000
## 149 0.000000 1.000000
## 150 0.011099 0.988901
## 151 0.000318 0.999682
## 152 0.000001 0.999999
## 153 0.000000 1.000000
## 154 0.000000 1.000000
## 155 0.000000 1.000000
## 156 0.000000 1.000000
## 157 0.000003 0.999997
## 158 0.000001 0.999999
## 159 0.000000 1.000000
## 160 0.000000 1.000000
## 161 0.000000 1.000000
## 162 0.000000 1.000000
## 163 0.000001 0.999999
## 164 0.024177 0.975823
## 165 0.000002 0.999998
## 166 0.000000 1.000000
## 167 0.000000 1.000000
## 168 0.000000 1.000000
## 169 0.000000 1.000000
## 170 0.999664 0.000336
## 171 0.999973 0.000027
## 172 0.999998 0.000002
## 173 0.999999 0.000001
## 174 0.999991 0.000009
## 175 0.000038 0.999962
## 176 0.000001 0.999999
## 177 0.000016 0.999984
## 178 0.000008 0.999992
## 179 0.000000 1.000000
## 180 0.000000 1.000000
## 181 0.000008 0.999992
## 182 0.000013 0.999987
## 183 0.000000 1.000000
## 184 0.000000 1.000000
## 185 0.000000 1.000000
## 186 0.000000 1.000000
## 187 0.000000 1.000000
## 188 0.000000 1.000000
## 189 0.000001 0.999999
## 190 0.999982 0.000018
## 191 0.999985 0.000015
## 192 0.999983 0.000017
## 193 0.999986 0.000014
## 194 0.999999 0.000001
## 195 0.999999 0.000001
## 196 0.999980 0.000020
## 197 0.999997 0.000003
## 198 1.000000 0.000000
## 199 0.999997 0.000003
## 200 0.999997 0.000003
## 201 0.999993 0.000007
## 202 0.999999 0.000001
## 203 0.999999 0.000001
## 204 1.000000 0.000000
## 205 0.999984 0.000016
## 206 0.999988 0.000012
## 207 1.000000 0.000000
## 208 0.999993 0.000007
## 209 0.999968 0.000032
## 210 0.999968 0.000032
## 211 0.999963 0.000037
## 212 0.999981 0.000019
## 213 1.000000 0.000000
## 214 1.000000 0.000000
## 215 0.999939 0.000061
## 216 0.999994 0.000006
## 217 0.999982 0.000018
## 218 0.999992 0.000008
## 219 0.999992 0.000008
## 220 0.999982 0.000018
## 221 1.000000 0.000000
## 222 0.999806 0.000194
## 223 0.999448 0.000552
##
## Estimated transition matrix,
## Q =
## No Yes
## No -1.922069 1.922069
## Yes 1.922069 -1.922069
##
## **Note that if Q is not symmetric the marginal
## reconstructions may be invalid.
##
## Log-likelihood = -18.955108
Internal nodes represent state probabilities/
Tips represent the final state.
plotTree(tree, fsize=0.4, ftype="i",lwd=1, offset = 1)
cols<-setNames(c("red","blue"),c("Does not Brood", "Broods"))
nodelabels(pie = fit$marginal.anc, piecol = cols, cex=0.3)
tiplabels(pie=to.matrix(brood.mode[tree$tip.label],
levels(brood.mode)),piecol=cols,cex=0.2)
add.simmap.legend(colors = cols,prompt=FALSE,x=0.9*par()$usr[1],
y=6,fsize=0.8)
Stochastic character mapping?