7.6: R Markdown para recrear análisis

Última actualización
Guardar como PDF

Page ID: 54479

\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \) \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)\(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\) \(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\)\(\newcommand{\AA}{\unicode[.8,0]{x212B}}\)

R rebajas para recrear análisis

Lectura en los archivos de datos

Primero leemos en los archivos de datos.

sqTree<-read.tree(text=getURL("https://raw.githubusercontent.com/lukejharmon/pcm/master/datafiles/squamate.phy"))
plot(sqTree)

sqData<-read.csv(text=getURL("https://raw.githubusercontent.com/lukejharmon/pcm/master/datafiles/brandley_table.csv"))

Simular carácter binario en el árbol

Este código genera gráficas como la Figura 7.4

qMatrix<-cbind(c(-1, 1), c(1, -1))*0.001
sh_slow<-sim.history(sqTree, qMatrix, anc="1")

## Done simulation(s).

plotSimmap(sh_slow, pts=F, ftype="off")

## no colors provided. using the following legend:
##       1       2 
## "black"   "red"

add.simmap.legend(leg=c("limbed", "limbless"), colors=c("black", "red"), x=0.024, y =23, prompt=F)

qMatrix<-cbind(c(-1, 1), c(1, -1))*0.01
sh_fast<-sim.history(sqTree, qMatrix, anc="1")

## Done simulation(s).

plotSimmap(sh_fast, pts=F, ftype="off")

## no colors provided. using the following legend:
##       1       2 
## "black"   "red"

qMatrix<-cbind(c(-0.02, 0.02), c(0.005, -0.005))
sh_asy<-sim.history(sqTree, qMatrix, anc="1")

## Note - the rate of substitution from i->j should be given by Q[j,i].
## Done simulation(s).

plotSimmap(sh_asy, pts=F, ftype="off")

## no colors provided. using the following legend:
##       1       2 
## "black"   "red"

Encuentra las especies sin extremidades

Los datos de Brandley et al. tienen mediciones de extremidades. Obtendremos nuestro carácter discreto contando especies con extremidades anteriores y traseras de longitud cero como sin fisuras. Esto es diferente del análisis original de Brandley et al., que cuenta cosas como espolones como “extremidades”, por lo que nuestros resultados pueden diferir un poco de los suyos.

limbless<-as.numeric(sqData[,"FLL"]==0 & sqData[,"HLL"]==0)
sum(limbless)

## [1] 51

# get names that match
nn<-sqData[,1]
nn2<-sub(" ", "_", nn)
names(limbless)<-nn2

Modelo Fit Mk

Podemos ajustar un modelo Mk simétrico a estos datos usando tanto verosimilitud como MCMC

# likelihood
td<-treedata(sqTree, limbless)

## Warning in treedata(sqTree, limbless): The following tips were not found in 'phy' and were dropped from 'data':
##  Gonatodes_albogularis
##  Lepidophyma_flavimaculatum
##  Trachyboa_boulengeri

dModel<-fitDiscrete(td$phy, td$data)

# MCMC
mk_diversitree<-make.mk2(force.ultrametric(td$phy), td$data[,1])
simplemk<-constrain(mk_diversitree, q01~q10)
er_bayes<-mcmc(simplemk, x.init=0.1, nsteps=10000, w=0.01)

## 1: {0.0189} -> -141.93605
## 2: {0.0162} -> -135.45732