No children

Version 0.4
iugo-cafe home

R script

The zip file for the workshop (upper navbar) has all the R scripts and data for the case studies. This is the file Case_Study_1.r from that zip file.

###################################################
### chunk number 1: RUNFIRST
###################################################
source("KalmanEM.R") # Load the function KalmanEM into work environment by sourcing this script
source("riskfigure.R")
source("TMUfigure.R")

###################################################
### chunk number 6: Cs1_Exercise1
###################################################
par(mfrow=c(3,3))
sim.u = -0.05 # growth rate
sim.Q = 0.01 # process error variance
sim.R = 0.05 # non-process error variance
nYr= 30 # number of years of data to generate
fracmissing = 0.1 # fraction of years that are missing
init = 7 # log of initial pop abundance (~1100 individuals)

years = seq(1:nYr) # col of years
for(i in 1:9){ # for loop over the 9 panels in graph
x = rep(NA,nYr) # creates vector for time series w/o measurement error
y = rep(NA,nYr) # creates vector for the time series w/ measurement error

x=init
for(t in 2:nYr) x[t] = x[t-1]+ sim.u + rnorm(1, mean=0, sd=sqrt(sim.Q))
for(t in 1:nYr) y[t]= x[t] + rnorm(1,mean=0,sd=sqrt(sim.R))
missYears = sample(years[2:(nYr-1)], floor(fracmissing*nYr), replace = FALSE) # 1st and last yrs must be there
y[missYears]=-99
plot(years[y!=-99], y[y!=-99],xlab="",ylab="log abundance",lwd=2,bty="l")
lines(years,x,type="l",lwd=2,lty=2)
title(paste("simulation ",i) )
}
legend("topright", c("Observed","True"),lty = c(-1, 2), pch = c(1, -1))

###################################################
### chunk number 16: Cs1_Exercise2
###################################################
fake.u = -0.05 # growth rate
fake.Q = 0.01 # process error variance
fake.R = 0.05 # non-process error variance
nYr= 30 # number of years of data to generate
fracmissing = 0.1 # fraction of years that are missing
init = 7 # log of initial pop abundance (~1100 individuals)
nsim = 9
years = seq(1:nYr) # col of years
params = matrix(NA, nrow=11, ncol=5, dimnames=list(c(paste("sim",1:9),"mean sim","true"),c("kem.U","den91.U","kem.Q","kem.R", "den91.Q")))
for(i in 1:nsim){
x.ts = matrix(NA,nrow=nsim,ncol=nYr) # creates vector for time series w/o measurement error
y.ts = matrix(NA,nrow=nsim,ncol=nYr) # creates vector for the time series w/ measurement error

x.ts[i,1]=init
for(t in 2:nYr) x.ts[i,t] = x.ts[i,t-1]+ fake.u + rnorm(1, mean=0, sd=sqrt(fake.Q))
for(t in 1:nYr) y.ts[i,t]= x.ts[i,t] + rnorm(1,mean=0,sd=sqrt(fake.R))
missYears = sample(years[2:(nYr-1)], floor(fracmissing*nYr), replace = FALSE) # 1st and last yrs must be there
y.ts[i,missYears]=-99

#Kalman-EM estimates
kem = KalmanEM(y.ts[i,],silent=T)
params[i,c(1,3,4)] = c(kem\$U,kem\$Q,kem\$R)

#Dennis et al 1991 estimates
den.years = years[y.ts[i,]!=-99] # the non missing years
den.yts = y.ts[i,y.ts[i,]!=-99] # the non missing counts
den.n.yts = length(den.years)
delta.pop = rep(NA, den.n.yts-1 ) # create a vector to store transitions
tau = rep(NA, den.n.yts-1 ) # create a vector of time step sizes
for (t in 2:den.n.yts ){
delta.pop[t-1] = den.yts[t] - den.yts[t-1] # store each transition
tau[t-1] = den.years[t]-den.years[t-1] # the transitions
} # end i loop
den91

Created on Oct 21, 2009 at 11:30:13 AM by eli