Step 1 - Load up packages and read in data

# we'll use 'car::Anova' for 'type II' model comparisons
library(car) 
## Loading required package: carData
# load in the data
sim_data <- read.csv(file = "./soay_demog_data.csv")

Step 2 - Fit the vital rate functions

# Reproduction e.g. showing order matters
repr_data <- subset(sim_data, Surv==1 & a>0)
mod_repr <- glm(Repr~a+z, family=binomial, data=repr_data)

Anova(mod_repr) # Both matter
anova(mod_repr, test="LRT") # size matters, age does not. 

mod_repr_a <- glm(Repr~a, family=binomial, data=repr_data)
mod_repr_z <- glm(Repr~z, family=binomial, data=repr_data)
anova(mod_repr, mod_repr_a, test="LRT")
anova(mod_repr, mod_repr_z, test="LRT")

# Now do the real work...

mod_surv <- glm(Surv ~ z + a, family = binomial, data = sim_data)
summary(mod_surv)
Anova(mod_surv)

grow_data <- subset(sim_data, !is.na(z1))
mod_grow <- lm(z1 ~ z + a, data = grow_data)
summary(mod_grow)
Anova(mod_grow)

repr_data <- subset(sim_data, Surv==1 & a>0)
mod_repr <- glm(Repr ~ z + a, family = binomial, data = repr_data)
summary(mod_repr)
Anova(mod_repr)

recr_data <- subset(sim_data, Surv==1 & Repr==1)
mod_recr <- glm(Recr ~ z + a, family = binomial, data = recr_data)
summary(mod_recr)
Anova(mod_recr)

rcsz_data <- subset(sim_data, !is.na(Rcsz))
mod_rcsz <- lm(Rcsz ~ z + a, data=rcsz_data)
summary(mod_rcsz)
Anova(mod_rcsz)

# refit the recruitment model, dropping size
mod_recr <- glm(Recr ~ a, family = binomial, data = recr_data)

Step 3 - Store the parameters

m_par <- c(
  # survival
  surv      = coef(mod_surv),
  # growth 
  grow      =  coef(mod_grow),
  grow_sd   =  summary(mod_grow)$sigma,
  # reproduce or not
  repr      =  coef(mod_repr),
  # recruit or not
  recr      =  coef(mod_recr),
  # recruit size
  rcsz      =  coef(mod_rcsz),
  rcsz_sd   =  summary(mod_rcsz)$sigma)

names(m_par) <- 
  sub(pattern = "(Intercept)", replace = "int", 
      x = names(m_par), fixed = TRUE)

names(m_par) <- 
  sub(pattern = ".", replace = "_", 
      x = names(m_par), fixed = TRUE)

Step 4 - Define the IPM vital rate functions

# Growth function, given you are size z and age a now returns the 
# pdf of size z1 next time
g_z1z <- function(z1, z, a, m_par){
  mean <- m_par["grow_int"] + 
    m_par["grow_z"] * z + m_par["grow_a"] * a
  sd <- m_par["grow_sd"]
  p_den_grow <- dnorm(z1, mean = mean, sd = sd)
  return(p_den_grow)
}

# Survival function, logistic regression
s_z <- function(z, a, m_par){
  linear_p <- m_par["surv_int"] + 
    m_par["surv_z"] * z + m_par["surv_a"] * a 
  p <- 1/(1+exp(-linear_p))                           
  return(p)
}

# Reproduction function, logistic regression
pb_z <- function(z, a, m_par){
  if (a==0) {
    p <- 0
  } else {
    linear_p <- m_par["repr_int"] + 
      m_par["repr_z"] * z + m_par["repr_a"] * a
    p <- 1/(1+exp(-linear_p))
  }
  return(p)
}

# Recruitment function, logistic regression
pr_z <- function(a, m_par) {
  linear_p <- m_par["recr_int"] + m_par["recr_a"] * a
  p <- 1/(1+exp(-linear_p))
  return(p)
}

# Recruit size function, given you are size z and age a now returns 
# the pdf of size z1 recruits next summer
c_z1z <- function(z1, z, a, m_par){
  mean <- m_par["rcsz_int"] + 
    m_par["rcsz_a"] * a + m_par["rcsz_z"] * z
  sd <- m_par["rcsz_sd"]                          
  p_den_rcsz <- dnorm(z1, mean = mean, sd = sd)   
  return(p_den_rcsz)
}

Step 5 - Define the IPM component kernel functions

# Define the survival-growth kernel
P_z1z <- function (z1, z, a, m_par) {
  return( s_z(z, a, m_par) * g_z1z(z1, z, a, m_par) )
}

# Define the reproduction kernel
F_z1z <- function (z1, z, a, m_par) {
  return( s_z(z, a, m_par) * pb_z(z, a, m_par) * 
            (1/2) * pr_z(a, m_par) * c_z1z(z1, z, a, m_par) )
}

Define functions to numerically implement the IPM

# Calculate the mesh points, mesh width and store with upper/lower 
# bounds and max age
mk_intpar <- function(m, L, U, M) {
  h <- (U - L) / m
  meshpts  <-  L + ((1:m) - 1/2) * h
  na <- M + 2
  return( list(meshpts = meshpts, M = M, na = na, h = h, m = m) )
}


# make initial state distribution (uniform size, age 0 only)
mk_init_nt0 <- function(i_par) {
  nt <- lapply(seq_len(i_par$na), function(x) rep(0, i_par$m))
  nt[[1]] <- rep(1, i_par$m) / (i_par$h * i_par$m)
  return(nt)
}

# Build the list of age/process specific kernels + store the 
# integration parameters in the same list
mk_age_IPM <- function(i_par, m_par) {
  within(i_par, {
    F <- P <- list()
    for (ia in seq_len(na)) {
      F[[ia]] <- outer(meshpts, meshpts, F_z1z, a = ia-1, m_par = m_par) * h
      P[[ia]] <- outer(meshpts, meshpts, P_z1z, a = ia-1, m_par = m_par) * h
    }
    rm(ia)
  })
}

Step 7 - Define functions to iterate the IPM

# 'right' iterate one time step to project the pop dynamics. 
# 'x' is the list of age-specific size dists
r_iter <- function(x, F, P) {
  na <- length(x)
  xnew <- list(0)
  for (ia in seq_len(na)) {
    xnew[[1]] <- (xnew[[1]] + F[[ia]] %*% x[[ia]])[,,drop=TRUE]
  }
  for (ia in seq_len(na-1)) {
    xnew[[ia+1]] <- (P[[ia]] %*% x[[ia]])[,,drop=TRUE]
  }
  xnew[[na]] <- xnew[[na]] + (P[[na]] %*% x[[na]])[,,drop=TRUE]
  return(xnew)
}

# 'left' iteration one time step to project reproductive value. 
# 'x' is the list of age-specific size dists
l_iter <- function(x, F, P) {
  na <- length(x)
  xnew <- list(0)
  for (ia in seq_len(na-1)) {
    xnew[[ia]] <- (x[[ia+1]] %*% P[[ia]]  + x[[1]] %*% F[[ia]])[,,drop=TRUE]
  }
  xnew[[na]] <- (x[[na]] %*% P[[na]] + x[[1]] %*% F[[na]])[,,drop=TRUE]
  return(xnew)
}