geiger versus multivariate normal

geiger.R

# Confirming that I understand how geiger's variance-covariance matrices work

library(geiger)
library(mvtnorm)

geo = get(data("geospiza"))


# Simulate one trait using geiger::sim.char
sim = t(
  sim.char(
    drop.tip(geo$phy, "olivacea"), 
    matrix(1),
    nsim = 1E5,
    root = 0
  )[ , 1, ]
)


# Simulate one trait using multivariate normal
sim2 = rmvnorm(1E5, sigma = vcv(drop.tip(geo$phy, "olivacea")))

# Simulate again to assess sampling error
sim3 = rmvnorm(1E5, sigma = vcv(drop.tip(geo$phy, "olivacea")))

plot(
  round(cov(sim), 3),
  round(cov(sim2), 3)
)
abline(0,1)


plot(
  round(cov(sim2), 3),
  round(cov(sim3), 3)
)
abline(0,1)