comparative.rmd

# Comparative methods

## Correlated evolutionary changes

PGLS [@symonds2014primer].


```{r}

# Scale edge lengths so that root has specified age
scale_phy = function(phy, age=1){
  node_root = phy$edge[(!phy$edge[,1] %in% phy$edge[,2]),1] %>% unique()
  age_root = max( dist.nodes(phy)[,node_root] )
  phy$edge.length = phy$edge.length / age_root * age
  phy
}

# Transform edge lengths to squish length into root edges or other edges.
# Also rescaled age of root to 1
# x = 1 all length in tip edges
# x = 0 no transform
# x = -1 all length in root edges
squish_phy = function( phy, x ){
  node_root = phy$edge[(!phy$edge[,1] %in% phy$edge[,2]),1] %>% unique()
  n_tips = length(phy$tip.label)
  
  mult_deep = 1 - x
  mult_shallow = 1 + x   

  edges_root = which(phy$edge[,1]==node_root)
  phy$edge.length[edges_root] = phy$edge.length[edges_root] * mult_deep
  
  edges_tip = which( phy$edge[,2] <= n_tips )
  phy$edge.length[edges_tip] = phy$edge.length[edges_tip] * mult_shallow
  
  edges = 1:nrow(phy$edge)
  edges_mid = which( ! edges %in% c(edges_root, edges_tip)  )
  
  if(x>0){
    phy$edge.length[edges_mid] = phy$edge.length[edges_mid] * mult_deep
  } else if (x<0){
    phy$edge.length[edges_mid] = phy$edge.length[edges_mid] * mult_shallow
  }

  
  tip_ages = dist.nodes(phy)[1:n_tips,node_root]
  trim = max(tip_ages) - min(tip_ages)
  tips_long = which(tip_ages > min(tip_ages))
  tips_edges_long = phy$edge[,2] %in% tips_long
  phy$edge.length[tips_edges_long] = phy$edge.length[tips_edges_long] - trim
  
  
  # Scale root age to one
  phy = scale_phy(phy)
  
  return(phy)
}

tree_text = "((((A:1,B:1):1,(C:1,D:1):1):1,((E:1,F:1):1,(G:1,H:1):1):1):1,(((I:1,J:1):1,(K:1,L:1):1):1,((M:1,N:1):1,(O:1,P:1):1):1):1);"
phy = read.tree(text=tree_text)
phy = scale_phy( phy )
# plot( phy )

phy_star = squish_phy(phy, 1)
phy_deep = squish_phy(phy, -0.98 )

label_offset = 0.05
tree_width = 1.5
ggtree_deep = ggtree( phy_deep ) + 
  geom_tiplab( fontface = "italic", offset=label_offset ) +
  xlim(0, tree_width) +
  geom_text2(aes(label=node), col="red", nudge_x=label_offset/2 )

plot(phy_star)
plot(phy_deep)

```



```{r}
set.seed(123456)

cov_matrix <- function(diag_term=1, cov_term=0) {
  G=2
  trueCovariance = matrix(cov_term,G,G)
  diag(trueCovariance) = diag_term

  trueCovariance
}

plot_matrix <- function(m, ... ) {
	# Basic plotting
	nr <- nrow(m)
	nc <- ncol(m)
	image(1:nc, 1:nr, t(m[nr:1, ]), axes=F,xlab="", ylab="", ... )
}

a_off_diag = 0
b_off_diag = 0.5
c_off_diag = 1
cov_a = cov_matrix(1, a_off_diag)
cov_b = cov_matrix(1, b_off_diag)
cov_c = cov_matrix(1, c_off_diag)

cov_a
cov_b
cov_c
```
## Inferring covariance in the absence of phylogenetic structure


```{r}

sim_char = function( phy, cov_sim ){
  x = sim.char(phy, cov_sim, nsim = 1, model = "BM", root = 0)
  
  # Convert array to matrix
  n_rows = dim(x)[1]
  x = matrix(x, nrow = n_rows)
  
  x
}

Z_star_a = sim_char(phy_star, cov_a)
Z_star_b = sim_char(phy_star, cov_b)
Z_star_c = sim_char(phy_star, cov_c)

Z_deep_a = sim_char(phy_deep, cov_a)
Z_deep_b = sim_char(phy_deep, cov_b)
Z_deep_c = sim_char(phy_deep, cov_c)


D = rbind(
  data.frame(Z_star_a) %>% mutate(tree="star", cov = a_off_diag, node=1:nrow(Z_star_a), name = phy$tip.label ),
  data.frame(Z_star_b) %>% mutate(tree="star", cov = b_off_diag, node=1:nrow(Z_star_b), name = phy$tip.label ),
  data.frame(Z_star_c) %>% mutate(tree="star", cov = c_off_diag, node=1:nrow(Z_star_c), name = phy$tip.label ),
  data.frame(Z_deep_a) %>% mutate(tree="deep", cov = a_off_diag, node=1:nrow(Z_deep_a), name = phy$tip.label ),
  data.frame(Z_deep_b) %>% mutate(tree="deep", cov = b_off_diag, node=1:nrow(Z_deep_b), name = phy$tip.label ),
  data.frame(Z_deep_c) %>% mutate(tree="deep", cov = c_off_diag, node=1:nrow(Z_deep_c), name = phy$tip.label )
)




D %>% 
  filter(cov==0) %>%
  ggplot(aes(x=X1, y=X2, color=cov )) +
  #geom_text( aes(label=name), color="red", vjust = 0, nudge_y = 0.01 ) +
  geom_text_repel(aes(label=name), color="black") +
  geom_point() +
  facet_wrap(~tree)


```

## Phylogenetic Independent Contrasts

```{r}

  

```



## Phylogenetic structure is the expected covariance structure

```{r}

```