class: title-slide, middle <style type="text/css"> .title-slide { background-image: url('assets/img/bg.jpg'); background-color: #23373B; background-size: contain; border: 0px; background-position: 600px 0; line-height: 1; } </style> # Maximum likelihood <hr width="65%" align="left" size="0.3" color="orange"></hr> ## Technical problem <hr width="65%" align="left" size="0.3" color="orange" style="margin-bottom:40px;" alt="@Martin Sanchez"></hr> .instructors[ **ECL707/807** - Dominique Gravel ] <img src="assets/img/logo.png" width="25%" style="margin-top:20px;"></img> <img src="assets/img/Bios2_reverse.png" width="23%" style="margin-top:20px;margin-left:35px"></img> --- # Distribution of biomes in Québec <div style='text-align:center;'> <img src="assets/img/bioclim.png" height="500px"></img> </div> --- # A metacommunity model of forest dynamics <div style='text-align:center;'> <img src="assets/img/modele_vissault.png" height="500px"></img> </div> --- # States - 'B' if there are only boreal tree species - 'T' if there are only temperate tree species - 'M' for mixtures - 'R' if there are no trees (regeneration) --- # Transition matrix between states $$ `\begin{bmatrix} P(B_{t+1}|B_t) & P(B_{t+1}|T_t) & P(B_{t+1}|M_t) & P(B_{t+1}|R_t) \\ P(T_{t+1}|B_t) & P(T_{t+1}|T_t) & P(T_{t+1}|M_t) & P(T_{t+1}|R_t) \\ P(M_{t+1}|B_t) & P(M_{t+1}|T_t) & P(M_{t+1}|M_t) & P(M_{t+1}|R_t) \\ P(R_{t+1}|B_t) & P(R_{t+1}|T_t) & P(R_{t+1}|M_t) & P(R_{t+1}|R_t) \\ \end{bmatrix}` $$ Where rows are states at time `\(t+1\)` and columns states at time `\(t\)`. A column is therefore a vector of probabilities corresponding to a multinomial distribution. Each column must therefore sum to 1. --- # Transition matrix between states $$ `\begin{bmatrix} P(B_{t+1}|B_t) & P(B_{t+1}|T_t) & P(B_{t+1}|M_t) & P(B_{t+1}|R_t) \\ P(T_{t+1}|B_t) & P(T_{t+1}|T_t) & P(T_{t+1}|M_t) & P(T_{t+1}|R_t) \\ P(M_{t+1}|B_t) & P(M_{t+1}|T_t) & P(M_{t+1}|M_t) & P(M_{t+1}|R_t) \\ P(R_{t+1}|B_t) & P(R_{t+1}|T_t) & P(R_{t+1}|M_t) & P(R_{t+1}|R_t) \\ \end{bmatrix}` $$ A single entry reads as follows : `\(P(M_{t+1} | T_t)\)` is the probability that a cell that was in state `\(T\)` at time `\(t\)` becomes a cell in state `\(M\)` at time `\(t+1\)`. Ecologically speaking, this would indicate a colonization by a `\(B\)` tree. --- # What we will do - Evaluate the transition matrix (each row is a multinomial model) - Use a grid search to evaluate paramaters - Putting some constraints on parameters to answer ecological questions --- # The problem The multinomial distribution is a generalization of the binomial to `\(n\)` discrete states. It is therefore more complicated to evaluate the different probabilities because the sum for the different states must equal to 1. You can't throw as many values as there are states when you try candidate parameters, you need to find a way that the sum is always 1. --- # Steps 1. Try first the transition from the state `\(B\)` (you will have to subset the data) 2. Write a likelihood function that has 4 parameters (the transition probabilities) 3. Test your function with candidate parameter values 4. Propose alternative candidate parameter values and compare 5. Redo the whole thing with a grid search to find your maximum likelihood estimates --- # Scientific questions 1. Are all transition probabilities `\(> 0\)` ? 2. Compare a model where transitions `\(P(R_{t+1}|B_t)\)`, `\(P(R_{t+1}|M_t)\)` and `\(P(R_{t+1}|T_t)\)` are the same. Is the likelihood lower or higher ? 3. Divide the data in three categories according to annual average temperature and recompute the transition matrices for each section. Then after, compute the likelihood for the entire dataset, with specific models for the different subsections. Is the model better ? --- # Playing with the models Markov chains are commonly used in forest ecology to simulate dynamics over the long term. You can run the model to get at the equilibrium distribution with a simple R code, where `\(\mathbf{T}\)` is a transition matrix and `\(\mathbf{X_t}\)` is a vector of initial conditions (relative abundance of the different states). You run the following : ```r Xt1 <- T %*% Xt ``` Once you are done with evaluation of the global and the local transition matrices, you can run the model for a few time steps (you can use pretty much any initial condition you want, as long as the vector sums to 1) and compare the solutions.