IPM_FRIDAY 15TH_SEPT_2017
OJO PAUL
September 14, 2017
library(grid)
library(Matrix)
library(MASS)
library(nlme)
library(dotCall64)
library(fields)## Loading required package: spam## 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: mapslibrary(IPMpack)
library(spam)
library(maps)
## Load data
getwd()## [1] "C:/Users/Paul/Documents/Rwork"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)
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")
####################################################
####################################################
####################################################
## 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.29676784 -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")
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.