SecSSE_Reproduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SecSSE Reproduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "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",
    "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": [
    "### Optional step (cleaning 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)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "states2<-states[!names(states) %in% setdiff(names(states), phy$tip.label)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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": [
    "### Computing sampling fractions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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": [
    "### 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))\n",
    "spec_matrix <- rbind(spec_matrix, c(3, 3, 3, 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, 3),\n",
    "                                          num_concealed_states = 4,\n",
    "                                          transition_matrix = spec_matrix,\n",
    "                                          model = \"CR\")\n",
    "lambda_list_cr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "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",
    "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, 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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "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",
    "q_matrix_cr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "\n",
    "idparsopt <- c(1:4) # 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],2))\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 = 4,\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 <- 4\n",
    "saveRDS(model, \"SecSSE_Results/Reproduction_results/reproduction_model_try1.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "\n",
    "idparsopt <- c(1:4) # 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],2))\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 = 4,\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 <- 4\n",
    "saveRDS(model, \"SecSSE_Results/Reproduction_results/reproduction_model_try2.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "\n",
    "idparsopt <- c(1:4) # 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],2))\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 = 4,\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 <- 4\n",
    "saveRDS(model, \"SecSSE_Results/Reproduction_results/reproduction_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))\n",
    "spec_matrix <- rbind(spec_matrix, c(3, 3, 3, 4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "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",
    "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, 3),\n",
    "                                   num_concealed_states = 4,\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, 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))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "q_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Speciation rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_vec_etd_sp<-rep(5, length(mu_vec_etd))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "idparsopt <- c(1:5, 9:10) # 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],4), rep(strt_mu[1],1), rep(strt_q[1],2))\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 = 4,\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/Reproduction_results/reproduction_examined_trait_sp_try1.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1:5, 9:10) # 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],4), rep(strt_mu[2],1), rep(strt_q[2],2))\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 = 4,\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/Reproduction_results/reproduction_examined_trait_sp_try2.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1:5, 9:10) # 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],4), rep(strt_mu[3],1), rep(strt_q[3],2))\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 = 4,\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/Reproduction_results/reproduction_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, 5:10) # 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],4), rep(strt_q[1],2))\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 = 4,\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/Reproduction_results/reproduction_examined_trait_mu_try1.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1, 5:10) # 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],4), rep(strt_q[2],2))\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 = 4,\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/Reproduction_results/reproduction_examined_trait_mu_try2.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1, 5:10) # 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],4), rep(strt_q[3],2))\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 = 4,\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/Reproduction_results/reproduction_examined_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:10) # 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],4), rep(strt_mu[1],4), rep(strt_q[1],2))\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 = 4,\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 <- 10\n",
    "saveRDS(model, \"SecSSE_Results/Reproduction_results/reproduction_examined_trait_net_diversification_try1.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1:10) # 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],4), rep(strt_mu[2],4), rep(strt_q[2],2))\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 = 4,\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 <- 10\n",
    "saveRDS(model, \"SecSSE_Results/Reproduction_results/reproduction_examined_trait_net_diversification_try2.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1:10) # 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],4), rep(strt_mu[3],4), rep(strt_q[3],2))\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 = 4,\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 <- 10\n",
    "saveRDS(model, \"SecSSE_Results/Reproduction_results/reproduction_examined_trait_net_diversification_try3.rds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conceal trait"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "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, 3),\n",
    "                                   num_concealed_states = 4,\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(5, length(mu_vec_ctd))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "idparsopt <- c(1:5, 9:10) # 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],4), rep(strt_mu[1],1), rep(strt_q[1],2))\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 = 4,\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/Reproduction_results/reproduction_concealed_trait_sp_try1.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1:5, 9:10) # 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],4), rep(strt_mu[2],1), rep(strt_q[2],2))\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 = 4,\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/Reproduction_results/reproduction_concealed_trait_sp_try2.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1:5, 9:10) # 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],4), rep(strt_mu[3],1), rep(strt_q[3],2))\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 = 4,\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/Reproduction_results/reproduction_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, 5:10) # 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],4), rep(strt_q[1],2))\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 = 4,\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/Reproduction_results/reproduction_concealed_trait_mu_try1.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1, 5:10) # 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],4), rep(strt_q[2],2))\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 = 4,\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/Reproduction_results/reproduction_concealed_trait_mu_try2.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1, 5:10) # 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],4), rep(strt_q[3],2))\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 = 4,\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/Reproduction_results/reproduction_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:10) # 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],4), rep(strt_mu[1],4), rep(strt_q[1],2))\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 = 4,\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 <- 10\n",
    "saveRDS(model, \"SecSSE_Results/Reproduction_results/reproduction_concealed_trait_net_diversification_try1.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1:10) # 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],4), rep(strt_mu[2],4), rep(strt_q[2],2))\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",
    "model <- secsse::cla_secsse_ml(phy = phy,\n",
    "                              traits = traits,\n",
    "                              num_concealed_states = 4,\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 <- 10\n",
    "saveRDS(model, \"SecSSE_Results/Reproduction_results/reproduction_concealed_trait_net_diversification_try2.rds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1:10) # 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],4), rep(strt_mu[3],4), rep(strt_q[3],2))\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 = 4,\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 <- 10\n",
    "saveRDS(model, \"SecSSE_Results/Reproduction_results/reproduction_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
}