SecSSE is an R package designed for multistate data sets under a concealed state and speciation (‘hisse’) framework. In this sense, it is parallel to the ‘MuSSE’ functionality implemented in ‘diversitree’, but it accounts for finding possible spurious relationships between traits and diversification rates (‘false positives’, Rabosky & Goldberg 2015) by testing against a ‘hidden trait’ (Beaulieu et al. 2013), which is responsible for more variation in diversification rates than the trait being investigated.
A good practice is always remove all the objects in memory and then load SecSSE:
rm(list = ls())
library(secsse)Similar to the ‘diversitree’ (Fitzjohn et al. 2012) and ‘hisse’ (Beaulieu & O’Meara 2016) packages, SecSSE uses two input files: a rooted, ultrametric tree in nexus format (for conversion of other formats to nexus, we refer to the documentation in package ‘ape’) and a data file with two columns, the first containing taxa names and the second a numeric code for trait state with a header (usually 0,1,2,3, etc., but notice that ‘NA’ is a valid code too, if you are not sure what trait state to assign to a taxon). A comma-separated value file (.csv) generated in MsExcel works particularly well. The *.csv file can be loaded into R using the read.csv() function. and should look like this:
data(traitinfo)
trait <- traitinfo
tail(trait)## species states
## 171 out_171 2
## 172 out_172 3
## 173 out_173 2
## 174 out_174 2
## 175 out_175 3
## 176 out_176 1This data set (here we see only the bottom lines of the data frame) has three character states labeled as 1, 2 and 3. Notice that unless you want to assign ambiguity to some but not all states (see below), the third column in your data file should be empty. Ambiguity about trait state (you are not sure which trait state to assign a taxon too, or you have no data on trait state for a particular taxon), can be assigned using ‘NA’. SecSSE handles ‘NA’ differently from a full trait state, in that it assigns probabilities to all trait states for a taxon demarcated with ‘NA’.
The second object we need is an ultrametric phylogenetic tree, that is rooted and has labeled tips. One can load it in R by using read.nexus(). In our example we load a prepared phylogeny named “phylo_Vign”:
data("phylo_Vign")For running SecSSE it is important that tree tip labels agree with taxon names in the data file, but also that these are in the same order. For this purpose, we run the following piece of code prior to any analysis:
traits <- sortingtraits(trait,phylo_Vign)If there is a mismatch in the number of taxa between data and tree file, you will receive an error message. However, to then identify which taxa are causing issues and if they are in the tree or data file, you can use the name.check function in the ‘geiger’(Harmon et al. 2008) package:
library(geiger)## Loading required package: ape#making sure that the first line is identified as containing header info:
rownames(trait) <- trait[,1]
#pick out all elements that do not agree between tree and data
mismat <- name.check(phylo_Vign,trait)
#this will call all taxa that are in the tree, but not the data file
#mismat$tree_not_data
#and conversely,
#mismat$data_not_treeIf you have taxa in your tree file that do not appear in your trait file, it is worth adding them with value ‘NA’ for trait state. After you are done properly setting up your data, you can proceed to setting parameters and constraints.
##Parameter settings and constraints SecSSE allows for the implementation of different models of evolution, and just as in ‘diversitree’ and ‘hisse’, parameters can be fixed at certain values (if prior information is known on particular values) or made to be equal to each other. Initial parameter values can also be supplied, to start off the maximum likelihood search with. The main function in the SecSSE package is secsse_ml, which performs a maximum likelihood search and uses as input a set of speciation rate parameters (lambda), a set of extinction rate parameters (mu), and a matrix composed of transition rates (q) between the various states. The identifiers of the parameters are broadly the same as those used in ‘hisse’, and numbers indicate examined state, whereas letters denote concealed state, 2A for example being in examined state 2, and concealed state A.
Both the speciation and extinction parameters are supplied as vectors, and the transition rates are supplied as a matrix, joined in a list.
The function secsse_ml takes the following arguments, PHY, TRAITS, NUM_CONCEALED_STATES, IDPARSLIST, INITPARSOPT, IDPARSOPT, IDPARSFIX, PARSFIX, COND,WEIGHTTRAITS, SAMPLING_FRACTION, TOL, METHODE, OPTIMMETHOD, and bigtree. These are best declared outside of the secsse_ml function, then called in the function. We discuss these here chronologically:
PHY: a user-supplied phylogenetic tree of class ‘phylo’ (see above)
TRAITS: user-supplied trait data of class ‘data frame’ (see above)
NUM_CONCEALED_STATES: In general, we recommend this value to be equal to the number of examined states in your data set (that way they have the same parametric complexity), however , this may or may not be computationally tractable depending on the size of the tree. An alternative is to set this value to 3, an advantage of having just three concealed states is that data interpretation gets a lot easier. Notice this value needs to be specified also under id_paramPos.
IDPARSLIST: a list of parameters to be supplied to the function. This list contains information on the number of parameters, as well as on how parameters interact. E.g., if we would like all speciation rates to behave similarly or if we want two transition rates to be identical, we can set this up here. Setting the parameters for this argument is the bulk of the work for setting up the model, especially when the number of states is relatively high.
The following is a visual example of the input parameters of SecSSE, this list is composed of three elements, one containing the number of parameters for lambda, one for mu, and a matrix of transition rates. Notice that this list contains the set-up for a model in which all parameters are free, every parameter having a unique value, indicating that each parameter is optimized separately. The diagonal of the q matrix is always set to NA, because transitions within a state are not possible. The dimensions of the transition matrix follow the following rule: (3n)2, where n is the number of observed states. Needless to say, we have not tried running SecSSE with n>10, both for computational and practical reasons, and neither should you probably, especially in combination with large trees.
#First we have to define idparslist, as well as, again, a user-specified value for the number of concealed states to be assessed by SecSSE.
idparslist <- id_paramPos(traits, num_concealed_states = 3)
#Let's take a look at the full all-free model by now simply typing
idparslist## $lambdas
## 1A 2A 3A 1B 2B 3B 1C 2C 3C
## 1 2 3 4 5 6 7 8 9
##
## $mus
## 1A 2A 3A 1B 2B 3B 1C 2C 3C
## 10 11 12 13 14 15 16 17 18
##
## $Q
## 1A 2A 3A 1B 2B 3B 1C 2C 3C
## 1A NA 19 20 21 22 23 24 25 26
## 2A 27 NA 28 29 30 31 32 33 34
## 3A 35 36 NA 37 38 39 40 41 42
## 1B 43 44 45 NA 46 47 48 49 50
## 2B 51 52 53 54 NA 55 56 57 58
## 3B 59 60 61 62 63 NA 64 65 66
## 1C 67 68 69 70 71 72 NA 73 74
## 2C 75 76 77 78 79 80 81 NA 82
## 3C 83 84 85 86 87 88 89 90 NAIf we would like the speciation rate of states 1B and 2B to be the same, we can do this as follows:
#idparslist[[1]][c(5,6)] <- 5 Notice that if one were to set extinction parameters to be the same, the numbering used to identify parameters is not the same as that in idparslist, but rather consecutive numbering referring to the elements within the extinction parameters component in idparslist:
#idparslist[[2]][c(1:9)] <- 7There are also several things we can do to improve the rate matrix and reduce its computational complexity. First of all, we should leave transitions between the same state out of the calculations with a simple command that orders all values on the diagonal of the matrix not to be calculated. This is included as a default within idparslist, but after modifying the q matrix in any way, it is a good idea to ensure that the diagonals are still not included in the calculations:
diag(idparslist[[3]]) <- NAAdditionally, we would like to set all dual transitions (so for example from state 0 to 1 AND from concealed state A to B) to 0, as these are unlikely to occur. It is a bit of a matter of personal preference whether or not you should do this, but we follow Beaulieu & O’Meara (2016) here and set dual transitions to zero. One good reason for doing so is simply to reduce computational burden.
idparslist[[3]][1,c(5,6,8,9)] <- 0
idparslist[[3]][2,c(4,6,7,9)] <- 0
idparslist[[3]][3,c(4,5,7,8)] <- 0
idparslist[[3]][4,c(2,3,8,9)] <- 0
idparslist[[3]][5,c(1,3,7,9)] <- 0
idparslist[[3]][6,c(1,2,7,8)] <- 0
idparslist[[3]][7,c(2,3,5,6)] <- 0
idparslist[[3]][8,c(1,3,4,6)] <- 0
idparslist[[3]][9,c(1,2,4,5)] <- 0These three actions together then yield the following:
idparslist## $lambdas
## 1A 2A 3A 1B 2B 3B 1C 2C 3C
## 1 2 3 4 5 6 7 8 9
##
## $mus
## 1A 2A 3A 1B 2B 3B 1C 2C 3C
## 10 11 12 13 14 15 16 17 18
##
## $Q
## 1A 2A 3A 1B 2B 3B 1C 2C 3C
## 1A NA 19 20 21 0 0 24 0 0
## 2A 27 NA 28 0 30 0 0 33 0
## 3A 35 36 NA 0 0 39 0 0 42
## 1B 43 0 0 NA 46 47 48 0 0
## 2B 0 52 0 54 NA 55 0 57 0
## 3B 0 0 61 62 63 NA 0 0 66
## 1C 67 0 0 70 0 0 NA 73 74
## 2C 0 76 0 0 79 0 81 NA 82
## 3C 0 0 85 0 0 88 89 90 NANotice that all entries in the lambda and mu vectors, as well as the rate matrix should be supplied to either idparsopt or idparsfix, including the zeros that represent dual transitions (which are supplied to idparsfix and set to zero under parsfix.
Numbers in all elements of the list can be skipped without a problem, as long as they are supplied correctly to other arguments. When Q-matrices get larger, it can be good to specify all values in the matrices separately and consecutively (no matter how laborious), for reasons of intuition. This can facilitate setting up idparsopt, idparsfix, and initparsopt, as well as help in setting up different models (using various combinations of parameter constraints) along the way. Here is a piece of code that can be copied for a 3-state analysis:
idparslist[[3]][1,c(2)] <- 19
idparslist[[3]][1,c(3)] <- 20
idparslist[[3]][1,c(4)] <- 21
idparslist[[3]][1,c(7)] <- 22
idparslist[[3]][1,c(5,6,8,9)] <- 0
idparslist[[3]][2,c(1)] <- 23
idparslist[[3]][2,c(3)] <- 24
idparslist[[3]][2,c(5)] <- 25
idparslist[[3]][2,c(8)] <- 26
idparslist[[3]][2,c(4,6,7,9)] <- 0
idparslist[[3]][3,c(1)] <- 27
idparslist[[3]][3,c(2)] <- 28
idparslist[[3]][3,c(6)] <- 29
idparslist[[3]][3,c(9)] <- 30
idparslist[[3]][3,c(4,5,7,8)] <- 0
idparslist[[3]][4,c(1)] <- 31
idparslist[[3]][4,c(5)] <- 32
idparslist[[3]][4,c(6)] <- 33
idparslist[[3]][4,c(7)] <- 34
idparslist[[3]][4,c(2,3,8,9)] <- 0
idparslist[[3]][5,c(2)] <- 35
idparslist[[3]][5,c(4)] <- 36
idparslist[[3]][5,c(6)] <- 37
idparslist[[3]][5,c(8)] <- 38
idparslist[[3]][5,c(1,3,7,9)] <- 0
idparslist[[3]][6,c(3)] <- 39
idparslist[[3]][6,c(4)] <- 40
idparslist[[3]][6,c(5)] <- 41
idparslist[[3]][6,c(9)] <- 42
idparslist[[3]][6,c(1,2,7,8)] <- 0
idparslist[[3]][7,c(1)] <- 43
idparslist[[3]][7,c(4)] <- 44
idparslist[[3]][7,c(8)] <- 45
idparslist[[3]][7,c(9)] <- 46
idparslist[[3]][7,c(2,3,5,6)] <- 0
idparslist[[3]][8,c(2)] <- 47
idparslist[[3]][8,c(5)] <- 48
idparslist[[3]][8,c(7)] <- 49
idparslist[[3]][8,c(9)] <- 50
idparslist[[3]][8,c(1,3,4,6)] <- 0
idparslist[[3]][9,c(3)] <- 51
idparslist[[3]][9,c(6)] <- 52
idparslist[[3]][9,c(7)] <- 53
idparslist[[3]][9,c(8)] <- 54
idparslist[[3]][9,c(1,2,4,5)] <- 0
diag(idparslist[[3]]) <- NAThis yields the following data setup:
idparslist## $lambdas
## 1A 2A 3A 1B 2B 3B 1C 2C 3C
## 1 2 3 4 5 6 7 8 9
##
## $mus
## 1A 2A 3A 1B 2B 3B 1C 2C 3C
## 10 11 12 13 14 15 16 17 18
##
## $Q
## 1A 2A 3A 1B 2B 3B 1C 2C 3C
## 1A NA 19 20 21 0 0 22 0 0
## 2A 23 NA 24 0 25 0 0 26 0
## 3A 27 28 NA 0 0 29 0 0 30
## 1B 31 0 0 NA 32 33 34 0 0
## 2B 0 35 0 36 NA 37 0 38 0
## 3B 0 0 39 40 41 NA 0 0 42
## 1C 43 0 0 44 0 0 NA 45 46
## 2C 0 47 0 0 48 0 49 NA 50
## 3C 0 0 51 0 0 52 53 54 NAINITPARSOPT: user-supplied values of parameters, a vector of values of lambda, mu, and q that should agree in number with the number of parameters specified in the model. If values are known beforehand, they can be specified as follows for the case of the above defined parameter set, where there are 5 lambda’s (two equal), 6 mu’s (all free), and q’s (all free, but no dual transitions):
initparsopt <- c(rep(1.2,9), rep(0.1,9), rep(0.25,36))IDPARSOPT: the id’s of the parameters we want to optimize (versus those that are to be fixed). The id’s should correspond to those specified under idparslist. For example, if we take our previously defined idparslist:
idparslist## $lambdas
## 1A 2A 3A 1B 2B 3B 1C 2C 3C
## 1 2 3 4 5 6 7 8 9
##
## $mus
## 1A 2A 3A 1B 2B 3B 1C 2C 3C
## 10 11 12 13 14 15 16 17 18
##
## $Q
## 1A 2A 3A 1B 2B 3B 1C 2C 3C
## 1A NA 19 20 21 0 0 22 0 0
## 2A 23 NA 24 0 25 0 0 26 0
## 3A 27 28 NA 0 0 29 0 0 30
## 1B 31 0 0 NA 32 33 34 0 0
## 2B 0 35 0 36 NA 37 0 38 0
## 3B 0 0 39 40 41 NA 0 0 42
## 1C 43 0 0 44 0 0 NA 45 46
## 2C 0 47 0 0 48 0 49 NA 50
## 3C 0 0 51 0 0 52 53 54 NAAnd we want to optimize only speciation rate parameters, while keeping the rest fixed, we specify the following:
idparsopt <- c(1:9)In this case, values must be provided for the extinction parameters and transition rate matrix under parsfix, and their corresponding numbers must be identified under idparsfix.
Another example:
#this would optimize speciation and extinction in the above setup
#idparsopt <- c(1:18)Often what we will want to do is to make all transition rates equal. Or define that all extinctions are the same. We first define our parameter list as follows:
idparslist[[2]][] <- 10
idparslist[[3]][1,c(2,3,4,7)] <- 11
idparslist[[3]][1,c(5,6,8,9)] <- 0
idparslist[[3]][2,c(1,3,5,8)] <- 11
idparslist[[3]][2,c(4,6,7,9)] <- 0
idparslist[[3]][3,c(1,2,6,9)] <- 11
idparslist[[3]][3,c(4,5,7,8)] <- 0
idparslist[[3]][4,c(1,5,6,7)] <- 11
idparslist[[3]][4,c(2,3,8,9)] <- 0
idparslist[[3]][5,c(2,4,6,8)] <- 11
idparslist[[3]][5,c(1,3,7,9)] <- 0
idparslist[[3]][6,c(3,4,5,9)] <- 11
idparslist[[3]][6,c(1,2,7,8)] <- 0
idparslist[[3]][7,c(1,4,8,9)] <- 11
idparslist[[3]][7,c(2,3,5,6)] <- 0
idparslist[[3]][8,c(2,5,7,9)] <- 11
idparslist[[3]][8,c(1,3,4,6)] <- 0
idparslist[[3]][9,c(3,6,7,8)] <- 11
idparslist[[3]][9,c(1,2,4,5)] <- 0
diag(idparslist[[3]]) <- NAThen we will optimize speciation and the single transition rate:
idparsopt <- c(1:9,11)IDPARSFIX: the id’s of parameters we want fixed at a certain value (including zero).Notice that 0 in idparslist is just another ID. Parallel to idparsopt, the following statement would fix all parameters associated with extinction rates:
idparsfix <- c(0,10)Notice that if dual transitions were set to zero under idparslist, we should do this here too.
PARSFIX: specifies at which values the parameters identified under idparsfix should be set. Should have the same number of entries as idparsfix (same order too). In this example, the first zero means that all those entries in idparslist with ID 0 will be fixed to zero. The second zero means that all the entries in idparslist with ID 10, will be fixed to 0.0001.
parsfix <- c(0,0.0001)One can also estimate initial lambda and mu values from the tree using a simple birth-death model that does not take into account trait states. Here we do this with the bd_ML function from the DDD package. A good starting point for q is lambda/5:
library(DDD)
startingpoint <- bd_ML(brts = ape::branching.times(phylo_Vign))## You are optimizing lambda0 mu0
## You are fixing lambda1 mu1
## Optimizing the likelihood - this may take a while.
## The loglikelihood for the initial parameter values is -657.65168907669.
##
## Maximum likelihood parameter estimates: lambda0: 0.066492, mu0: 0.000115, lambda1: 0.000000, mu1: 0.000000:
## Maximum loglikelihood: -645.684269intGuessLamba <- startingpoint$lambda0
intGuessMu <- startingpoint$mu0
#Make sure that the dimensions of initparsopt agree with those of idparslist and idparsopt, especially in the case of the initial guesses for rates supplied to the Q matrix. Rules of thumb are that if n=number of examined states, both intGuessLamba and intGuessMu should be replicated 2n times, and (intGuessLamba/5) should be replicated (2n)2/2 times:
initparsopt <- c(rep(intGuessLamba,9), rep((intGuessLamba/5),1))COND: conditioning on the state of the root. Set to “maddison_cond” if you want conditioning as done in other -SSE packages, or “proper_cond” if you want to use our new improved conditioning.
root_state_weight: SecSSe offers to methods to weigh the probabilities of states at the root:“proper_weights” and “maddison_weights”. In the accompanying paper you can read the differences between them.
SAMPLING_FRACTION: include a sampling fraction. Sampling.f always has as many elements as there are examined states, so a SecSSE analysis with 3 states could have the following sampling_fraction = c(0.5,0.25,0.75), in which half of taxa in state 1 are sampled, a quarter in state two, and three quarters in state three. If 100% of known taxa in each state are sampled, sampling_fraction=c(1,1,1). If only an overall value is known (for example, we know we sampled 80% of all taxa, but we do not know how they are distributed across states), we assign this value to each state: sampling_fraction = c(0.8,0.8,0.8). Sampling.f is always placed after the ‘cond’ statement.
TOL: basically, a range of values between which samples in the ML chain will be accepted or not. Typically, the value of tol = c(1e-04, 1e-05, 1e-07) is generally best.
METHODE: method for integration of likelihood values along branches, generally we recommend “ode45”.
OPTIMMETHOD: optimization method, generally we recommend “simplex”.
RUN_PARALLEL: this specifies whether or not to use the SecSSE tree-breaking function. If you have a large tree, this tree can be broken into two pieces so that computation of likelihood along branches can take place simultaneously on the two pieces, yielding a gain in computation time. The size of the two pieces is established by SecSSE, and depends on how balanced the tree is; a better-balanced tree yields two pieces of relatively equal size and results in relatively larger gain in computation time. With large trees (say, n>1000), it is our experience that even two chunks of tree of unequal size yield a time advantage. Needless to say, your computational setup needs to be able to accommodate parallel computation (multiple cores, nodes).
##Running the likelihood function
After we have defined all of the necessary parameters for running secsse_ml, we can start running our analysis:
#secsse_ml(phylo_Vign,traits, num_concealed_states=3, idparslist, idparsopt, initparsopt, idparsfix, parsfix, cond="maddison_cond",root_state_weight = "maddison_weights", tol = c(1e-04, 1e-05, 1e-07), sampling_fraction=c(1,1,1), maxiter = 1000 * round((1.25)^length(idparsopt)), use_fortran=TRUE,methode="ode45", optimmethod = "simplex", num_cycles = 1, run_parallel=FALSE)We can save output to an R data file, for example, here called output.RDS:
#out<-secsse_ml(phylo_Vign,traits, num_concealed_states=3, idparslist, idparsopt, initparsopt, idparsfix, parsfix, cond="maddison_cond",root_state_weight = "maddison_weights", sampling_fraction=c(1,1,1), tol = c(1e-04, 1e-05, 1e-07), maxiter = 1000 * round((1.25)^length(idparsopt)), use_fortran=TRUE,methode="ode45", optimmethod = "simplex", num_cycles = 1, run_parallel=FALSE)
#saveRDS(out, file="output.RDS")Later on, we can retrieve the data in this file, simply by entering:
#readRDS("output.RDS")The following is sample output, with two concealed states, notice in this case all transition rates, including dual rates, were set to the fixed value of 0.01:
#$MLpars[[1]]
# 1A 2A 3A 1B 2B 3B
#4.842634e-16 1.080409e-01 7.843821e-02 4.029147e-09 3.018863e-02 3.018863e-02
#$MLpars[[2]]
# 1A 2A 3A 1B 2B 3B
#0.002000000 0.002000109 0.002734071 0.001988593 0.002169052 0.003969142
#$MLpars[[3]]
# 1A 2A 3A 1B 2B 3B
#1A NA 0.01 0.01 0.01 0.01 0.01
#2A 0.01 NA 0.01 0.01 0.01 0.01
#3A 0.01 0.01 NA 0.01 0.01 0.01
#1B 0.01 0.01 0.01 NA 0.01 0.01
#2B 0.01 0.01 0.01 0.01 NA 0.01
#3B 0.01 0.01 0.01 0.01 0.01 NA
#$ML
#[1] -848.0895The maximum likelihood value at the bottom of the output can be used in model comparison.
##SecSSE tool to facilitate composition of q matrices Often, q matrices can get quite large and complicated, the more states you are analyzing. We have devised a tool to more easily put together q matrices. This tool starts from the so-called ‘masterBlock’, the basic matrix in which we only find information on transitions between examined states. The information contained in this ‘masterBlock’ is then automatically mimicked for inclusion in the full matrix, to ensure that the same complexity in examined state transitions is also found in concealed states. The use of the ‘masterBlock’ implies that you are using the same number of concealed as examined states. Here, we are generating a ‘masterBlock’ that yields a 3-state q matrix.
The ‘masterBlock’ can be declared as follows:
masterBlock<-matrix(99,ncol=3,nrow=3,byrow=T) in which ‘99’ is an example value you can use to populate the matrix at first, to be replaced by values you specify. If you make this value conspicuously different from others, you can ensure that you are not skipping the specification of values, as any non-specified rates will take this value. ‘Ncol’ and ‘nrow’ will need to reflect the number of states you are analyzing.
We first declare all values on the diagonal to be ‘NA’, then we specify values for the ‘masterBlock’. The values have a row and column indicator, so that e.g. ‘[2,7]’ refers to position 7 in row 2, or to a transition from 2A to 7A more specifically.
diag(masterBlock) <- NA
masterBlock[1,2] <- 6
masterBlock[1,3] <- 7
masterBlock[2,1] <- 8
masterBlock[2,3] <- 9
masterBlock[3,1] <- 10
masterBlock[3,2] <- 11After completing the declaration of the ‘masterBlock’, we will need to specify whether or not we want the variation in examined states to be exactly the same as in the concealed state (so that e.g. the transition 1A->3A takes the same value as 5A->5C), or if we want the concealed state to have additional variation to account for type I error in transition rates (so that the total amount of transition parameters between concealed states is the same as between examined states, but the values are different). This is done by:
diff.conceal <- FALSEFinally, we need to make sure the ‘masterBlock’ is used as a baseline for building the transition matrix in IDPARSLIST:
myQ<-q_doubletrans(traits,masterBlock,diff.conceal)
idparslist[[3]] <- myQWhich makes our final q matrix look as follows:
idparslist[[3]]## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] NA 6 7 6 0 0 7 0 0
## [2,] 8 NA 9 0 6 0 0 7 0
## [3,] 10 11 NA 0 0 6 0 0 7
## [4,] 8 0 0 NA 6 7 9 0 0
## [5,] 0 8 0 8 NA 9 0 9 0
## [6,] 0 0 8 10 11 NA 0 0 9
## [7,] 10 0 0 11 0 0 NA 6 7
## [8,] 0 10 0 0 11 0 8 NA 9
## [9,] 0 0 10 0 0 11 10 11 NAMatching the amount of variation in rates between the concealed states, yields the following:
diff.conceal <- TRUE
myQ <- q_doubletrans(traits,masterBlock,diff.conceal)
idparslist[[3]] <- myQ
idparslist[[3]]## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] NA 6 7 12 0 0 13 0 0
## [2,] 8 NA 9 0 12 0 0 13 0
## [3,] 10 11 NA 0 0 12 0 0 13
## [4,] 14 0 0 NA 6 7 15 0 0
## [5,] 0 14 0 8 NA 9 0 15 0
## [6,] 0 0 14 10 11 NA 0 0 15
## [7,] 16 0 0 17 0 0 NA 6 7
## [8,] 0 16 0 0 17 0 8 NA 9
## [9,] 0 0 16 0 0 17 10 11 NANote that in this case, the number of transition rate parameters doubles compared to the previous example, a change that will need to be applied also to IDPARSOPT and IDPARSFIX.
##SecSSE function to reduce number of transition rate parameters by including multiplicative factors SecSSE has the capability of reducing computational burden by decreasing the number of transition rate parameters through the inclusion of multiplicative factors. Factors can also be used to disentangle complex patterns of trait-dependent diversification when multiple traits are included.
Suppose you are running an analysis with a large number of transition rate parameters, but you suspect there are linear relationships between some of them. If the transition between lobed (L) and palmate (P) feet is twice as infrequent as that between palmate and semi-palmate (S) feet, and could say that P->L is 2(P->S). The reverse would also be true: L->P is 2(S->P). By applying these factors, we are reducing the transition matrix from a 6 parameters to 4, and in models where transitions between concealed states are allowed, we are reducing our parameters from 12 to 8. Of course, the inclusion of these factors comes with a loss of resolution, and is therefore best done with parameters where exact estimation is not essential.
In SecSSE, the factors are represented in a function separate from secsse_ml, and the setup of this function is very similar to secsse_ml, but requires the addition of two parameters, SHAREFACTORS and INITFACTORS.
SHAREFACTORS: these are the identifiers of the factors you want to specify. In the above example, we have two factors, one governing transitions from P->S and one from S->P. Transitions in opposite directions are better not fixed to the same multiplicative factor, so that at least two are needed here. In this case these are specified as follows:
#shareFactors <- c(.1,.2)INITFACTORS: Since these shared factors need initial parameter estimates, just as other transition parameters in the model do, we need to specify these. The initial guesses are best set to 1, so they behave similar to the parameters they are ‘tied’ to, unless we have very good evidence (e.g. from a previous run) that these are bigger or smaller:
#initFactors <- c(1,1)Aside from setting these two parameters, we need to specify in our rate matrix which rate parameters we want to be governed by which factors. Imagine we have a 3-state matrix, where 1 refers to lobed feet, 2 to semi-palmate and 3 to palmate:
# diag(masterBlock) <- NA
# masterBlock[1,2] <- 6
# masterBlock[1,3] <- 6.1 #factor 1: lobed to palmate
#
# masterBlock[2,1] <- 7
# masterBlock[2,3] <- 8
#
# masterBlock[3,1] <- 7.2 #factor 2: palmate to lobed
# masterBlock[3,2] <- 9Finally, we run the function secsse_ml_struc instead of secsse_ml, and make sure that both new parameters are included.
#secsse_ml_struc(phylo_Vign..., shareFactors, initFactors)Multiplicative factors can also be used in connection with lambdas or mus, in the same way as they are used for transition rates. Note that in such case the factors will need to be unique across the entire dataset, so that both speciation- and transition-related factors have unique values for shareFactors. They can also be used to disentangle complex patterns of diversification when multiple traits are taken into account. Assume that aside from foot shape (the above example), we are also looking at the presence or absence of a spur, and we would like to know how the two traits interact to influence diversification. In such a case, presence or absence of spur can be used as a multiplicative factor, and models can be run where presence or absence is coded as the same multiplicative factor (.1), and where presence or absence are coded as two different factors (.1,.2).
##Note on assigning ambiguity to taxon trait states If the user wishes to assign a taxon to multiple trait states, because he/she is unsure which state best describes the taxon, he/she can use ‘NA’. ‘NA’ is used when there is no information on possible state at all; for example when a state was not measured or a taxon is unavailable for inspection. ‘NA’ means a taxon is equally likely to pertain to any state. In case the user does have some information, for example if a taxon can pertain to multiple states, or if there is uncertainty regarding state but one or multiple states can with certainty be excluded, SecSSE offers flexibility to handle ambiguity. In this case, the user only needs to supply a trait file, with at least four columns, one for the taxon name, and three for trait state. Below, we show an example of what the trait info should be like (the column with species’ names has been removed).If a taxon may pertain to trait state 1 or 3, but not to 2, the three columns should have at least the values 1 and a 3, but never 2 (species in the third row). On the other hand, the species in the fifth row can pertain to all states: the first column would have a 1, the second a 2, the third a 3 (although if you only have this type of ambiguity, it is easier to assign ‘NA’ and use a single-column data file).
# traits traits traits
# [1,] 2 2 2
# [2,] 1 1 1
# [3,] 2 2 2
# [4,] 3 1 1
# [5,] 1 2 3##Do you feel SecSSE? If not, please feel free to e-mail the authors. For help with this R package only.
##References Beaulieu, J. M., O’meara, B. C., & Donoghue, M. J. (2013). Identifying hidden rate changes in the evolution of a binary morphological character: the evolution of plant habit in campanulid angiosperms. Systematic biology, 62(5), 725-737.
Beaulieu, J. M., & O’Meara, B. C. (2016). Detecting hidden diversification shifts in models of trait-dependent speciation and extinction. Systematic biology, 65(4), 583-601.
FitzJohn, R. G. (2012). Diversitree: comparative phylogenetic analyses of diversification in R. Methods in Ecology and Evolution, 3(6), 1084-1092.
Harmon, L. J., Weir, J. T., Brock, C. D., Glor, R. E., & Challenger, W. (2008). GEIGER: investigating evolutionary radiations. Bioinformatics, 24(1), 129-131.
Rabosky, D. L., & Goldberg, E. E. (2015). Model inadequacy and mistaken inferences of trait-dependent speciation. Systematic Biology, 64(2), 340-355.