IPMpack Example
Goals
In this lab, you will go through the steps of using an integral projection model and get acquainted with the R package IPMpack. The point of this experience is to become more conversant and comfortable with data-based modeling in general, and with structured population models in particular.
IPMpack is an R package designed to help you learn about IPM’s in general, and also for suphisticated population modeling.
This work comes from the appendix of Merow et al. (2014) Advancing population ecology with integral projection models: a practical guide. Methods in Ecology and Evolution 5:99-110.
Study organism
“The Austrian Dragonhead (Dracocephalum austriacum L., Lamiaceae) is a long-lived perennial plant occurring from the Pyrenees to the Caucasus, with scattered populations in Western Europe (Bensettiti et al. 2002). The plant is 20 to 50 cm high, with hairy stems and large blue-violet insect-pollinated flowers (3.5 to 5 cm in length). Dracocephalum austriacum does not spread clonally and seeds lack adaptations for dispersal. It usually occurs in habitat patches where competition with other plants is low and competition with shrubs and trees may have a negative effect on performance (Olivier et al. 1995). Dracocephalum austriacum is listed as Vulnerable by the World Conservation Union (IUCN) and is protected under National, European and international legislation (Red List of French Endangered Flora, Habitats Directive, Bern Convention). Threats include pillaging, trampling and damage by grazing cattle (Bensettiti et al. 2002). Previous studies with this species have shown large population differences in adaptive genetic variation (Bonin et al. 2007) and positive correlations between population heterozygosity and stochastic population growth rate (Nicole 2005). The present field study was carried out in 7 of the 15 known French populations of D. austriacum. The study populations occur in the subalpine zone (1250 to 2000 m a.s.l.) of the Alps, on exposed, xeric, stony grasslands or heaths, preferentially on calcareous, thin soils (Lauber and Wagner 1998). Populations vary widely in size, density, habitat characteristics and soil composition.”
Dracocephalum
Let’s get ready. Load the package.
library(IPMpack)## Loading required package: Matrix## Loading required package: MASS## Loading required package: nlme## Load data
d <- read.csv("Appendix_A_Simple_data.csv")
## examine data
head(d)## size sizeNext surv fec.seed fec.flower
## 1 3.09 NA 0 NA NA
## 2 2.81 NA 0 NA NA
## 3 4.46 6.07 1 15 1
## 4 1.68 NA 0 NA NA
## 5 3.99 NA 0 NA NA
## 6 4.07 NA 0 NA NA## Make graphs of important demographic relationships.
par(mfrow=c(2,2),mar=c(4,4,2,1))
plot(d$size,jitter(d$surv),xlab="Size (t)", ylab="Survival to t+1") # jittered
plot(d$size,d$sizeNext,xlab="Size (t)",ylab="Size (t+1)")
plot(d$size,d$fec.seed,xlab="Size (t)",ylab="Seed Number")
hist(d$sizeNext[is.na(d$size)],main=,xlab="Recruit Size")
## make placeholder to store results from demographic relationships
params=data.frame(
surv.int=NA, # Intercept from logistic regression of survival
surv.slope=NA, # Slope from logistic regression of survival
growth.int=NA,
growth.slope=NA,
growth.sd=NA,
seed.int=NA,
seed.slope=NA,
recruit.size.mean=NA, # Mean recruit size
recruit.size.sd=NA, # Standard deviation of recruit size
establishment.prob=NA # Probability of establishment
# Intercept from linear regression of growth
# Slope from linear regression of growth
# Residual sd from the linear regression of growth
# Intercept from Poisson regression of seed number
# Slope from Poisson regression of seed number
)
#################################################
## Fit statistical models to demeographic relationships
##################################################
# 1. survival: logistic regression
surv.reg=glm(surv~size,data=d,family=binomial())
summary(surv.reg)##
## Call:
## glm(formula = surv ~ size, family = binomial(), data = d)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2275 -0.5815 0.2491 0.5597 2.0972
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.9646 0.4760 -8.329 <2e-16 ***
## size 1.2897 0.1338 9.642 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 488.69 on 399 degrees of freedom
## Residual deviance: 315.82 on 398 degrees of freedom
## (100 observations deleted due to missingness)
## AIC: 319.82
##
## Number of Fisher Scoring iterations: 5params$surv.int=coefficients(surv.reg)[1]
params$surv.slope=coefficients(surv.reg)[2]
# 2. growth: linear regression
growth.reg=lm(sizeNext~size,data=d)
summary(growth.reg)##
## Call:
## lm(formula = sizeNext ~ size, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.8037 -0.5434 0.0932 0.5741 2.6732
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.68135 0.19580 13.69 <2e-16 ***
## size 0.57922 0.03902 14.84 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8866 on 278 degrees of freedom
## (220 observations deleted due to missingness)
## Multiple R-squared: 0.4421, Adjusted R-squared: 0.4401
## F-statistic: 220.3 on 1 and 278 DF, p-value: < 2.2e-16params$growth.int=coefficients(growth.reg)[1]
params$growth.slope=coefficients(growth.reg)[2]
params$growth.sd=sd(resid(growth.reg))
# 3. seeds: Poisson regression
seed.reg=glm(fec.seed~size,data=d,family=poisson())
summary(seed.reg)##
## Call:
## glm(formula = fec.seed ~ size, family = poisson(), data = d)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6.5338 -1.8830 0.0529 1.5375 4.9690
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.33234 0.06370 20.92 <2e-16 ***
## size 0.30647 0.01164 26.34 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2598.3 on 279 degrees of freedom
## Residual deviance: 1840.1 on 278 degrees of freedom
## (220 observations deleted due to missingness)
## AIC: 2942.2
##
## Number of Fisher Scoring iterations: 5params$seed.int=coefficients(seed.reg)[1]
params$seed.slope=coefficients(seed.reg)[2]
# 4. size distribution of recruits
params$recruit.size.mean=mean(d$sizeNext[is.na(d$size)])
params$recruit.size.sd=sd(d$sizeNext[is.na(d$size)])
# 5. establishment probability
params$establishment.prob=sum(is.na(d$size))/sum(d$fec.seed,na.rm=TRUE)
## show how the statistical models conform to the data
par(mfrow=c(2,2))
xx=seq(0,8,by=.01) # sizes at which to evaluate predictions
plot(d$size,jitter(d$surv), xlab="Size (t)",ylab="Survival to t+1") # jittered
lines(xx,predict(surv.reg,
data.frame(size=xx),type="response"), col="red",lwd=3)
plot(d$size,d$sizeNext,xlab="Size (t)",ylab="Size (t+1)")
lines(xx,predict(growth.reg,data.frame(size=xx)),col="red",lwd=3)
plot(d$size,d$fec.seed,xlab="Size (t)",ylab="Seed Number (t)")
lines(xx,predict(seed.reg,
data.frame(size=xx),type='response'),col='red',lwd=3)
hist(d$sizeNext[is.na(d$size)],freq=FALSE,xlab="Recruit size",main='')
lines(xx,dnorm(xx,params$recruit.size.mean,
params$recruit.size.sd), col='red',lwd=3)
#######################################
## calculate demographic rates for the IPM using coefficients
## from the statistically fitted relationships
## and create the IPM Kernel
#########################################
# 1. survival probability function
s.x=function(x,params) {
u=exp(params$surv.int+params$surv.slope*x)
return(u/(1+u))
}
# 2. growth function
g.yx=function(xp,x,params) {
dnorm(xp,mean=params$growth.int+params$growth.slope*x,sd=params$growth.sd) }
# 3. reproduction function
f.yx=function(xp,x,params) {
params$establishment.prob*
dnorm(xp,mean=params$recruit.size.mean,sd=params$recruit.size.sd)*
exp(params$seed.int+params$seed.slope*x)
}
## create kernel
## create the details of the kernel matrix,
## including the size (100 x 100 matrix)
min.size=.9*min(c(d$size,d$sizeNext),na.rm=T)
max.size=1.1*max(c(d$size,d$sizeNext),na.rm=T)
n=100 # number of cells in the matrix
b=min.size+c(0:n)*(max.size-min.size)/n # boundary points
y=0.5*(b[1:n]+b[2:(n+1)]) # mesh points
h=y[2]-y[1] # step size
## Create the kernel
G=h*outer(y,y,g.yx,params=params) # growth matrix
S=s.x(y,params=params) # survival
F=h*outer(y,y,f.yx,params=params) # reproduction matrix
P=G # placeholder; redefine P on the next line
for(i in 1:n) P[,i]=G[,i]*S[i] # growth/survival matrix
K=P+F # full matrix
########################################
## Analyze the kernel
#######################################
# Lambda
(lam <- Re(eigen(K)$values[1]))## [1] 1.013391# Is the population likely to grow or shrink?
## Later on, we plot the following:
# Stable stage distribution
w.eigen <- Re(eigen(K)$vectors[,1])
stable.dist <- w.eigen/sum(w.eigen)
## Reproductive value
v.eigen <- Re(eigen(t(K))$vectors[,1])
repro.val <- v.eigen/v.eigen[1]
## The eigenvectors and lambda can be combined to obtain
## the sensitivity and elasticity matrices.
v.dot.w=sum(stable.dist*repro.val)*h
sens=outer(repro.val,stable.dist)/v.dot.w
elas=matrix(as.vector(sens)*as.vector(K)/lam,nrow=n)
## Plot summaries
library(fields)## Loading required package: spam## Loading required package: dotCall64## Loading required package: grid## Spam version 2.1-1 (2017-07-02) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.##
## Attaching package: 'spam'## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve## Loading required package: mapspar(mfrow=c(2,3),mar=c(4,5,2,2))
image.plot(y,y,t(K), xlab="Size (t)",ylab="Size (t+1)",
col=topo.colors(100), main="IPM matrix")
contour(y,y,t(K), add = TRUE, drawlabels = TRUE)
plot(y,stable.dist,xlab="Size",type="l",main="Stable size distribution")
plot(y,repro.val,xlab="Size",type="l",main="Reproductive values")
image.plot(y,y,t(elas),xlab="Size (t)",ylab="Size (t+1)",main="Elasticity")
image.plot(y,y,t(sens),xlab="Size (t)",ylab="Size (t+1)", main="Sensitivity")
####################################################
####################################################
####################################################
## EXTRA
## MOre accurately modeling Growth variability
#Size Dependent Variance in Growth
# Plot of the noise in the growth model.
plot(growth.reg$model$size,abs(resid(growth.reg)),
xlab='size', ylab='|residual|')
# in this regression we assume that the variation around the regression line
# increases exponentially with size.
growth.reg=gls(sizeNext~size, weights=varExp(),
na.action=na.omit, data=d)
summary(growth.reg)## Generalized least squares fit by REML
## Model: sizeNext ~ size
## Data: d
## AIC BIC logLik
## 711.298 725.8085 -351.649
##
## Variance function:
## Structure: Exponential of variance covariate
## Formula: ~fitted(.)
## Parameter estimates:
## expon
## 0.2935594
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 2.3280286 0.14227273 16.36314 0
## size 0.6585909 0.03301602 19.94762 0
##
## Correlation:
## (Intr)
## size -0.945
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.29676786 -0.68003669 0.07382739 0.69180794 3.35550523
##
## Residual standard error: 0.1664024
## Degrees of freedom: 280 total; 278 residualplot(d$size,d$sizeNext, main='Growth/Shrinkage/Stasis')
lines(xx,predict(growth.reg, data.frame(size=xx)),col='red',lwd=3)
# (B) Modify the params data frame (where we store all the coefficients
# used to build the IPM) to have coefficients called growth.sd.int and
# growth.sd.slope and set the values of those coefficients equal to the
# appropriate coefficients from part (A).
params$growth.int=coefficients(growth.reg)[1]
params$growth.slope=coefficients(growth.reg)[2]
params$growth.sigma2=summary(growth.reg)$sigma^2
params$growth.sigma2.exp=as.numeric(growth.reg$modelStruct$varStruct)
#(C) Modify the growth function, g.xy(), to allow the standard deviation
# argument, sd, to be a function of size. This will follow a similar
# pattern to the argument for mean. It’s probably easiest to redefine
# the function g.yx() so it works with code above to create IPMs.
g.yx=function(xp,x,params) {
dnorm(xp,mean=params$growth.int+params$growth.slope*x,
sd=sqrt(params$growth.sigma2*exp(2*params$growth.sigma2.exp*x)))
}
# (D) Rerun the code from sections 1.4.4 and 1.4.5 to obtain λ,
# sensitivities, elasticities, etc.
G=h*outer(y,y,g.yx,params=params) # growth matrix
S=s.x(y,params=params) # survival
P = G # placeholder
for(i in 1:n) P[,i]=G[,i]*S[i] # growth/survival matrix
F=h*outer(y,y,f.yx,params=params) # reproduction matrix
K=P+F # full matrix
(lam=Re(eigen(K)$values[1])) # new population growth rate## [1] 0.980247w.eigen=Re(eigen(K)$vectors[,1])
stable.dist=w.eigen/sum(w.eigen)
v.eigen=Re(eigen(t(K))$vectors[,1])
repro.val=v.eigen/v.eigen[1]
v.dot.w=sum(stable.dist*repro.val)*h
sens=outer(repro.val,stable.dist)/v.dot.w
elas=matrix(as.vector(sens)*as.vector(K)/lam,nrow=n)
par(mfrow=c(2,3),mar=c(4,5,2,2))
image.plot(y,y,t(K), xlab="Size (t)",ylab="Size (t+1)",
col=topo.colors(100), main="IPM matrix")
contour(y,y,t(K), add = TRUE, drawlabels = TRUE)
plot(y,stable.dist,xlab="Size",type="l",main="Stable size distribution")
plot(y,repro.val,xlab="Size",type="l",main="Reproductive values")
image.plot(y,y,t(elas),xlab="Size (t)",ylab="Size (t+1)",main="Elasticity")
image.plot(y,y,t(sens),xlab="Size (t)",ylab="Size (t+1)", main="Sensitivity")