SecSSE_ASE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ancestral state estimation SecSSE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Loading packages\n",
    "library(\"mclust\")\n",
    "library(\"Rtsne\")\n",
    "library(\"ape\")\n",
    "library(\"secsse\")\n",
    "library(\"DDD\")\n",
    "library(\"doMC\")\n",
    "library(\"tidyverse\")\n",
    "library(\"qgraph\")\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",
    "\n",
    "df<-read.csv(\"Species_data.tsv\", sep=\"\\t\") # omit sep =\"\\t\" for .csv files\n",
    "\n",
    "phy<-read.nexus(\"16FC_16C_374_sp.tree\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating AICc function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aicc<-function(ll,k){round((2*(-round(ll,4))+2*k+(2*k*(k+1))/(374-k-1)),3)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing column names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clnames<-c(\"Log Likelihood\", \"Parameters\", \"AICc\", \"Delta AICc\", \"AICcWt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximum Body-size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mb = Mclust(as.numeric(df$Body.size))\n",
    "summary(mb, parameters = TRUE)\n",
    "df$Body.size <- as.factor(mb$classification)\n",
    "\n",
    "# Optional step (cleaning data)\n",
    "states<-df$Body.size\n",
    "names(states)=str_replace((df$Species), \" \", \"_\")\n",
    "states2<-states[!names(states) %in% setdiff(names(states), phy$tip.label)]\n",
    "\n",
    "states2<-states[!names(states) %in% setdiff(names(states), phy$tip.label)]\n",
    "\n",
    "states3<-as.data.frame(cbind(names(states2), as.factor(states2)))\n",
    "colnames(states3)<-c(\"species\",\"states\")\n",
    "rownames(states3)<-NULL\n",
    "traits <- sortingtraits(states3, phy)\n",
    "traits<-traits-1\n",
    "\n",
    "# Computing sampling fractions\n",
    "f<-c(\n",
    "length(states2[states2==\"1\"])/length(df$Body.size[df$Body.size==\"1\"]),\n",
    "length(states2[states2==\"2\"])/length(df$Body.size[df$Body.size==\"2\"]),\n",
    "length(states2[states2==\"3\"])/length(df$Body.size[df$Body.size==\"3\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementing model structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_matrix_cr <- c()\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(0, 0, 0, 1))\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(1, 1, 1, 1))\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(2, 2, 2, 1))\n",
    "\n",
    "lambda_list_cr <- secsse::create_lambda_list(state_names = c(0, 1, 2),\n",
    "                                          num_concealed_states = 3,\n",
    "                                          transition_matrix = spec_matrix_cr,\n",
    "                                          model = \"CR\")\n",
    "mu_vec_cr <- secsse::create_mu_vector(state_names = c(0, 1, 2),\n",
    "                                   num_concealed_states = 3,\n",
    "                                   model = \"CR\",\n",
    "                                   lambda_list = lambda_list_cr)\n",
    "\n",
    "shift_matrix_cr  <- c()\n",
    "shift_matrix_cr  <- rbind(shift_matrix_cr , c(0, 1, 3))\n",
    "shift_matrix_cr  <- rbind(shift_matrix_cr , c(1, 0, 3))\n",
    "shift_matrix_cr  <- rbind(shift_matrix_cr , c(1, 2, 4))\n",
    "shift_matrix_cr  <- rbind(shift_matrix_cr , c(2, 1, 4))\n",
    "\n",
    "q_matrix_cr <- secsse::create_q_matrix(state_names = c(0, 1, 2),\n",
    "                                    num_concealed_states = 3,\n",
    "                                    shift_matrix = shift_matrix_cr,\n",
    "                                    diff.conceal = TRUE)\n",
    "q_matrix_cr[q_matrix_cr>=5]<-5\n",
    "\n",
    "spec_matrix <- c()\n",
    "spec_matrix <- rbind(spec_matrix, c(0, 0, 0, 1))\n",
    "spec_matrix <- rbind(spec_matrix, c(1, 1, 1, 2))\n",
    "spec_matrix <- rbind(spec_matrix, c(2, 2, 2, 3))\n",
    "\n",
    "lambda_list_etd <- secsse::create_lambda_list(state_names = c(0, 1, 2),\n",
    "                                          num_concealed_states = 3,\n",
    "                                          transition_matrix = spec_matrix,\n",
    "                                          model = \"ETD\")\n",
    "\n",
    "mu_vec_etd <- secsse::create_mu_vector(state_names = c(0, 1, 2),\n",
    "                                   num_concealed_states = 3,\n",
    "                                   model = \"ETD\",\n",
    "                                   lambda_list = lambda_list_etd)\n",
    "\n",
    "mu_vec_etd_sp<-rep(4, length(mu_vec_etd))\n",
    "\n",
    "lambda_list_ctd <- secsse::create_lambda_list(state_names = c(0, 1, 2),\n",
    "                                          num_concealed_states = 3,\n",
    "                                          transition_matrix = spec_matrix,\n",
    "                                          model = \"CTD\")\n",
    "\n",
    "mu_vec_ctd <- secsse::create_mu_vector(state_names = c(0, 1, 2),\n",
    "                                   num_concealed_states = 3,\n",
    "                                   model = \"CTD\",\n",
    "                                   lambda_list = lambda_list_ctd)\n",
    "\n",
    "### Speciation rates\n",
    "\n",
    "mu_vec_ctd_sp<-rep(4, length(mu_vec_ctd))\n",
    "\n",
    "shift_matrix <- c()\n",
    "shift_matrix <- rbind(shift_matrix, c(0, 1, 7))\n",
    "shift_matrix <- rbind(shift_matrix, c(1, 0, 7))\n",
    "shift_matrix <- rbind(shift_matrix, c(1, 2, 8))\n",
    "shift_matrix <- rbind(shift_matrix, c(2, 1, 8))\n",
    "\n",
    "\n",
    "q_matrix <- secsse::create_q_matrix(state_names = c(0, 1, 2),\n",
    "                                    num_concealed_states = 3,\n",
    "                                    shift_matrix = shift_matrix,\n",
    "                                    diff.conceal = TRUE)\n",
    "q_matrix[q_matrix>=9]<-9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extracting diversification and transition rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "vec_AICc<-c()\n",
    "vec_k<-c()\n",
    "vec_L<-c()\n",
    "vec_name<-c()\n",
    "for (file in list.files(\"SecSSE_Results/Size_results/\", full.names = TRUE, pattern = \"size\")){\n",
    "    vec_AICc<-c(vec_AICc, aicc(readRDS(file)$ML,readRDS(file)$k))\n",
    "    vec_k<-c(vec_k, readRDS(file)$k)\n",
    "    vec_L<-c(vec_L, readRDS(file)$ML)\n",
    "    vec_name<-c(vec_name, file)\n",
    "}\n",
    "\n",
    "data_max_Body_size<-cbind(vec_L,vec_k,vec_AICc, vec_name)\n",
    "\n",
    "vec_AICc2<-c()\n",
    "vec_k2<-c()\n",
    "vec_L2<-c()\n",
    "vec_name2<-c()\n",
    "\n",
    "for (i in 1:nrow(data_max_Body_size)){\n",
    "    x<-as.numeric(data_max_Body_size[i,3])\n",
    "    if (i%%3 == 0){\n",
    "        minval<-which.min(c(data_max_Body_size[i-2,3], data_max_Body_size[i-1,3], data_max_Body_size[i,3]))       \n",
    "        minval<-i+minval-3\n",
    "        vec_AICc2<-c(vec_AICc2, vec_AICc[minval])\n",
    "        vec_k2<-c(vec_k2, vec_k[minval])\n",
    "        vec_L2<-c(vec_L2, vec_L[minval])\n",
    "        vec_name2<-c(vec_name2, vec_name[minval])\n",
    "            \n",
    "        if (grepl(\"_model_\",vec_name[minval])){\n",
    "            idparslist <- list()\n",
    "            idparslist[[1]] <- lambda_list_cr\n",
    "            idparslist[[2]] <- mu_vec_cr\n",
    "            idparslist[[3]] <- q_matrix_cr\n",
    "        }\n",
    "\n",
    "        else{\n",
    "            idparslist <- list()\n",
    "            if (grepl(\"_examined_\",vec_name[minval])){\n",
    "                idparslist <- list()\n",
    "                idparslist[[1]] <- lambda_list_etd\n",
    "                idparslist[[2]] <- mu_vec_etd\n",
    "                if (grepl(\"_sp_\",vec_name[minval])){\n",
    "                    idparslist[[2]] <- mu_vec_etd_sp\n",
    "                }\n",
    "                if (grepl(\"_ex_\",vec_name[minval])){\n",
    "                    idparslist[[1]] <- lambda_list_cr\n",
    "                }\n",
    "            }\n",
    "            if (grepl(\"_concealed_\",vec_name[minval])){\n",
    "                idparslist <- list()\n",
    "                idparslist[[1]] <- lambda_list_ctd\n",
    "                idparslist[[2]] <- mu_vec_ctd\n",
    "                if (grepl(\"_sp_\",vec_name[minval])){\n",
    "                    idparslist[[2]] <- mu_vec_ctd_sp\n",
    "                }\n",
    "                if (grepl(\"_ex_\",vec_name[minval])){\n",
    "                    idparslist[[1]] <- lambda_list_cr\n",
    "                }\n",
    "            }\n",
    "\n",
    "            idparslist[[3]] <- q_matrix\n",
    "        }\n",
    "        if((min(vec_AICc2)==vec_AICc[minval])){\n",
    "            print(vec_name[minval])\n",
    "            temp_list<-idparslist\n",
    "        }\n",
    "        print(round(secsse::extract_par_vals(idparslist, readRDS(vec_name[minval])$MLpars), 3))\n",
    "    }\n",
    "}\n",
    "\n",
    "\n",
    "\n",
    "data_max_Body_size<-data.frame(cbind(round(vec_L2,2),vec_k2,round(vec_AICc2,2),round(akaike.weights(vec_AICc2)$deltaAIC, 2), round(akaike.weights(vec_AICc2)$weights,2)))\n",
    "\n",
    "colnames(data_max_Body_size)<-clnames\n",
    "\n",
    "write.table(data_max_Body_size[order(data_max_Body_size[,3]),], \"SecSSE_Results/data_max_Body_size.tsv\", sep =\"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ancestral state estimation on the best-fitting model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "model<-readRDS(\"SecSSE_Results/Size_results/size_concealed_trait_sp_try3.rds\")\n",
    "\n",
    "CTD_par <- secsse::extract_par_vals(temp_list, model$MLpars)\n",
    "\n",
    "CTD_par\n",
    "\n",
    "spec_rates <- CTD_par[1:3]\n",
    "ext_rates <- CTD_par[4:6]\n",
    "Q_Examined <- CTD_par[7:8]\n",
    "Q_Concealed <- CTD_par[9]\n",
    "\n",
    "model$MLpars[[3]]\n",
    "\n",
    "params <- secsse::id_paramPos(c(0, 1, 2), 3)\n",
    "\n",
    "params[[1]][1:3]<-spec_rates[1]\n",
    "params[[1]][4:6]<-spec_rates[2]\n",
    "params[[1]][7:9]<-spec_rates[3]\n",
    "params[[2]][1:3]<-ext_rates[1]\n",
    "params[[2]][4:6]<-ext_rates[1]\n",
    "params[[2]][7:9]<-ext_rates[1]\n",
    "params[[3]]<-model$MLpars[[3]]\n",
    "\n",
    "ll <- secsse::secsse_loglik(parameter = params,\n",
    "                             phy = phy,\n",
    "                             traits = traits,\n",
    "                             num_concealed_states = 3,\n",
    "                             see_ancestral_states = TRUE,\n",
    "                             sampling_fraction = f,\n",
    "                             num_threads =  8)\n",
    "\n",
    "saveRDS(ll$states, \"ASE/bds_SecSSE.rds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reproduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optional step (cleaning data)\n",
    "states<-df$Reproduction\n",
    "names(states)=str_replace((df$Species), \" \", \"_\")\n",
    "states2<-states[!names(states) %in% setdiff(names(states), phy$tip.label)]\n",
    "\n",
    "states2<-states[!names(states) %in% setdiff(names(states), phy$tip.label)]\n",
    "\n",
    "states3<-as.data.frame(cbind(names(states2), as.factor(states2)))\n",
    "colnames(states3)<-c(\"species\",\"states\")\n",
    "rownames(states3)<-NULL\n",
    "traits <- sortingtraits(states3, phy)\n",
    "traits<-traits-1\n",
    "\n",
    "# Computing sampling fractions\n",
    "f<-c(\n",
    "length(na.omit(states2[states2==\"O\"]))/length(na.omit(df$Reproduction[df$Reproduction==\"O\"])),\n",
    "length(na.omit(states2[states2==\"Oo\"]))/length(na.omit(df$Reproduction[df$Reproduction==\"Oo\"])),\n",
    "length(na.omit(states2[states2==\"PV\"]))/length(na.omit(df$Reproduction[df$Reproduction==\"PV\"])),\n",
    "length(na.omit(states2[states2==\"YV\"]))/length(na.omit(df$Reproduction[df$Reproduction==\"YV\"])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementing model structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_matrix_cr <- c()\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(0, 0, 0, 1))\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(1, 1, 1, 1))\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(2, 2, 2, 1))\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(3, 3, 3, 1))\n",
    "\n",
    "lambda_list_cr <- secsse::create_lambda_list(state_names = c(0, 1, 2, 3),\n",
    "                                          num_concealed_states = 4,\n",
    "                                          transition_matrix = spec_matrix_cr,\n",
    "                                          model = \"CR\")\n",
    "\n",
    "mu_vec_cr <- secsse::create_mu_vector(state_names = c(0, 1, 2, 3),\n",
    "                                   num_concealed_states = 4,\n",
    "                                   model = \"CR\",\n",
    "                                   lambda_list = lambda_list_cr)\n",
    "\n",
    "shift_matrix_cr <- c()\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(0, 3, 3))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(3, 0, 3))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(1, 3, 3))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(3, 1, 3))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(2, 3, 3))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(3, 2, 3))\n",
    "\n",
    "q_matrix_cr <- secsse::create_q_matrix(state_names = c(0, 1, 2, 3),\n",
    "                                    num_concealed_states = 4,\n",
    "                                    shift_matrix = shift_matrix_cr,\n",
    "                                    diff.conceal = TRUE)\n",
    "q_matrix_cr[q_matrix_cr>=4]<-4\n",
    "\n",
    "spec_matrix <- c()\n",
    "spec_matrix <- rbind(spec_matrix, c(0, 0, 0, 1))\n",
    "spec_matrix <- rbind(spec_matrix, c(1, 1, 1, 2))\n",
    "spec_matrix <- rbind(spec_matrix, c(2, 2, 2, 3))\n",
    "spec_matrix <- rbind(spec_matrix, c(3, 3, 3, 4))\n",
    "\n",
    "lambda_list_etd <- secsse::create_lambda_list(state_names = c(0, 1, 2, 3),\n",
    "                                          num_concealed_states = 4,\n",
    "                                          transition_matrix = spec_matrix,\n",
    "                                          model = \"ETD\")\n",
    "\n",
    "mu_vec_etd <- secsse::create_mu_vector(state_names = c(0, 1, 2, 3),\n",
    "                                   num_concealed_states = 4,\n",
    "                                   model = \"ETD\",\n",
    "                                   lambda_list = lambda_list_etd)\n",
    "\n",
    "shift_matrix <- c()\n",
    "shift_matrix <- rbind(shift_matrix, c(0, 3, 9))\n",
    "shift_matrix <- rbind(shift_matrix, c(3, 0, 9))\n",
    "shift_matrix <- rbind(shift_matrix, c(1, 3, 9))\n",
    "shift_matrix <- rbind(shift_matrix, c(3, 1, 9))\n",
    "shift_matrix <- rbind(shift_matrix, c(2, 3, 9))\n",
    "shift_matrix <- rbind(shift_matrix, c(3, 2, 9))\n",
    "\n",
    "q_matrix <- secsse::create_q_matrix(state_names = c(0, 1, 2, 3),\n",
    "                                    num_concealed_states = 4,\n",
    "                                    shift_matrix = shift_matrix,\n",
    "                                    diff.conceal = TRUE)\n",
    "q_matrix[q_matrix>=10]<-10\n",
    "\n",
    "mu_vec_etd_sp<-rep(5, length(mu_vec_etd))\n",
    "\n",
    "lambda_list_ctd <- secsse::create_lambda_list(state_names = c(0, 1, 2, 3),\n",
    "                                          num_concealed_states = 4,\n",
    "                                          transition_matrix = spec_matrix,\n",
    "                                          model = \"CTD\")\n",
    "\n",
    "mu_vec_ctd <- secsse::create_mu_vector(state_names = c(0, 1, 2, 3),\n",
    "                                   num_concealed_states = 4,\n",
    "                                   model = \"CTD\",\n",
    "                                   lambda_list = lambda_list_ctd)\n",
    "\n",
    "mu_vec_ctd_sp<-rep(5, length(mu_vec_ctd))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extracting diversification and transition rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vec_AICc<-c()\n",
    "vec_k<-c()\n",
    "vec_L<-c()\n",
    "vec_name<-c()\n",
    "for (file in list.files(\"SecSSE_Results/Reproduction_results/\", full.names = TRUE, pattern = \"reproduction\")){\n",
    "    vec_AICc<-c(vec_AICc, aicc(readRDS(file)$ML,readRDS(file)$k))\n",
    "    vec_k<-c(vec_k, readRDS(file)$k)\n",
    "    vec_L<-c(vec_L, readRDS(file)$ML)\n",
    "    vec_name<-c(vec_name, file)\n",
    "}\n",
    "\n",
    "data_max_Reproduction<-cbind(vec_L,vec_k,vec_AICc, vec_name)\n",
    "\n",
    "vec_AICc2<-c()\n",
    "vec_k2<-c()\n",
    "vec_L2<-c()\n",
    "vec_name2<-c()\n",
    "\n",
    "for (i in 1:nrow(data_max_Reproduction)){\n",
    "    x<-as.numeric(data_max_Reproduction[i,3])\n",
    "    if (i%%3 == 0){\n",
    "        minval<-which.min(c(data_max_Reproduction[i-2,3], data_max_Reproduction[i-1,3], data_max_Reproduction[i,3]))       \n",
    "        minval<-i+minval-3\n",
    "        vec_AICc2<-c(vec_AICc2, vec_AICc[minval])\n",
    "        vec_k2<-c(vec_k2, vec_k[minval])\n",
    "        vec_L2<-c(vec_L2, vec_L[minval])\n",
    "        vec_name2<-c(vec_name2, vec_name[minval])\n",
    "            \n",
    "        if (grepl(\"_model_\",vec_name[minval])){\n",
    "            idparslist <- list()\n",
    "            idparslist[[1]] <- lambda_list_cr\n",
    "            idparslist[[2]] <- mu_vec_cr\n",
    "            idparslist[[3]] <- q_matrix_cr\n",
    "        }\n",
    "\n",
    "        else{\n",
    "            idparslist <- list()\n",
    "            if (grepl(\"_examined_\",vec_name[minval])){\n",
    "                idparslist <- list()\n",
    "                idparslist[[1]] <- lambda_list_etd\n",
    "                idparslist[[2]] <- mu_vec_etd\n",
    "                if (grepl(\"_sp_\",vec_name[minval])){\n",
    "                    idparslist[[2]] <- mu_vec_etd_sp\n",
    "                }\n",
    "                if (grepl(\"_ex_\",vec_name[minval])){\n",
    "                    idparslist[[1]] <- lambda_list_cr\n",
    "                }\n",
    "            }\n",
    "            if (grepl(\"_concealed_\",vec_name[minval])){\n",
    "                idparslist <- list()\n",
    "                idparslist[[1]] <- lambda_list_ctd\n",
    "                idparslist[[2]] <- mu_vec_ctd\n",
    "                if (grepl(\"_sp_\",vec_name[minval])){\n",
    "                    idparslist[[2]] <- mu_vec_ctd_sp\n",
    "                }\n",
    "                if (grepl(\"_ex_\",vec_name[minval])){\n",
    "                    idparslist[[1]] <- lambda_list_cr\n",
    "                }\n",
    "            }\n",
    "\n",
    "            idparslist[[3]] <- q_matrix\n",
    "        }\n",
    "        if((min(vec_AICc2)==vec_AICc[minval])){\n",
    "            print(vec_name[minval])\n",
    "            temp_list<-idparslist\n",
    "        }\n",
    "        print(round(secsse::extract_par_vals(idparslist, readRDS(vec_name[minval])$MLpars),3))\n",
    "    }\n",
    "}\n",
    "\n",
    "\n",
    "\n",
    "data_max_Reproduction<-data.frame(cbind(round(vec_L2,2),vec_k2,round(vec_AICc2,2),round(akaike.weights(vec_AICc2)$deltaAIC, 2), round(akaike.weights(vec_AICc2)$weights,2)))\n",
    "\n",
    "colnames(data_max_Reproduction)<-clnames\n",
    "\n",
    "write.table(data_max_Reproduction[order(data_max_Reproduction[,3]),], \"SecSSE_Results/data_max_Reproduction.tsv\", sep =\"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ancestral state estimation on the best-fitting model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model<-readRDS(\"SecSSE_Results/Reproduction_results/reproduction_concealed_trait_sp_try1.rds\")\n",
    "\n",
    "CTD_par <- secsse::extract_par_vals(temp_list, model$MLpars)\n",
    "\n",
    "CTD_par\n",
    "\n",
    "spec_rates <- CTD_par[1:4]\n",
    "ext_rates <- CTD_par[5:8]\n",
    "\n",
    "model$MLpars[[3]]\n",
    "\n",
    "params <- secsse::id_paramPos(c(0, 1, 2, 3), 4)\n",
    "\n",
    "params[[1]][1:4]<-spec_rates[1]\n",
    "params[[1]][5:8]<-spec_rates[2]\n",
    "params[[1]][9:12]<-spec_rates[3]\n",
    "params[[1]][13:16]<-spec_rates[4]\n",
    "params[[2]][1:4]<-ext_rates[1]\n",
    "params[[2]][5:8]<-ext_rates[1]\n",
    "params[[2]][9:12]<-ext_rates[1]\n",
    "params[[2]][13:16]<-ext_rates[1]\n",
    "params[[3]]<-model$MLpars[[3]]\n",
    "\n",
    "ll <- secsse::secsse_loglik(parameter = params,\n",
    "                             phy = phy,\n",
    "                             traits = traits,\n",
    "                             num_concealed_states = 4,\n",
    "                             see_ancestral_states = TRUE,\n",
    "                             sampling_fraction = f,\n",
    "                             num_threads =  8)\n",
    "ll\n",
    "\n",
    "saveRDS(ll$states, \"ASE/rp_SecSSE.rds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Habitat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optional step (cleaning data)\n",
    "states<-df$Habitat\n",
    "names(states)=str_replace((df$Species), \" \", \"_\")\n",
    "states2<-states[!names(states) %in% setdiff(names(states), phy$tip.label)]\n",
    "setdiff(phy$tip.label,names(states2))\n",
    "states3<-as.data.frame(cbind(names(states2), as.factor(states2)))\n",
    "colnames(states3)<-c(\"species\",\"states\")\n",
    "rownames(states3)<-NULL\n",
    "traits <- sortingtraits(states3, phy)\n",
    "traits<-traits-1\n",
    "# Computing sampling fractions\n",
    "f<-c(\n",
    "length(na.omit(states2[states2==\"D\"]))/length(na.omit(df$Habitat[df$Habitat==\"D\"])),\n",
    "length(na.omit(states2[states2==\"IS\"]))/length(na.omit(df$Habitat[df$Habitat==\"IS\"])),\n",
    "length(na.omit(states2[states2==\"O\"]))/length(na.omit(df$Habitat[df$Habitat==\"O\"])),\n",
    "length(na.omit(states2[states2==\"OS\"]))/length(na.omit(df$Habitat[df$Habitat==\"OS\"])),\n",
    "length(na.omit(states2[states2==\"R\"]))/length(na.omit(df$Habitat[df$Habitat==\"R\"])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementing model structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_matrix_cr <- c()\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(0, 0, 0, 1))\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(1, 1, 1, 1))\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(2, 2, 2, 1))\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(3, 3, 3, 1))\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(4, 4, 4, 1))\n",
    "lambda_list_cr <- secsse::create_lambda_list(state_names = c(0, 1, 2, 3, 4),\n",
    "                                          num_concealed_states = 5,\n",
    "                                          transition_matrix = spec_matrix_cr,\n",
    "                                          model = \"CR\")\n",
    "mu_vec_cr <- secsse::create_mu_vector(state_names = c(0, 1, 2, 3, 4),\n",
    "                                   num_concealed_states = 5,\n",
    "                                   model = \"CR\",\n",
    "                                   lambda_list = lambda_list_cr)\n",
    "shift_matrix_cr <- c()\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(0, 3, 3))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(1, 2, 4))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(1, 3, 5))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(1, 4, 6))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(2, 3, 7))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(2, 4, 8))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(3, 0, 9))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(3, 1, 10))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(4, 1, 11))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(4, 2, 12))\n",
    "shift_matrix_cr <- rbind(shift_matrix_cr, c(4, 3, 13))\n",
    "q_matrix_cr <- secsse::create_q_matrix(state_names = c(0, 1, 2, 3, 4),\n",
    "                                    num_concealed_states = 5,\n",
    "                                    shift_matrix = shift_matrix_cr,\n",
    "                                    diff.conceal = TRUE)\n",
    "q_matrix_cr[q_matrix_cr>=14]<-14\n",
    "spec_matrix <- c()\n",
    "spec_matrix <- rbind(spec_matrix, c(0, 0, 0, 1))\n",
    "spec_matrix <- rbind(spec_matrix, c(1, 1, 1, 2))\n",
    "spec_matrix <- rbind(spec_matrix, c(2, 2, 2, 3))\n",
    "spec_matrix <- rbind(spec_matrix, c(3, 3, 3, 4))\n",
    "spec_matrix <- rbind(spec_matrix, c(4, 4, 4, 5))\n",
    "lambda_list_etd <- secsse::create_lambda_list(state_names = c(0, 1, 2, 3, 4),\n",
    "                                          num_concealed_states = 5,\n",
    "                                          transition_matrix = spec_matrix,\n",
    "                                          model = \"ETD\")\n",
    "mu_vec_etd <- secsse::create_mu_vector(state_names = c(0, 1, 2, 3, 4),\n",
    "                                   num_concealed_states = 5,\n",
    "                                   model = \"ETD\",\n",
    "                                   lambda_list = lambda_list_etd)\n",
    "shift_matrix <- c()\n",
    "shift_matrix <- c()\n",
    "shift_matrix <- rbind(shift_matrix, c(0, 3, 11))\n",
    "shift_matrix <- rbind(shift_matrix, c(1, 2, 12))\n",
    "shift_matrix <- rbind(shift_matrix, c(1, 3, 13))\n",
    "shift_matrix <- rbind(shift_matrix, c(1, 4, 14))\n",
    "shift_matrix <- rbind(shift_matrix, c(2, 3, 15))\n",
    "shift_matrix <- rbind(shift_matrix, c(2, 4, 16))\n",
    "shift_matrix <- rbind(shift_matrix, c(3, 0, 17))\n",
    "shift_matrix <- rbind(shift_matrix, c(3, 1, 18))\n",
    "shift_matrix <- rbind(shift_matrix, c(4, 1, 19))\n",
    "shift_matrix <- rbind(shift_matrix, c(4, 2, 20))\n",
    "shift_matrix <- rbind(shift_matrix, c(4, 3, 21))\n",
    "q_matrix <- secsse::create_q_matrix(state_names = c(0, 1, 2, 3, 4),\n",
    "                                    num_concealed_states = 5,\n",
    "                                    shift_matrix = shift_matrix,\n",
    "                                    diff.conceal = TRUE)\n",
    "q_matrix[q_matrix>=22]<-22\n",
    "mu_vec_etd_sp<-rep(6, length(mu_vec_etd))\n",
    "lambda_list_ctd <- secsse::create_lambda_list(state_names = c(0, 1, 2, 3, 4),\n",
    "                                          num_concealed_states = 5,\n",
    "                                          transition_matrix = spec_matrix,\n",
    "                                          model = \"CTD\")\n",
    "mu_vec_ctd <- secsse::create_mu_vector(state_names = c(0, 1, 2, 3, 4),\n",
    "                                   num_concealed_states = 5,\n",
    "                                   model = \"CTD\",\n",
    "                                   lambda_list = lambda_list_ctd)\n",
    "mu_vec_ctd_sp<-rep(6, length(mu_vec_ctd))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extracting diversification and transition rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "vec_AICc<-c()\n",
    "vec_k<-c()\n",
    "vec_L<-c()\n",
    "vec_name<-c()\n",
    "for (file in list.files(\"SecSSE_Results/Habitat_results/\", full.names = TRUE, pattern = \"habitat\")){\n",
    "    vec_AICc<-c(vec_AICc, aicc(readRDS(file)$ML,readRDS(file)$k))\n",
    "    vec_k<-c(vec_k, readRDS(file)$k)\n",
    "    vec_L<-c(vec_L, readRDS(file)$ML)\n",
    "    vec_name<-c(vec_name, file)\n",
    "}\n",
    "\n",
    "data_max_Habitat<-cbind(vec_L,vec_k,vec_AICc, vec_name)\n",
    "\n",
    "vec_AICc2<-c()\n",
    "vec_k2<-c()\n",
    "vec_L2<-c()\n",
    "vec_name2<-c()\n",
    "\n",
    "for (i in 1:nrow(data_max_Habitat)){\n",
    "    x<-as.numeric(data_max_Habitat[i,3])\n",
    "    if (i%%3 == 0){\n",
    "        minval<-which.min(c(data_max_Habitat[i-2,3], data_max_Habitat[i-1,3], data_max_Habitat[i,3]))       \n",
    "        minval<-i+minval-3\n",
    "        vec_AICc2<-c(vec_AICc2, vec_AICc[minval])\n",
    "        vec_k2<-c(vec_k2, vec_k[minval])\n",
    "        vec_L2<-c(vec_L2, vec_L[minval])\n",
    "        vec_name2<-c(vec_name2, vec_name[minval])\n",
    "            \n",
    "        if (grepl(\"_model_\",vec_name[minval])){\n",
    "            idparslist <- list()\n",
    "            idparslist[[1]] <- lambda_list_cr\n",
    "            idparslist[[2]] <- mu_vec_cr\n",
    "            idparslist[[3]] <- q_matrix_cr\n",
    "        }\n",
    "\n",
    "        else{\n",
    "            idparslist <- list()\n",
    "            if (grepl(\"_examined_\",vec_name[minval])){\n",
    "                idparslist <- list()\n",
    "                idparslist[[1]] <- lambda_list_etd\n",
    "                idparslist[[2]] <- mu_vec_etd\n",
    "                if (grepl(\"_sp_\",vec_name[minval])){\n",
    "                    idparslist[[2]] <- mu_vec_etd_sp\n",
    "                }\n",
    "                if (grepl(\"_ex_\",vec_name[minval])){\n",
    "                    idparslist[[1]] <- lambda_list_cr\n",
    "                }\n",
    "            }\n",
    "            if (grepl(\"_concealed_\",vec_name[minval])){\n",
    "                idparslist <- list()\n",
    "                idparslist[[1]] <- lambda_list_ctd\n",
    "                idparslist[[2]] <- mu_vec_ctd\n",
    "                if (grepl(\"_sp_\",vec_name[minval])){\n",
    "                    idparslist[[2]] <- mu_vec_ctd_sp\n",
    "                }\n",
    "                if (grepl(\"_ex_\",vec_name[minval])){\n",
    "                    idparslist[[1]] <- lambda_list_cr\n",
    "                }\n",
    "            }\n",
    "\n",
    "            idparslist[[3]] <- q_matrix\n",
    "        }\n",
    "        if((min(vec_AICc2)==vec_AICc[minval])){\n",
    "            print(vec_name[minval])\n",
    "            temp_list<-idparslist\n",
    "        }\n",
    "        print(round(secsse::extract_par_vals(idparslist, readRDS(vec_name[minval])$MLpars), 3))\n",
    "    }\n",
    "}\n",
    "\n",
    "\n",
    "\n",
    "data_max_Habitat<-data.frame(cbind(round(vec_L2,2),vec_k2,round(vec_AICc2,2),round(akaike.weights(vec_AICc2)$deltaAIC, 2), round(akaike.weights(vec_AICc2)$weights,2)))\n",
    "\n",
    "colnames(data_max_Habitat)<-clnames\n",
    "\n",
    "write.table(data_max_Habitat[order(data_max_Habitat[,3]),], \"SecSSE_Results/data_max_Habitat.tsv\", sep =\"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ancestral state estimation on the best-fitting model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "model<-readRDS(\"SecSSE_Results/Habitat_results/habitat_concealed_trait_sp_try3.rds\")\n",
    "\n",
    "CTD_par <- secsse::extract_par_vals(temp_list, model$MLpars)\n",
    "\n",
    "CTD_par\n",
    "\n",
    "spec_rates <- CTD_par[1:5]\n",
    "ext_rates <- CTD_par[6:10]\n",
    "\n",
    "model$MLpars[[3]]\n",
    "\n",
    "params <- secsse::id_paramPos(c(0, 1, 2, 3, 4), 5)\n",
    "\n",
    "params[[1]][1:5]<-spec_rates[1]\n",
    "params[[1]][6:10]<-spec_rates[2]\n",
    "params[[1]][11:15]<-spec_rates[3]\n",
    "params[[1]][16:20]<-spec_rates[4]\n",
    "params[[1]][21:25]<-spec_rates[5]\n",
    "params[[2]][1:5]<-ext_rates[1]\n",
    "params[[2]][6:10]<-ext_rates[1]\n",
    "params[[2]][11:16]<-ext_rates[1]\n",
    "params[[2]][16:20]<-ext_rates[1]\n",
    "params[[1]][21:25]<-spec_rates[1]\n",
    "params[[3]]<-model$MLpars[[3]]\n",
    "\n",
    "ll <- secsse::secsse_loglik(parameter = params,\n",
    "                             phy = phy,\n",
    "                             traits = traits,\n",
    "                             num_concealed_states = 5,\n",
    "                             see_ancestral_states = TRUE,\n",
    "                             sampling_fraction = f,\n",
    "                             num_threads =  8)\n",
    "\n",
    "saveRDS(ll$states, \"ASE/ht_SecSSE.rds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df<-read.csv(\"table_diet.tsv\", sep=\"\\t\")\n",
    "\n",
    "df$Species<-str_replace((df$Species), \" \", \"_\")\n",
    "states<-cbind(c(str_replace((df$Species), \" \", \"_\"), setdiff(phy$tip.label, str_replace((df$Species), \" \", \"_\"))), c(df$Diet, rep(NA, 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",
    "colnames(states_traits)<-c(\"Species\", \"Diet\")\n",
    "\n",
    "states<- states_traits$Diet\n",
    "names(states)= states_traits$Species\n",
    "\n",
    "states2<-states\n",
    "\n",
    "# Computing sampling fractions\n",
    "f<-c(\n",
    "(length(states2[states2==\"Invertebrate-feeders\"])-length(states2[is.na(states2)]))/length(df$Diet[df$Diet==\"Invertebrate-feeders\"]),\n",
    "(length(states2[states2==\"Mesopredator\"])-length(states2[is.na(states2)]))/length(df$Diet[df$Diet==\"Mesopredator\"]),\n",
    "(length(states2[states2==\"Macropredator\"])-length(states2[is.na(states2)]))/length(df$Diet[df$Diet==\"Macropredator\"])\n",
    ")\n",
    "\n",
    "\n",
    "states3<-as.data.frame(cbind(names(states2), as.factor(states2)))\n",
    "colnames(states3)<-c(\"species\",\"states\")\n",
    "rownames(states3)<-NULL\n",
    "traits <- sortingtraits(states3, phy)\n",
    "traits<-traits-1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementing model structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_matrix_cr <- c()\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(0, 0, 0, 1))\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(1, 1, 1, 1))\n",
    "spec_matrix_cr <- rbind(spec_matrix_cr, c(2, 2, 2, 1))\n",
    "\n",
    "lambda_list_cr <- secsse::create_lambda_list(state_names = c(0, 1, 2),\n",
    "                                          num_concealed_states = 3,\n",
    "                                          transition_matrix = spec_matrix_cr,\n",
    "                                          model = \"CR\")\n",
    "mu_vec_cr <- secsse::create_mu_vector(state_names = c(0, 1, 2),\n",
    "                                   num_concealed_states = 3,\n",
    "                                   model = \"CR\",\n",
    "                                   lambda_list = lambda_list_cr)\n",
    "\n",
    "shift_matrix_cr  <- c()\n",
    "shift_matrix_cr  <- rbind(shift_matrix_cr , c(0, 2, 3))\n",
    "shift_matrix_cr  <- rbind(shift_matrix_cr , c(2, 0, 4))\n",
    "shift_matrix_cr  <- rbind(shift_matrix_cr , c(1, 2, 5))\n",
    "shift_matrix_cr  <- rbind(shift_matrix_cr , c(2, 1, 6))\n",
    "\n",
    "q_matrix_cr <- secsse::create_q_matrix(state_names = c(0, 1, 2),\n",
    "                                    num_concealed_states = 3,\n",
    "                                    shift_matrix = shift_matrix_cr,\n",
    "                                    diff.conceal = TRUE)\n",
    "q_matrix_cr[q_matrix_cr>=7]<-7\n",
    "\n",
    "spec_matrix <- c()\n",
    "spec_matrix <- rbind(spec_matrix, c(0, 0, 0, 1))\n",
    "spec_matrix <- rbind(spec_matrix, c(1, 1, 1, 2))\n",
    "spec_matrix <- rbind(spec_matrix, c(2, 2, 2, 3))\n",
    "\n",
    "lambda_list_etd <- secsse::create_lambda_list(state_names = c(0, 1, 2),\n",
    "                                          num_concealed_states = 3,\n",
    "                                          transition_matrix = spec_matrix,\n",
    "                                          model = \"ETD\")\n",
    "\n",
    "mu_vec_etd <- secsse::create_mu_vector(state_names = c(0, 1, 2),\n",
    "                                   num_concealed_states = 3,\n",
    "                                   model = \"ETD\",\n",
    "                                   lambda_list = lambda_list_etd)\n",
    "\n",
    "mu_vec_etd_sp<-rep(4, length(mu_vec_etd))\n",
    "\n",
    "lambda_list_ctd <- secsse::create_lambda_list(state_names = c(0, 1, 2),\n",
    "                                          num_concealed_states = 3,\n",
    "                                          transition_matrix = spec_matrix,\n",
    "                                          model = \"CTD\")\n",
    "\n",
    "mu_vec_ctd <- secsse::create_mu_vector(state_names = c(0, 1, 2),\n",
    "                                   num_concealed_states = 3,\n",
    "                                   model = \"CTD\",\n",
    "                                   lambda_list = lambda_list_ctd)\n",
    "\n",
    "### Speciation rates\n",
    "\n",
    "mu_vec_ctd_sp<-rep(4, length(mu_vec_ctd))\n",
    "\n",
    "shift_matrix <- c()\n",
    "shift_matrix <- rbind(shift_matrix, c(0, 2, 7))\n",
    "shift_matrix <- rbind(shift_matrix, c(2, 0, 8))\n",
    "shift_matrix <- rbind(shift_matrix, c(1, 2, 9))\n",
    "shift_matrix <- rbind(shift_matrix, c(2, 1, 10))\n",
    "\n",
    "\n",
    "q_matrix <- secsse::create_q_matrix(state_names = c(0, 1, 2),\n",
    "                                    num_concealed_states = 3,\n",
    "                                    shift_matrix = shift_matrix,\n",
    "                                    diff.conceal = TRUE)\n",
    "q_matrix[q_matrix>=11]<-11"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extracting diversification and transition rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "vec_AICc<-c()\n",
    "vec_k<-c()\n",
    "vec_L<-c()\n",
    "vec_name<-c()\n",
    "for (file in list.files(\"SecSSE_Results/Diet_results/\", full.names = TRUE, pattern = \"diet\")){\n",
    "    vec_AICc<-c(vec_AICc, aicc(readRDS(file)$ML,readRDS(file)$k))\n",
    "    vec_k<-c(vec_k, readRDS(file)$k)\n",
    "    vec_L<-c(vec_L, readRDS(file)$ML)\n",
    "    vec_name<-c(vec_name, file)\n",
    "}\n",
    "\n",
    "data_Diet<-cbind(vec_L,vec_k,vec_AICc, vec_name)\n",
    "\n",
    "vec_AICc2<-c()\n",
    "vec_k2<-c()\n",
    "vec_L2<-c()\n",
    "vec_name2<-c()\n",
    "\n",
    "for (i in 1:nrow(data_Diet)){\n",
    "    x<-as.numeric(data_Diet[i,3])\n",
    "    if (i%%3 == 0){\n",
    "        minval<-which.min(c(data_Diet[i-2,3], data_Diet[i-1,3], data_Diet[i,3]))       \n",
    "        minval<-i+minval-3\n",
    "        vec_AICc2<-c(vec_AICc2, vec_AICc[minval])\n",
    "        vec_k2<-c(vec_k2, vec_k[minval])\n",
    "        vec_L2<-c(vec_L2, vec_L[minval])\n",
    "        vec_name2<-c(vec_name2, vec_name[minval])\n",
    "            \n",
    "        if (grepl(\"_model_\",vec_name[minval])){\n",
    "            idparslist <- list()\n",
    "            idparslist[[1]] <- lambda_list_cr\n",
    "            idparslist[[2]] <- mu_vec_cr\n",
    "            idparslist[[3]] <- q_matrix_cr\n",
    "        }\n",
    "\n",
    "        else{\n",
    "            idparslist <- list()\n",
    "            if (grepl(\"_examined_\",vec_name[minval])){\n",
    "                idparslist <- list()\n",
    "                idparslist[[1]] <- lambda_list_etd\n",
    "                idparslist[[2]] <- mu_vec_etd\n",
    "                if (grepl(\"_sp_\",vec_name[minval])){\n",
    "                    idparslist[[2]] <- mu_vec_etd_sp\n",
    "                }\n",
    "                if (grepl(\"_ex_\",vec_name[minval])){\n",
    "                    idparslist[[1]] <- lambda_list_cr\n",
    "                }\n",
    "            }\n",
    "            if (grepl(\"_concealed_\",vec_name[minval])){\n",
    "                idparslist <- list()\n",
    "                idparslist[[1]] <- lambda_list_ctd\n",
    "                idparslist[[2]] <- mu_vec_ctd\n",
    "                if (grepl(\"_sp_\",vec_name[minval])){\n",
    "                    idparslist[[2]] <- mu_vec_ctd_sp\n",
    "                }\n",
    "                if (grepl(\"_ex_\",vec_name[minval])){\n",
    "                    idparslist[[1]] <- lambda_list_cr\n",
    "                }\n",
    "            }\n",
    "\n",
    "            idparslist[[3]] <- q_matrix\n",
    "        }\n",
    "        if((min(vec_AICc2)==vec_AICc[minval])){\n",
    "            print(vec_name[minval])\n",
    "            temp_list<-idparslist\n",
    "        }\n",
    "        print(round(secsse::extract_par_vals(idparslist, readRDS(vec_name[minval])$MLpars), 3))\n",
    "    }\n",
    "}\n",
    "\n",
    "\n",
    "\n",
    "data_Diet<-data.frame(cbind(round(vec_L2,2),vec_k2,round(vec_AICc2,2),round(akaike.weights(vec_AICc2)$deltaAIC, 2), round(akaike.weights(vec_AICc2)$weights,2)))\n",
    "\n",
    "colnames(data_Diet)<-clnames\n",
    "\n",
    "write.table(data_Diet[order(data_Diet[,3]),], \"SecSSE_Results/data_Diet.tsv\", sep =\"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ancestral state estimation on the best-fitting model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model<-readRDS(\"SecSSE_Results/Diet_results/diet_concealed_trait_sp_try1.rds\")\n",
    "\n",
    "CTD_par <- secsse::extract_par_vals(temp_list, model$MLpars)\n",
    "\n",
    "CTD_par\n",
    "\n",
    "spec_rates <- CTD_par[1:3]\n",
    "ext_rates <- CTD_par[4:6]\n",
    "Q_Examined <- CTD_par[7:10]\n",
    "Q_Concealed <- CTD_par[11]\n",
    "\n",
    "model$MLpars[[3]]\n",
    "\n",
    "params <- secsse::id_paramPos(c(0, 1, 2), 3)\n",
    "\n",
    "params[[1]][1:3]<-spec_rates[1]\n",
    "params[[1]][4:6]<-spec_rates[2]\n",
    "params[[1]][7:9]<-spec_rates[3]\n",
    "params[[2]][1:3]<-ext_rates[1]\n",
    "params[[2]][4:6]<-ext_rates[1]\n",
    "params[[2]][7:9]<-ext_rates[1]\n",
    "params[[3]]<-model$MLpars[[3]]\n",
    "\n",
    "ll <- secsse::secsse_loglik(parameter = params,\n",
    "                             phy = phy,\n",
    "                             traits = traits,\n",
    "                             num_concealed_states = 3,\n",
    "                             see_ancestral_states = TRUE,\n",
    "                             sampling_fraction = f,\n",
    "                             num_threads =  8)\n",
    "\n",
    "saveRDS(ll$states, \"ASE_diet.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
}