Direct effects of a non-native invader erode native plant fitness in the forest understory

Authors | Lalasia Bialic-Murphy^1 , Nathan Brouwer^2 , Susan Kalisz^1

1 Department of Ecology and Evolutionary Biology, University of Tennessee Knoxville, Knoxville, Tennessee, USA. 2 Department of Conservation and Field Research, The National Aviary, Pittsburgh, Pennsylvania, USA.

Code written by | Lalasia Bialic-Murphy, 2019


Overview | This file uses parameter coefficients calculated in a separate R script (not provided in our supplement) and builds annual integral projection models from 2010-2016 for four treatments: A) deer access and garlic mustard ambient, B) deer access and garlic mustard weeded, C) deer exclusion and garlic mustard ambient, and D) deer exclusion and garlic mustard weeded (i.e., removal of both stressors). Using the component matrices of the discretized kernel matrix \(K=P+F\), this script generates the annual population growth rates, the cumulative growth rates from 2010 to 2016 (λc), and the time-averaged geometric mean growth rate, which caputures the time-averaged population growth rate over a single transition year (λper year = sixth root of λc).

This document executes R code and was composed using R Markdown. The required demographic parameter estimates can be uploaded from the associated Rda files in the archived Dryad Digital Repository http://doi:10.5061/dryad.4k7d5n7 (Bialic-Murphy et al. 2019).

R Packages | Required libraries containing necessary functions for this script.

library(popbio)      # for matrix analysis functions

Load the working directory where the Rda data files are stored

rm(list=ls(all=T))
setwd("~/")

I. Vital rate functions and integral projection model Load the vital rate coefficients contained in the P.array(year_i).Rda files. The P.arrays contain 19 coefficients for treatments A (column 1,named D G), B (column 2, named ND G), C (column 3, named D NG), and D (column 4, named ND NG). The first P.array in the associated annual Rda file contains the coefficients based on the raw data and the following 2000 p.arrays are based on the permutation analysis to test for significant differences between the garlic mustard ambient treatment and the garlic mustard weeded treatment (Caswell 2001), where individuals within each plot were randomly assigned to either the garlic mustard ambient treatment or the garlic mustard weeded treatment.

Note: Only the first two p.arrays for each treatment and each year are loaded here to make the code run faster. To see the full output, including the results for the permuation test, change nran to 2000.

load('Data_Dryad/P.array2010.Rda')
p.2010 <- p.array
p.2010 <- p.array[1:19, 1:4, 1:2]

load('Data_Dryad/P.array2011.Rda')
p.2011 <- p.array 
p.2011 <- p.array[1:19, 1:4, 1:2]

load('Data_Dryad/P.array2012.Rda')
p.2012 <- p.array
p.2012 <- p.array[1:19, 1:4, 1:2]

load('Data_Dryad/P.array2013.Rda')
p.2013 <- p.array
p.2013 <- p.array[1:19, 1:4, 1:2]

load('Data_Dryad/P.array2014.Rda')
p.2014 <- p.array
p.2014 <- p.array[1:19, 1:4, 1:2]

load('Data_Dryad/P.array2015.Rda')
p.2015 <- p.array
p.2015 <- p.array[1:19, 1:4, 1:2]
str(p.2015)

Compute meshpoints to match the discretized kernel matrix (\(K\)).

s0 <- c(0.16408, 4.41405) # Size range
bigM <- 150 # Discretized kernel dimensions
min.sz<-0.9*min((range(s0,na.rm=T)))
max.sz<-1.1*max((range(s0,na.rm=T)))
h=(max.sz-min.sz)/(bigM+1); 
x=(seq(min.sz, max.sz, length=bigM) + seq(min.sz+h, max.sz+h, length=bigM))/2; 

nran <- 2 # Input the number of random samples
N_stressors <- 4

Set-up Kernel components and compile the kernel \(K\). The pvec = the vital rate coefficents 1-19 and the stressor.n = the associated treatments (A,B,C, and D)

1 The survival-growth function \(p(y,x)\)

pyx<-function(y,x, pvec, stressor.n) { 
  p<-sx(x, pvec, stressor.n)*gyx(y, x, pvec, stressor.n)
  return(p) 
}

1.1 Survival probability function \(s(x)\)

sx<-function(x, pvec, stressor.n) {
  xbeta<-pvec[1,stressor.n] + pvec[2,stressor.n]*x; ## pr of survival
  s<-exp(xbeta)/(1+exp(xbeta))
  return(s);
}

1.2 Growth function \(g(y,x)\)

gyx<-function(y, x, pvec, stressor.n) {
  mux<-pvec[3, stressor.n] + pvec[4, stressor.n]*x ## Change in size from time t to t+1
  sigmax2<-pvec[5, stressor.n]*exp(2*pvec[6, stressor.n]*mux) 
  ## sigmax2 accounts for the size dependent variance of growth, weight=varExp
  sigmax<-sqrt(sigmax2)
  g<-dnorm(y, mux, sigmax)
  return(g)
}

2 Fertility function \(f(y,x)\), The fertility function is composed of a series of regression models and three discrete life-stage transitions. These life-stage transitions are represented as a series of equations. Refer to Appendix S1 of the associated JoE manuscript (Lalasia Bialic-Murphy,2019) for a description of the transitions represented by the mathematical symbols annotated in the below R code ‘chunk’.

fyx<-function(x = seq(0.1,4,length.out = 150), 
              pvec, 
              stressor.n, 
              second.fruit = TRUE,
              print.betas = FALSE,
              plot.fec = FALSE
) {
  
  if(print.betas == TRUE){
    print(pvec[c(7,8,9,10,11,12,19,20), stressor.n])
  }
  a <-exp(pvec[7, stressor.n]  + pvec[8,  stressor.n]*x)/(1+exp(pvec[7,  stressor.n] + pvec[8,  stressor.n]*x)) 
      # a = pr of producing 1 fruit, (f_1)  
  b <-exp(pvec[9, stressor.n]  + pvec[10, stressor.n]*x)/(1+exp(pvec[9,  stressor.n] + pvec[10, stressor.n]*x)) 
      # b = pr of producing a 2nd fruit, conditional on flwr'ing (f_2)
  b[is.nan(b)] = 0 # Change NaN to zero
  c <-pvec[19, stressor.n]  # Herbivory of flowering plant (x,θ_D)
  d <-pvec[11, stressor.n]  # C. Mean number of seeds produced per fruit (c)
  e <-pvec[12, stressor.n]  # D. Mean germination rate, γ_(s-g)
  if(second.fruit == FALSE){
    f <- a*c*d*e
  }
  # Add 2nd fruiting stem
  if(second.fruit == TRUE){
    f <- a*c*d*e + a*b*c*d*e
  }
  return(f)
}

3 Transitions from the discrete life stages to the continuous size-dependent component of the IPM kernel

leaf_Pyx<-function(y, pvec, stressor.n) {
  a <- dnorm(y, pvec[17,stressor.n], pvec[18,stressor.n]) 
       # a = mean size and SD of 3 leaf stage,η_t (y)=N(μ_t,σ_t)  
  b <- pvec[16, stressor.n] 
       # b = mean single-leaf to 3 leaf stage transition, γ_(o-t) 
  c <- a*b
  return(c)
}

4 The kernel: \(K(y,x)= p(y,x) + f(y,x) + c(y,x)\)

Kyx<-function(y, x, pvec, stressor.n) {
  k<-pyx(y, x, pvec, stressor.n)+fyx(x, pvec, stressor.n,second.fruit = TRUE)
  return(k) 
}

5 Numerical integration function

This function can be used to compute a discretized matrix either for the kernel or its component parts: the survival-growth matrix \(P\), the fertility matrix \(F\).

bigmat<-function(bigM, pvec, stressor.n){
  # Set matrix size and convergence tolerance 
  min.sz<-min.sz
  max.sz<-min.sz
  
  # Compute meshpoints iteration matrix KD 
  h=(max.sz-min.sz)/(bigM+1); 
  y=(seq(min.sz,max.sz,length=bigM)+seq(min.sz+h,max.sz+h,length=bigM))/2;  
  
  K=outer(y,y,Kyx, pvec, stressor.n);
  KD=h*K;
  return(KD);        
}

6 Add discrete stages

add.discrete.stages <- function(stressor.n,pvec,bigM) {
  dim_mat = bigM+3
  min.sz<-min.sz # Set matrix size and convergence tolerance
  max.sz<-max.sz
 
  h=(max.sz-min.sz)/(bigM+1); # Compute meshpoints iteration matrix KD
  y=(seq(min.sz, max.sz, length=bigM) + seq(min.sz+h, max.sz+h, length=bigM))/2;  
  
  pyx.i <- h*outer(y,y,  # h = step size; y = mesh
                   pyx,  # Growth function
                   pvec, # Betas going into growth function g.yx
                   stressor.n)
  
  fyx.i <- fyx(y, pvec, stressor.n) 

  # Creating matrix so that we can add discrete stages (i.e., germ-coty, seedling-1leaf, 1 leaf-1 leaf)
  full_mat <- matrix(0, nrow = dim_mat,ncol = dim_mat)
  # Add p matrix to dummy matrix
  full_mat[-c(1:3),-c(1:3)] <- pyx.i  
  # Add fec vect dummy matrix
  full_mat[1,-c(1:3)] <- fyx.i
  
  # Add 3L recruit vector to size dist
  R.vector <- leaf_Pyx(y,pvec,stressor.n)
  full_mat[-c(1:3),3] <- R.vector
  
  # Add additional transitions
  full_mat[2,1] <- pvec[13] # Germ to cotyledon stage (germ to cot), γ_(g-c)
  full_mat[3,2] <- pvec[14] # Cotyledon to single leaf stage, γ_(c-o) 
  full_mat[3,3] <- pvec[15] # Statsis of 1 leaf stage, 1L to 1L (O to O), γ_(o-o)
  
  output <-list(full.mat = full_mat,
                pyx = pyx.i,
                fxy = fyx.i
  )
  return(output)
}

II. Evaluate the integral projection model

Combine all kernel components and calcualte annual population growth rate for each treatment and the cumulative growth rate, λc, from 2010 to 2016

s1 <- add.discrete.stages(stressor.n = 1,p.array[,,1],bigM)
s2 <- add.discrete.stages(stressor.n = 2,p.array[,,1],bigM)
s3 <- add.discrete.stages(stressor.n = 3,p.array[,,1],bigM)
s4 <- add.discrete.stages(stressor.n = 4,p.array[,,1],bigM)

lambdas <- c(eigen.analysis(s1$full.mat)$lambda1
             ,eigen.analysis(s2$full.mat)$lambda1
             ,eigen.analysis(s3$full.mat)$lambda1
             ,eigen.analysis(s4$full.mat)$lambda1)

# For the cumulative growth rate, we multiple the sequence of matrices. 
eigen.analysis(s4$full.mat%*%s3$full.mat%*%s2$full.mat%*%s1$full.mat)$lambda1
eigen.analysis(s1$full.mat%*%s2$full.mat%*%s3$full.mat%*%s4$full.mat)$lambda1
lambdas
lambdas[1]*lambdas[2]*lambdas[3]*lambdas[4]
# Lambda 1 = treatment A, lambda 2 = treatment B, lambda 2 = treatment C, lambda 4 = treatment D

# Create empty arrays to receive output
lam.2010 <- lam.2011 <- lam.2012 <- lam.2013 <- lam.2014 <- lam.2015 <- cum.lams <- matrix(NA,nrow=nran,ncol=N_stressors)
dimnames(cum.lams) <- list(dimnames(p.2011)[[3]],dimnames(p.2011)[[2]])
dimnames(lam.2010) <-dimnames(lam.2011) <- dimnames(lam.2012) <- dimnames(lam.2013) <- dimnames(lam.2014) <-dimnames(lam.2015) <-dimnames(cum.lams)

A.2010 <- A.2011 <- A.2012 <- A.2013 <- A.2014 <- A.2014 <-A.2015 <- A.cum <- array(NA,dim=c(bigM+3,bigM+3,nran,N_stressors))

start.time <- Sys.time()
for (i in 1:nran) { # Create matrix and lambda for each random permutation
  for (j in 1:N_stressors) {
    A.2010[,,i,j] <- add.discrete.stages(stressor.n = j,p.2010[,,i],bigM)$full.mat
    lam.2010[i,j] <- eigen.analysis(A.2010[,,i,j])$lambda1
    A.2011[,,i,j] <- add.discrete.stages(stressor.n = j,p.2011[,,i],bigM)$full.mat
    lam.2011[i,j] <- eigen.analysis(A.2011[,,i,j])$lambda1
    A.2012[,,i,j] <- add.discrete.stages(stressor.n = j,p.2012[,,i],bigM)$full.mat
    lam.2012[i,j] <- eigen.analysis(A.2012[,,i,j])$lambda1
    A.2013[,,i,j] <- add.discrete.stages(stressor.n = j,p.2013[,,i],bigM)$full.mat
    lam.2013[i,j] <- eigen.analysis(A.2013[,,i,j])$lambda1
    A.2014[,,i,j] <- add.discrete.stages(stressor.n = j,p.2014[,,i],bigM)$full.mat
    lam.2014[i,j] <- eigen.analysis(A.2014[,,i,j])$lambda1
    A.2015[,,i,j] <- add.discrete.stages(stressor.n = j,p.2015[,,i],bigM)$full.mat
    lam.2015[i,j] <- eigen.analysis(A.2015[,,i,j])$lambda1
    A.cum[,,i,j]  <- A.2015[,,i,j] %*% A.2014[,,i,j] %*%A.2013[,,i,j] %*% A.2012[,,i,j] %*% 
                     A.2011[,,i,j]%*% A.2010[,,i,j]
    cum.lams[i,j] <- eigen.analysis(A.cum[,,i,j])$lambda1
  }
}
end.time <- Sys.time(); end.time - start.time