SecSSE_Diet.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SecSSE Diet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Loading packages\n",
    "library(\"Rtsne\")\n",
    "library(\"ape\")\n",
    "library(\"secsse\")\n",
    "library(\"DDD\")\n",
    "library(\"doMC\")\n",
    "library(\"tidyverse\")\n",
    "library(\"qgraph\")\n",
    "library(\"stringr\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading data\n",
    "\n",
    "phy<-read.nexus(\"16FC_16C_374_sp.tree\")\n",
    "df<-read.csv(\"table_diet.tsv\", sep=\"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optional step (cleaning data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Optional step (cleaning data)\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",
    "states_traits<-as.data.frame(states[!states[,1] %in% setdiff(states[,1], phy$tip.label),])\n",
    "colnames(states_traits)<-c(\"Species\", \"Diet\")\n",
    "states<- states_traits$Diet\n",
    "names(states)= states_traits$Species\n",
    "\n",
    "states2<-states"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing sampling fractions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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",
    "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "write.table(cbind(phy$tip.label,traits), \"table_data_diet.tsv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Starting values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BD_model <- bd_ML(brts = ape::branching.times(phy))\n",
    "BD_lambda <- BD_model$lambda0\n",
    "BD_mu <- BD_model$mu0\n",
    "\n",
    "strt_lambda<-c(BD_lambda, BD_lambda/2, BD_lambda*2)\n",
    "strt_mu <- c(BD_mu, BD_mu/2, BD_mu*2)\n",
    "strt_q <- strt_lambda/3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constant Rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_matrix <- c()\n",
    "spec_matrix <- rbind(spec_matrix, c(0, 0, 0, 1))\n",
    "spec_matrix <- rbind(spec_matrix, c(1, 1, 1, 1))\n",
    "spec_matrix <- rbind(spec_matrix, c(2, 2, 2, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "lambda_list_cr <- secsse::create_lambda_list(state_names = c(0, 1, 2),\n",
    "                                          num_concealed_states = 3,\n",
    "                                          transition_matrix = spec_matrix,\n",
    "                                          model = \"CR\")\n",
    "lambda_list_cr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "mu_vec_cr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "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",
    "q_matrix_cr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### Try 1\n",
    "\n",
    "idparsopt <- c(1:7) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(strt_lambda[1], strt_mu[1], rep(strt_q[1],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_cr\n",
    "idparslist[[2]] <- mu_vec_cr\n",
    "idparslist[[3]] <- q_matrix_cr\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 7\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_model_try1.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "\n",
    "idparsopt <- c(1:7) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(strt_lambda[2], strt_mu[2], rep(strt_q[2],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_cr\n",
    "idparslist[[2]] <- mu_vec_cr\n",
    "idparslist[[3]] <- q_matrix_cr\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 7\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_model_try2.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "\n",
    "idparsopt <- c(1:7) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(strt_lambda[3], strt_mu[3], rep(strt_q[3],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_cr\n",
    "idparslist[[2]] <- mu_vec_cr\n",
    "idparslist[[3]] <- q_matrix_cr\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 7\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_model_try3.rds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examined trait "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "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",
    "lambda_list_etd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "mu_vec_etd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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\n",
    "q_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Speciation rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_vec_etd_sp<-rep(4, length(mu_vec_etd))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "idparsopt <- c(1:4, 7:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[1],3), rep(strt_mu[1],1), rep(strt_q[1],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_etd\n",
    "idparslist[[2]] <- mu_vec_etd_sp\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 9\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_examined_trait_sp_try1.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1:4, 7:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[2],3), rep(strt_mu[2],1), rep(strt_q[2],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_etd\n",
    "idparslist[[2]] <- mu_vec_etd_sp\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 9\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_examined_trait_sp_try2.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1:4, 7:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[3],3), rep(strt_mu[3],1), rep(strt_q[3],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_etd\n",
    "idparslist[[2]] <- mu_vec_etd_sp\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 9\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_examined_trait_sp_try3.rds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extinction rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "idparsopt <- c(1, 4:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[1],1), rep(strt_mu[1],3), rep(strt_q[1],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_cr\n",
    "idparslist[[2]] <- mu_vec_etd\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 9\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_examined_trait_mu_try1.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1, 4:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[2],1), rep(strt_mu[2],3), rep(strt_q[2],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_cr\n",
    "idparslist[[2]] <- mu_vec_etd\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 9\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_examined_trait_mu_try2.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1, 4:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[3],1), rep(strt_mu[3],3), rep(strt_q[3],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_cr\n",
    "idparslist[[2]] <- mu_vec_etd\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 9\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_examined_trait_mu_try3.rds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Diversification rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "idparsopt <- c(1:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[1],3), rep(strt_mu[1],3), rep(strt_q[1],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_etd\n",
    "idparslist[[2]] <- mu_vec_etd\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 11\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_examined_trait_net_diversification_try1.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[2],3), rep(strt_mu[2],3), rep(strt_q[2],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_etd\n",
    "idparslist[[2]] <- mu_vec_etd\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 11\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_examined_trait_net_diversification_try2.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[3],3), rep(strt_mu[3],3), rep(strt_q[3],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_etd\n",
    "idparslist[[2]] <- mu_vec_etd\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 11\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_examined_trait_net_diversification_try3.rds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conceal trait"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "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",
    "lambda_list_ctd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "mu_vec_ctd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Speciation rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_vec_ctd_sp<-rep(4, length(mu_vec_ctd))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "idparsopt <- c(1:4, 7:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[1],3), rep(strt_mu[1],1), rep(strt_q[1],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_ctd\n",
    "idparslist[[2]] <- mu_vec_ctd_sp\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 9\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_concealed_trait_sp_try1.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1:4, 7:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[2],3), rep(strt_mu[2],1), rep(strt_q[2],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_ctd\n",
    "idparslist[[2]] <- mu_vec_ctd_sp\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 9\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_concealed_trait_sp_try2.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1:4, 7:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[3],3), rep(strt_mu[3],1), rep(strt_q[3],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_ctd\n",
    "idparslist[[2]] <- mu_vec_ctd_sp\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 9\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_concealed_trait_sp_try3.rds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extinction rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "idparsopt <- c(1, 4:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[1],1), rep(strt_mu[1],3), rep(strt_q[1],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_cr\n",
    "idparslist[[2]] <- mu_vec_ctd\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 9\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_concealed_trait_mu_try1.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "idparsopt <- c(1, 4:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[2],1), rep(strt_mu[2],3), rep(strt_q[2],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_cr\n",
    "idparslist[[2]] <- mu_vec_ctd\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 9\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_concealed_trait_mu_try2.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1, 4:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[3],1), rep(strt_mu[3],3), rep(strt_q[3],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_cr\n",
    "idparslist[[2]] <- mu_vec_ctd\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 9\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_concealed_trait_mu_try3.rds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Net diversification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "idparsopt <- c(1:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[1],3), rep(strt_mu[1],3), rep(strt_q[1],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_ctd\n",
    "idparslist[[2]] <- mu_vec_ctd\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 11\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_concealed_trait_net_diversification_try1.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[2],3), rep(strt_mu[2],3), rep(strt_q[2],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_ctd\n",
    "idparslist[[2]] <- mu_vec_ctd\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 11\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_concealed_trait_net_diversification_try2.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1:11) # our maximum rate parameter was 3\n",
    "idparsfix <- c(0) # we want to keep all zeros at zero\n",
    "initparsopt <- c(rep(strt_lambda[3],3), rep(strt_mu[3],3), rep(strt_q[3],5))\n",
    "initparsfix <- c(0.0) # all zeros remain at zero.\n",
    "idparslist <- list()\n",
    "idparslist[[1]] <- lambda_list_ctd\n",
    "idparslist[[2]] <- mu_vec_ctd\n",
    "idparslist[[3]] <- q_matrix\n",
    "\n",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 3,\n",
    "                              idparslist = idparslist,\n",
    "                              idparsopt = idparsopt,\n",
    "                              initparsopt = initparsopt,\n",
    "                              idparsfix = idparsfix,\n",
    "                              parsfix = initparsfix,\n",
    "                              sampling_fraction = f,\n",
    "                              verbose = FALSE,\n",
    "                              num_threads =  8)\n",
    "model$k <- 11\n",
    "saveRDS(model, \"SecSSE_Results/Diet_results/diet_concealed_trait_net_diversification_try3.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
}