corHMM_Diet.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Testing and optimizing transition between states"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "library(\"corHMM\")\n",
    "library(\"mclust\")\n",
    "library(\"stringr\")\n",
    "library(\"qpcR\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading data\n",
    "phy<-read.nexus(\"16FC_16C_374_sp.tree\")\n",
    "df<-read.table(\"table_data_diet.tsv\", header = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cleaning and preparing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df$Species<-str_replace((df$Species), \" \", \"_\")\n",
    "\n",
    "states<-cbind(c(str_replace((df$Species), \" \", \"_\"), setdiff(phy$tip.label, str_replace((df$Species), \" \", \"_\"))), c(df$traits, rep(\"?\", length(setdiff(phy$tip.label, str_replace((df$Species), \" \", \"_\"))))))\n",
    "\n",
    "states_traits<-as.data.frame(states[!states[,1] %in% setdiff(states[,1], phy$tip.label),])\n",
    "\n",
    "colnames(states_traits)<-c(\"Species\", \"Diet\")\n",
    "\n",
    "states_traits<-states_traits[match(phy$tip.label,states_traits[,1]),]\n",
    "\n",
    "states_traits[is.na(states_traits)]<-\"?\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing the evolution of diet across the tree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cleaning data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "states_traits<-states_traits[,c(1,2)]\n",
    "LegendAndRateMat <- getStateMat4Dat(states_traits)\n",
    "RateMat <- LegendAndRateMat$rate.mat\n",
    "RateMat_trans <- dropStateMatPars(RateMat, c(1,3))\n",
    "pars2equal <- list(c(1,3), c(2,4))\n",
    "RateMat_trans_sym <- equateStateMatPars(RateMat_trans, pars2equal)\n",
    "pars2equal <- list(c(1:4))\n",
    "RateMat_trans_eq <- equateStateMatPars(RateMat_trans, pars2equal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One rate models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 1 : equal rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_1_eq<-corHMM(phy, states_traits, rate.cat = 1, rate.mat=NULL, model = \"ER\", node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 2 symetric rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_1_sym<-corHMM(phy, states_traits, rate.cat = 1, rate.mat=NULL, model = \"SYM\", node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 3  all rates differ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_1_ard<-corHMM(phy, states_traits, rate.cat = 1, rate.mat=NULL, model = \"ARD\", node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 4 : no transition between Invertebrate-feeders and Macropredator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_1_tran_I_M<-corHMM(phy, states_traits, rate.cat = 1, rate.mat=RateMat_trans, node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 5 : no transition between Invertebrate-feeders and Macropredator, symetric rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_1_tran_I_M_sym<-corHMM(phy, states_traits, rate.cat = 1, rate.mat=RateMat_trans_sym, node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 6 : no transition between Invertebrate-feeders and Macropredator, equal rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_1_tran_I_M_eq<-corHMM(phy, states_traits, rate.cat = 1, rate.mat=RateMat_trans_eq, node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Two rates models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RateMat_trans_2_rates<-list(RateMat_trans, RateMat_trans)\n",
    "RateMat_trans_sym_2_rates<-list(RateMat_trans_sym, RateMat_trans_sym)\n",
    "RateMat_trans_eq_2_rates<-list(RateMat_trans_eq, RateMat_trans_eq)\n",
    "RateClassMat <- getRateCatMat(2) \n",
    "RateMat_trans_2_rates <- getFullMat(RateMat_trans_2_rates, RateClassMat)\n",
    "RateMat_trans_sym_2_rates <- getFullMat(RateMat_trans_sym_2_rates, RateClassMat)\n",
    "RateMat_trans_eq_2_rates <- getFullMat(RateMat_trans_eq_2_rates, RateClassMat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 1 : equal rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_2_eq<-corHMM(phy, states_traits, rate.cat = 2, rate.mat=NULL, model = \"ER\", node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 2 : symmetric rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_2_sym<-corHMM(phy, states_traits, rate.cat = 2, rate.mat=NULL, model = \"SYM\", node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 3 : all rates differ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_2_ard<-corHMM(phy, states_traits, rate.cat = 2, rate.mat=NULL, model = \"ARD\", node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 4 : no transition between Invertebrate-feeders and Macropredator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_2_tran_I_M<-corHMM(phy, states_traits, rate.cat = 2, rate.mat= RateMat_trans_2_rates, node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 5 : no transition between Invertebrate-feeders and Macropredator symmetric rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_2_tran_I_M_sym<-corHMM(phy, states_traits, rate.cat = 2, rate.mat=RateMat_trans_sym_2_rates, node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 6 : no transition between Invertebrate-feeders and Macropredator equal rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_2_tran_I_M_eq<-corHMM(phy, states_traits, rate.cat = 2, rate.mat=RateMat_trans_eq_2_rates, node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Three rates models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RateMat_trans_3_rates<-list(RateMat_trans, RateMat_trans, RateMat_trans)\n",
    "RateMat_trans_sym_3_rates<-list(RateMat_trans_sym, RateMat_trans_sym, RateMat_trans_sym)\n",
    "RateMat_trans_eq_3_rates<-list(RateMat_trans_eq, RateMat_trans_eq, RateMat_trans_eq)\n",
    "RateClassMat <- getRateCatMat(3) \n",
    "RateMat_trans_3_rates <- getFullMat(RateMat_trans_3_rates, RateClassMat)\n",
    "RateMat_trans_sym_3_rates <- getFullMat(RateMat_trans_sym_3_rates, RateClassMat)\n",
    "RateMat_trans_eq_3_rates <- getFullMat(RateMat_trans_eq_3_rates, RateClassMat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 1 : equal rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_3_eq<-corHMM(phy, states_traits, rate.cat = 3, rate.mat=NULL, model = \"ER\", node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 2 : symmetric rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_3_sym<-corHMM(phy, states_traits, rate.cat = 3, rate.mat=NULL, model = \"SYM\", node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 3 : all rates differ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_3_ard<-corHMM(phy, states_traits, rate.cat = 3, rate.mat=NULL, model = \"ARD\", node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 4  : no transition between Invertebrate-feeders and Macropredator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "diet_3_tran_I_M<-corHMM(phy, states_traits, rate.cat = 3, rate.mat=RateMat_trans_3_rates, node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 5 : no transition between Invertebrate-feeders and Macropredator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_3_tran_I_M_sym<-corHMM(phy, states_traits, rate.cat = 3, rate.mat=RateMat_trans_sym_3_rates, node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model 6 : no transition between Invertebrate-feeders and Macropredator equal rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diet_3_tran_I_M_eq<-corHMM(phy, states_traits, rate.cat = 3, rate.mat=RateMat_trans_eq_3_rates, node.states = \"marginal\",\n",
    "fixed.nodes=FALSE, p=NULL, root.p=\"yang\", ip=NULL, nstarts=0, n.cores=5,\n",
    "get.tip.states = FALSE, lewis.asc.bias = FALSE, collapse = TRUE, lower.bound = 1e-9,\n",
    "upper.bound = 100, opts=NULL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assessing the best-fit model across all combinations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_diet<-data.frame(cbind(c(diet_1_eq$loglik, diet_1_sym$loglik, diet_1_ard$loglik, diet_1_tran_I_M$loglik, diet_1_tran_I_M_sym$loglik, diet_1_tran_I_M_eq$loglik,\n",
    "                          diet_2_eq$loglik, diet_2_sym$loglik, diet_2_ard$loglik, diet_2_tran_I_M$loglik, diet_2_tran_I_M_sym$loglik, diet_2_tran_I_M_eq$loglik,\n",
    "                          diet_3_eq$loglik, diet_3_sym$loglik, diet_3_ard$loglik, diet_3_tran_I_M$loglik, diet_3_tran_I_M_sym$loglik, diet_3_tran_I_M_eq$loglik),\n",
    "                        c(diet_1_eq$AICc, diet_1_sym$AICc, diet_1_ard$AICc, diet_1_tran_I_M$AICc, diet_1_tran_I_M_sym$AICc, diet_1_tran_I_M_eq$AICc,\n",
    "                          diet_2_eq$AICc, diet_2_sym$AICc, diet_2_ard$AICc, diet_2_tran_I_M$AICc, diet_2_tran_I_M_sym$AICc, diet_2_tran_I_M_eq$AICc,\n",
    "                          diet_3_eq$AICc, diet_3_sym$AICc, diet_3_ard$AICc, diet_3_tran_I_M$AICc, diet_3_tran_I_M_sym$AICc, diet_3_tran_I_M_eq$AICc),\n",
    "                akaike.weights(c(diet_1_eq$AICc, diet_1_sym$AICc, diet_1_ard$AICc, diet_1_tran_I_M$AICc, diet_1_tran_I_M_sym$AICc, diet_1_tran_I_M_eq$AICc,\n",
    "                          diet_2_eq$AICc, diet_2_sym$AICc, diet_2_ard$AICc, diet_2_tran_I_M$AICc, diet_2_tran_I_M_sym$AICc, diet_2_tran_I_M_eq$AICc,\n",
    "                          diet_3_eq$AICc, diet_3_sym$AICc, diet_3_ard$AICc, diet_3_tran_I_M$AICc, diet_3_tran_I_M_sym$AICc, diet_3_tran_I_M_eq$AICc))$deltaAIC,\n",
    "                akaike.weights(c(diet_1_eq$AICc, diet_1_sym$AICc, diet_1_ard$AICc, diet_1_tran_I_M$AICc, diet_1_tran_I_M_sym$AICc, diet_1_tran_I_M_eq$AICc,\n",
    "                          diet_2_eq$AICc, diet_2_sym$AICc, diet_2_ard$AICc, diet_2_tran_I_M$AICc, diet_2_tran_I_M_sym$AICc, diet_2_tran_I_M_eq$AICc,\n",
    "                          diet_3_eq$AICc, diet_3_sym$AICc, diet_3_ard$AICc, diet_3_tran_I_M$AICc, diet_3_tran_I_M_sym$AICc, diet_3_tran_I_M_eq$AICc))$weights,\n",
    "                c((max(as.vector(diet_1_eq$index.mat)[!is.na(as.vector(diet_1_eq$index.mat))])), (max(as.vector(diet_1_sym$index.mat)[!is.na(as.vector(diet_1_sym$index.mat))])), (max(as.vector(diet_1_ard$index.mat)[!is.na(as.vector(diet_1_ard$index.mat))])), (max(as.vector(diet_1_tran_I_M$index.mat)[!is.na(as.vector(diet_1_tran_I_M$index.mat))])), (max(as.vector(diet_1_tran_I_M_sym$index.mat)[!is.na(as.vector(diet_1_tran_I_M_sym$index.mat))])), (max(as.vector(diet_1_tran_I_M_eq$index.mat)[!is.na(as.vector(diet_1_tran_I_M_eq$index.mat))])),\n",
    "(max(as.vector(diet_2_eq$index.mat)[!is.na(as.vector(diet_2_eq$index.mat))])), (max(as.vector(diet_2_sym$index.mat)[!is.na(as.vector(diet_2_sym$index.mat))])), (max(as.vector(diet_2_ard$index.mat)[!is.na(as.vector(diet_2_ard$index.mat))])), (max(as.vector(diet_2_tran_I_M$index.mat)[!is.na(as.vector(diet_2_tran_I_M$index.mat))])), (max(as.vector(diet_2_tran_I_M_sym$index.mat)[!is.na(as.vector(diet_2_tran_I_M_sym$index.mat))])), (max(as.vector(diet_2_tran_I_M_eq$index.mat)[!is.na(as.vector(diet_2_tran_I_M_eq$index.mat))])),\n",
    "(max(as.vector(diet_3_eq$index.mat)[!is.na(as.vector(diet_3_eq$index.mat))])), (max(as.vector(diet_3_sym$index.mat)[!is.na(as.vector(diet_3_sym$index.mat))])), (max(as.vector(diet_3_ard$index.mat)[!is.na(as.vector(diet_3_ard$index.mat))])), (max(as.vector(diet_3_tran_I_M$index.mat)[!is.na(as.vector(diet_3_tran_I_M$index.mat))])), (max(as.vector(diet_3_tran_I_M_sym$index.mat)[!is.na(as.vector(diet_3_tran_I_M_sym$index.mat))])), (max(as.vector(diet_3_tran_I_M_eq$index.mat)[!is.na(as.vector(diet_3_tran_I_M_eq$index.mat))])))\n",
    "                ))\n",
    "rownames(df_diet)<-c(\"diet_1_eq\", \"diet_1_sym\", \"diet_1_ard\", \"diet_1_tran_I_M\", \"diet_1_tran_I_M_sym\", \"diet_1_tran_I_M_eq\",\n",
    "                          \"diet_2_eq\", \"diet_2_sym\", \"diet_2_ard\", \"diet_2_tran_I_M\", \"diet_2_tran_I_M_sym\", \"diet_2_tran_I_M_eq\",\n",
    "                          \"diet_3_eq\", \"diet_3_sym\", \"diet_3_ard\", \"diet_3_tran_I_M\", \"diet_3_tran_I_M_sym\", \"diet_3_tran_I_M_eq\")\n",
    "colnames(df_diet)<-c(\"loglik\", \"AICc\", \"Delta_AICc\", \"AICcWt\", \"K_rates\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "write.table(df_diet, \"Results/df_diet.tsv\", sep =\"\\t\")\n",
    "saveRDS(eval(parse(text = rownames(df_diet)[which.min(df_diet$AICc)])), \"ASE/diet_corHMM.rds\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.3.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}