library(readxl)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
plants <- read_excel("individualData_noStages.xlsx")



# ---- USER-DEFINED SIZE CLASS BREAKS ----

# Define the breakpoints for size classes
size_breaks <- c(0, 5, 10, Inf)

# Add sizeClass columns based on the defined breaks
plants <- plants %>%
  mutate(sizeClass_t1 = cut(
    size_t1,
    breaks = size_breaks,
    include.lowest = TRUE,
    labels = FALSE
  )) %>%
  mutate(sizeClass_t2 = cut(
    size_t2,
    breaks = size_breaks,
    include.lowest = TRUE,
    labels = FALSE
  ))


# ---- CALCULATE TRANSITION PROBABILITIES - FROM STAGE 1 ----

# Transition from stage 1-----
#Stage1: total at time 1
stage1Total <- plants %>%
  filter(sizeClass_t1 == 1) %>%
  nrow()

stage1Total
## [1] 296
#Stage1: Survive and remain in stage 1
stage1Remain <- plants %>%
  filter(sizeClass_t1 == 1) %>%
  filter(sizeClass_t2 == 1) %>%
  nrow()

#Calculate the probability of making that transition
(t1_1 <- stage1Remain / stage1Total)
## [1] 0.01351351
#Stage1: Survive and move to stage 2
stage1ToStage2 <- plants %>%
  filter(sizeClass_t1 == 1) %>%
  filter(sizeClass_t2 == 2) %>%
  nrow()

#Calculate the probability of making that transition
(t1_2 <- stage1ToStage2 / stage1Total)
## [1] 0.06081081
#Stage1: Survive and move to stage 3
stage1ToStage3 <- plants %>%
  filter(sizeClass_t1 == 1) %>%
  filter(sizeClass_t2 == 3) %>%
  nrow()

#Calculate the probability of making that transition
(t1_3 <- stage1ToStage3 / stage1Total)
## [1] 0
# ---- CALCULATE TRANSITION PROBABILITIES - FROM STAGE 2 ----
# Stage2: total at time 1
stage2Total <- plants %>%
  filter(sizeClass_t1 == 2) %>%
  nrow()
stage2Total
## [1] 172
# Stage2: Survive and remain in stage 2
stage2Remain <- plants %>%
  filter(sizeClass_t1 == 2) %>%
  filter(sizeClass_t2 == 2) %>%
  nrow()

#Calculate the probability of making that transition
(t2_2 <- stage2Remain / stage2Total)
## [1] 0.2790698
# Stage2: Survive and move to stage 1
stage2ToStage1 <- plants %>%
  filter(sizeClass_t1 == 2) %>%
  filter(sizeClass_t2 == 1) %>%
  nrow()

#Calculate the probability of making that transition
(t2_1 <- stage2ToStage1 / stage2Total)
## [1] 0
# Stage2: Survive and move to stage 3
stage2ToStage3 <- plants %>%
  filter(sizeClass_t1 == 2) %>%
  filter(sizeClass_t2 == 3) %>%
  nrow()

#Calculate the probability of making that transition
(t2_3 <- stage2ToStage3 / stage2Total)
## [1] 0.1569767
# Stage2: Survive and move to stage 1 (regression)
stage2ToStage1 <- plants %>%
  filter(sizeClass_t1 == 2) %>%
  filter(sizeClass_t2 == 1) %>%
  nrow()

#Calculate the probability of making that transition
(t2_1 <- stage2ToStage1 / stage2Total)
## [1] 0
# ---- CALCULATE TRANSITION PROBABILITIES - FROM STAGE 3 ----

stage3Total <- plants %>%
  filter(sizeClass_t1 == 3) %>%
  nrow()
stage3Total
## [1] 32
# Stage3: Survive and remain in stage 3
stage3Remain <- plants %>%
  filter(sizeClass_t1 == 3) %>%
  filter(sizeClass_t2 == 3) %>%
  nrow()
(t3_3 <- stage3Remain / stage3Total)
## [1] 0.96875
# Stage3: Survive and move to stage 2
stage3ToStage2 <- plants %>%
  filter(sizeClass_t1 == 3) %>%
  filter(sizeClass_t2 == 2) %>%
  nrow()
(t3_2 <- stage3ToStage2 / stage3Total)
## [1] 0
# Stage3: Survive and regress to stage 1
stage3ToStage1 <- plants %>%
  filter(sizeClass_t1 == 3) %>%
  filter(sizeClass_t2 == 1) %>%
  nrow()
(t3_1 <- stage3ToStage1 / stage3Total)
## [1] 0
matU <- matrix(rep(0, 9), nrow = 3)

matU <- matrix(
  c(t1_1, t2_1, t3_1, 
    t1_2, t2_2, t3_2, 
    t1_3, t2_3, t3_3),
  nrow = 3,
  byrow = TRUE
)

matU
##            [,1]      [,2]    [,3]
## [1,] 0.01351351 0.0000000 0.00000
## [2,] 0.06081081 0.2790698 0.00000
## [3,] 0.00000000 0.1569767 0.96875
Rage::plot_life_cycle(matU)
# REPRODUCTION (Seedlings) ----
# Approach 1. Distribute seedlings proportionally among stages based on number of seeds

Seedlings_t2 <- 409


# Step 1: Fruit per stage
(
  reprod <- plants %>%
    group_by(sizeClass_t1) %>%
    summarise(nIndivs = n(), totalSeeds = sum(n_seeds)) %>%
    mutate(prop_seeds = totalSeeds / sum(totalSeeds))
)
## # A tibble: 3 × 4
##   sizeClass_t1 nIndivs totalSeeds prop_seeds
##          <int>   <int>      <dbl>      <dbl>
## 1            1     296       1535      0.302
## 2            2     172       2443      0.481
## 3            3      32       1098      0.216
# Step 2: divide up seedlings among the stages, proportionally to fruit produced
reprod <- reprod %>%
  mutate(seedlingsPerStage = Seedlings_t2 * prop_seeds)

reprod
## # A tibble: 3 × 5
##   sizeClass_t1 nIndivs totalSeeds prop_seeds seedlingsPerStage
##          <int>   <int>      <dbl>      <dbl>             <dbl>
## 1            1     296       1535      0.302             124. 
## 2            2     172       2443      0.481             197. 
## 3            3      32       1098      0.216              88.5
# Step 3: convert to a "per individual" measure
reprod  <- reprod %>%
  mutate(seedlingsPerIndiv = seedlingsPerStage / nIndivs)
reprod
## # A tibble: 3 × 6
##   sizeClass_t1 nIndivs totalSeeds prop_seeds seedlingsPerStage seedlingsPerIndiv
##          <int>   <int>      <dbl>      <dbl>             <dbl>             <dbl>
## 1            1     296       1535      0.302             124.              0.418
## 2            2     172       2443      0.481             197.              1.14 
## 3            3      32       1098      0.216              88.5             2.76
#When you have parameterised the survival growth part of the MPM you can create the matrix for that part (matU)

matU <- matrix(
  c(t1_1, t2_1, t3_1, 
    t1_2, t2_2, t3_2, 
    t1_3, t2_3, t3_3),
  nrow = 3,
  byrow = TRUE
)

matU
##            [,1]      [,2]    [,3]
## [1,] 0.01351351 0.0000000 0.00000
## [2,] 0.06081081 0.2790698 0.00000
## [3,] 0.00000000 0.1569767 0.96875
Rage::plot_life_cycle(matU)
#The reproduction matrix can be made by first making a zero matrix and then adding the reproduction to the top row

matF <- matrix(rep(0, 9), nrow = 3)
matF[1, ] <- reprod$seedlingsPerIndiv


matF
##          [,1]    [,2]     [,3]
## [1,] 0.417848 1.14445 2.764738
## [2,] 0.000000 0.00000 0.000000
## [3,] 0.000000 0.00000 0.000000
#THe whole matrix (A) is the sum of the two matrices U and F


matA <- matU + matF
matA
##            [,1]      [,2]     [,3]
## [1,] 0.43136155 1.1444497 2.764738
## [2,] 0.06081081 0.2790698 0.000000
## [3,] 0.00000000 0.1569767 0.968750
Rage::plot_life_cycle(matA)