In brouwern/popcycles: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(popcycles)
library(magic)

Perigine falcon demography as vec-perm (Wooten and Bell 1992)

Dispersal

Dispersal is organized by the stage of the birds: the matrix M_1 will describe mi8gration for juveniles, while M_2 will describe movements of adults

(TODO: re-organize dispersal by patch? M.j and M.a? does this work?)

The dispersal matrix for juveniles (M1) is $$\mathbf{M_1} = \left[\begin{array} {rrr} 1-d & d \ d & 1-d \ \end{array}\right] $$ The probablity of moving is d, while the probablity of remianing at the site where they were born is 1-d. The probablity of dispersal is 0.27, so the probablity of remaining within the original site is 0.73.

We'll assign the dispersal probablity to the object d.

d  <- 0.27 # dispersal

To make the matrix. In R we can make this matrix like this.

M1 <- matrix(data = c(1-d, d,
                      d,   1-d),
             byrow= T,
             nrow = 2)

Adults do not disperse, so the dispersal matrix for adults (M2) is just an identity matrix

$$\mathbf{M_1} = \mathbf{I} = \left[\begin{array} {rrr} 1 & 0 \ 0 & 1 \ \end{array}\right] $$

M2 <- diag(2)

We can put these

M <- magic::adiag(M1,M2)

Demography

Demography within each geographic location is described by three parametrs

juvenile survival = p_i
adult survival = q_i
adult fecundity = f_i

Demography is described by a 2 x 2 matrix

$$\mathbf{B_i} = \left[\begin{array} {rrr} 0 & f_i \ p_i & q_i \ \end{array}\right] $$ The matrix is indexed B_i because there will be two matrices (B1, B2), one for each geographic region.

Combining the demographic matrices B1 and B2 into a block matrix gives you

$$\mathbf{B} = \left[\begin{array}{c|c} B_1 & 0 \ \hline 0 & B_2 \ \end{array}\right] $$

This block matrix represents both demographic matrices. We could write out the whole matrix like this

$$\mathbf{B} = \left[\begin{array}{cc|cc} 0 & f_1 & 0 & 0 \ p_1 & q_1 & 0 & 0 \ \hline 0 & 0 & 0 & f_2 \ 0 & 0 & p_2 & q_2 \ \end{array}\right] $$

Model parameters

Parameters for the model are taken from Wooten & x 199x

# Parameters - site 1
f1 <- 0.26
p1 <- 0.72
q1 <- 0.77

# Parameters - site 2
f2 <- 0.19
p2 <- 0.72
q2 <- 0.77

Demography matrices

B1 <- matrix(data = c(0,f1,p1,q1),
             byrow= T,
             nrow = 2)

B2 <- matrix(data = c(0,f2,p2,q2),
             byrow= T,
             nrow = 2)

Build block demography matrix

B <- magic::adiag(B1,B2)

Look at the block matrix

Vec-permutation matrix as given in paper

P <- matrix(data = c(1,0,0,0,
                     0,0,1,0,
                     0,1,0,0,
                     0,0,0,1),
            byrow = T,
            nrow =4)

The p-matrix (permutation matrix)

Eigen analysis of demographic matrix

The eigen values and all eigen vectors are given by the eigen() command.

The dominant eigenvalue is the population growth rate.

We can look at each subpopulation like this

This subpopulation has the highest growth rate

eigen(B[1:2,1:2]) # 0.9641589

This subpopulation has the lower growth rate

eigen(B[3:4,3:4]) # 0.9188773

(TODO: are these values correct?)

We do eigenanalysis on the full matrix, but I'm not sure exactly if/how to interpret this. Without dispersal connecting the subpopulations this may violate some key rules of matrix population models.

eigen(B)

Eigen vectors are used in sensitivity and elasticity analysis. I forget if the matrix of results is organized by row or column

TODO: I have a function called calc_lam() that returns just the biologically important information, need to find where I put it / re-make it

Eigen analysis of fulle system

Bt(P)M*P

The eigen values and all eigen vectors are given by the eigen() command.

The dominant eigenvalue is the population growth rate (TODO: is the value correct?)

Eigen vectors are used in sensitivity and elasticity analysis, but for models using vec-permutation I believe calculation of sensitivity is different

eigen(B%*%t(P)%*%M%*%P)

Build the vec-permutation by hand one matrix at a time

# First permutation marix
## matrix
E11 <- matrix(data = c(1,0,
                     0,0),
            byrow = T,
            nrow =2)

#
E1 <- kronecker(E11,t(E11))

E12 <- matrix(data = c(0,1,
                       0,0),
              byrow = T,
              nrow =2)

E2 <- kronecker(E12,t(E12))


E21 <- matrix(data = c(0,0,
                       1,0),
              byrow = T,
              nrow =2)

E3 <- kronecker(E21,t(E21))

E22 <- matrix(data = c(0,0,
                       0,1),
              byrow = T,
              nrow =2)

E4 <-kronecker(E22,t(E22))

E1+E2+E3+E4

Black-headed gull Lebreton 1996

B.prad1 <- matrix(data = c(0.00, 0.096, 0.160, 0.224, 0.320,
                          0.80, 0,     0,     0,     0,
                          0.00, 0.820, 0,     0,     0,
                          0.00, 0.000, 0.82,  0,     0,
                          0.00, 0.000, 0.00,  0.82,  0.82),
                 nrow = 5, byrow = T)
B.prad2 <- matrix(data = c(0.00, 0.100, 0.160, 0.200, 0.200,
                           0.80, 0,     0,     0,     0,
                           0.00, 0.820, 0,     0,     0,
                           0.00, 0.000, 0.82,  0,     0,
                           0.00, 0.000, 0.00,  0.82,  0.82),
                  nrow = 5, byrow = T)


eigen(B.prad2)
#calc.lam(B.prad2)

B.prad.diag <- adiag(B.prad1, B.prad2)


### Pradel migratio matries

#only 1st stage dispersales
M1 <- matrix(data = c(0.75, 0.375,
                      0.25, 0.625),
             nrow = 2, byrow =T)

M5 <- M4 <- M3 <- M2 <- diag(nrow = 2)

M.pradel.diag <- adiag(M1,M2,M3,M4,M5)



### Create permuation matrix
s <- dim(B.prad1)[1]
p <- dim(M1)[1]

E.blank.pradel <- matrix(data =0, 
                         ncol = p,nrow = s)

P.pradel <- kronecker(E.blank.pradel,t(E.blank.pradel))
for(i in 1:dim(E.blank.pradel)[1]){

  for(j in 1:dim(E.blank.pradel)[2]){
    E.blank.work <- E.blank.pradel

    E.blank.work[i,j] <- 1

    E.blank.work <- kronecker(E.blank.work,t(E.blank.work))
    P.pradel <- P.pradel +E.blank.work
  }

}



### create final matrix


P.pradel%*%B.prad.diag%*%t(P.pradel)%*%M.pradel.diag



















E.blank <- matrix(data =0, ncol = 2,nrow = 2)

P <- kronecker(E.blank,t(E.blank))
for(i in 1:dim(E.blank)[1]){


  for(j in 1:dim(E.blank)[2]){
    E.blank.work <- E.blank

    E.blank.work[i,j] <- 1

    E.blank.work <- kronecker(E.blank.work,t(E.blank.work))
    P <- P +E.blank.work
  }

}