SecSSE_Size_replicated.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SecSSE Size"
   ]
  },
  {
   "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\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading data\n",
    "\n",
    "args<-commandArgs(trailingOnly = TRUE)\n",
    "\n",
    "df<-read.csv(\"Species_data.tsv\", sep=\"\\t\") # omit sep =\"\\t\" for .csv files\n",
    "\n",
    "phy<-read.nexus(args[1])\n",
    "\n",
    "phy<-phy[[as.numeric(args[2])]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clusterizing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mb = Mclust(as.numeric(df$Body.size))\n",
    "summary(mb, parameters = TRUE)\n",
    "df$Body.size <- as.factor(mb$classification)"
   ]
  },
  {
   "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$Body.size\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(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": [
    "### 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, 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))"
   ]
  },
  {
   "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>=5]<-5\n",
    "q_matrix_cr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "### Try 1\n",
    "\n",
    "idparsopt <- c(1,2,3,4,5) # 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],3))\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 <- 5\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_model_try1.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "\n",
    "idparsopt <- c(1,2,3,4,5) # 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],3))\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 <- 5\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_model_try2.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "\n",
    "idparsopt <- c(1,2,3,4,5) # 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],3))\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 <- 5\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_model_try3.rds\", sep =\"\"))"
   ]
  },
  {
   "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, 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))"
   ]
  },
  {
   "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>=9]<-9\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:9) # 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],3))\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 <- 7\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_examined_trait_sp_try1.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1:4, 7:9) # 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],3))\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 <- 7\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_examined_trait_sp_try2.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1:4, 7:9) # 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],3))\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 <- 7\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_examined_trait_sp_try3.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extinction rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "idparsopt <- c(1, 4:9) # 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],3))\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 <- 7\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_examined_trait_mu_try1.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1, 4:9) # 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],3))\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 <- 7\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_examined_trait_mu_try2.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1, 4:9) # 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],3))\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 <- 7\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_examined_trait_mu_try3.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Diversification rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "idparsopt <- c(1:9) # 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],3))\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 <- 9\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_examined_trait_net_diversification_try1.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1:9) # 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],3))\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 <- 9\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_examined_trait_net_diversification_try2.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1:9) # 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],3))\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 <- 9\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_examined_trait_net_diversification_try3.rds\", sep =\"\"))"
   ]
  },
  {
   "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:9) # 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],3))\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 <- 7\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_concealed_trait_sp_try1.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1:4, 7:9) # 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],3))\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 <- 7\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_concealed_trait_sp_try2.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1:4, 7:9) # 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],3))\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 <- 7\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_concealed_trait_sp_try3.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extinction rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "idparsopt <- c(1, 4:9) # 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],3))\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 <- 7\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_concealed_trait_mu_try1.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "idparsopt <- c(1, 4:9) # 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],3))\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 <- 7\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_concealed_trait_mu_try2.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1, 4:9) # 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],3))\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 <- 7\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_concealed_trait_mu_try3.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Net diversification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 1\n",
    "idparsopt <- c(1:9) # 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],3))\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 <- 9\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_concealed_trait_net_diversification_try1.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 2\n",
    "idparsopt <- c(1:9) # 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],3))\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 <- 9\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_concealed_trait_net_diversification_try2.rds\", sep =\"\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Try 3\n",
    "idparsopt <- c(1:9) # 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],3))\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 <- 9\n",
    "saveRDS(model, paste(\"Results\", args[2], \"/Size_results/size_concealed_trait_net_diversification_try3.rds\", sep =\"\"))"
   ]
  }
 ],
 "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
}