-
Notifications
You must be signed in to change notification settings - Fork 0
/
continuous_anc.Rmd
71 lines (48 loc) · 2.31 KB
/
continuous_anc.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
---
title: "Continuous Ancestral State"
output: github_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
#devtools::install_github("arborworkflows/aRbor")
library(aRbor)
library(geiger)
library(knitr)
```
## Loading in sample data and running function
```{r}
tr <- read.tree("/Users/carrietribble/Desktop/arbor_hackathon_docs/heliconia/Heliconia_all_dated_mod.phy")
trait <- read.csv('/Users/carrietribble/Desktop/arbor_hackathon_docs/heliconia/Heliconia_phylogeny_matrix.csv')
td <- treeplyr::make.treedata(tr, trait)
ace_arbor <- aRbor::aceArbor( select_(td, "InflorL"), charType='continuous', aceType = 'marginal')
```
## Summary of Continuous Ancestral State Reconstruction
```{r include = FALSE}
trait_name <- attributes(ace_arbor)$names
orig_mean <- mean(as.matrix(attributes(ace_arbor)$td$dat))
anc_mean <- mean(ace_arbor$InflorL$estimate)
orig_range <- c(min(as.matrix(attributes(ace_arbor)$td$dat)),
max(as.matrix(attributes(ace_arbor)$td$dat)))
anc_range <- c(min(ace_arbor[[1]]$estimate),
max(ace_arbor[[1]]$estimate))
anc_range_CI <- c(min(ace_arbor[[1]]$lowerCI95),
max(ace_arbor[[1]]$upperCI95))
# ordering and binning confidence intervals by age
CI_age <- data.frame(age = branching.times(attributes(ace_arbor)$td$phy),
CI_width = ace_arbor[[1]]$upperCI95 - ace_arbor[[1]]$lowerCI95)
```
The mean value of `r paste(trait_name)` extant values is `r orig_mean`, the reconstructed nodes have a mean value of `r anc_mean`. The original range of your data values is `r orig_range[1]` to `r orig_range[2]`. The estimated node values range from `r anc_range[1]` to `r anc_range[2]`. When we include the 95% confidence invervals for those estimates, the estimates range from `r anc_range_CI[1]` to `r anc_range_CI[2]`. A portion of those estimates are displayed below.
```{r, echo = FALSE}
if (nrow(ace_arbor[[1]]) > 20) {
rows_to_include <- 20
} else {rows_to_include <- nrow(ace_arbor[[1]])}
kable(ace_arbor[[1]][1:rows_to_include,])
```
The plot below shows the relationship between the width of your 95% confidence intervals and time as indicated by node ages.
```{r echo = FALSE}
plot(CI_width~age, data = CI_age,
pch = 19,
main = trait_name,
xlab = "Node Age",
ylab = "Width of 95% CI for Estimated Values")
```