Time-series analysis of population monitoring data

Navigation

No children

Time-series analysis of population monitoring data

Login

Search


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_2.r from that zip file.

###################################################
### chunk number 1: RUNFIRST
###################################################
source("KalmanEM.R") # Load the function KalmanEM into work environment

###################################################
### chunk number 3: Cs2.readindata
###################################################
d = read.csv("Case_Study_2_data.csv",header=TRUE)
years = d[,1] #[,1] means all rows, column 1
dat = d[,2:ncol(d)]
n = ncol(dat)
#set up legend names
legendnames = (unlist(dimnames(d)[2]))[2:ncol(d)]
for(i in 1:length(legendnames)) cat(paste(i,legendnames[i],"\n",sep=" "))

###################################################
### chunk number 5: Cs2-analysis 1
###################################################
whichPop = c(1,1,1,1,1)
varcov.R = "diagonal"
kem1 = KalmanEM(dat, whichPop=c(1,1,1,1,1), varcov.R="diagonal")

###################################################
### chunk number 5: Cs2-Output parameters
###################################################
# to use a different analysis, replace kem1 with kem2 or kem3 ..
kem=kem1
kem$U
kem$Q
kem$R
kem$loglike
kem$K
kem$AIC

###################################################
### chunk number 10: Cs2-plot data and population estimates
###################################################
kem=kem1
matplot(years, dat,xlab="",ylab="index of log abundance",pch=c("1","2","3","4","5"),ylim=c(5,9),bty="n")
lines(years,kem$states-1.96*kem1$states.se,type="l",lwd=1,lty=2,col="red")
lines(years,kem$states+1.96*kem1$states.se,type="l",lwd=1,lty=2,col="red")
lines(years,kem$states,type="l",lwd=2)
title("Observations and total population estimate",cex.main=.9)

###################################################
### chunk number 20: plot residuals
###################################################
kem=kem1
plotdat = dat; plotdat[plotdat == -99] = NA;
matrix.of.biases = matrix(kem$A,nrow=nrow(dat),ncol=ncol(dat),byrow=T)
xs = matrix(kem$states,nrow=dim(plotdat)[1],ncol=dim(plotdat)[2],byrow=F)
resids = plotdat-matrix.of.biases-xs
par(mfrow=c(2,3))
for(i in 1:n){
plot(resids[!is.na(resids[,i]),i],ylab="residuals")
title(legendnames[i])
}
par(mfrow=c(1,1))

###################################################
### chunk number 17: Cs2-analysis 2
###################################################
kem2 = KalmanEM(dat, whichPop=c(1,1,1,1,1), varcov.R="diagonal", R.groups=c(1,1,1,1,1))

###################################################
### chunk number 19: Cs2-compare AICs
###################################################
c(kem1$AIC,kem2$AIC)


###################################################
### chunk number 21: Cs2-analysis 3
###################################################
Z = c(1,1,2,2,2)
U.groups = c(1,1)
Q.groups = c(1,1)
R.groups = c(1,1,1,1,1)
varcov.Q = "diagonal"
kem3 = KalmanEM(dat, whichPop=Z, varcov.Q="diagonal", varcov.R="diagonal", U.groups=U.groups, Q.groups=Q.groups, R.groups=R.groups)



###################################################
### For analyses 4 on, here is the form of the KalmanEM call
###################################################
## kem = KalmanEM(dat, whichPop=Z, varcov.Q="diagonal", varcov.R="diagonal", U.groups=U.groups, Q.groups=Q.groups)

## The Z, U.groups, and Q.groups for analyses 4-7 are below if you get stuck



###################################################
### chunk number 27: Cs2-case4
###################################################
Z=c(1,2,3,4,5)
U.groups=c(1,1,1,1,1)
Q.groups=c(1,1,1,1,1)


###################################################
### chunk number 28: Cs2-case5 eval=FALSE
###################################################
Z=c(1,1,2,2,2)
U.groups=c(1,2)
Q.groups=c(1,2)


###################################################
### chunk number 29: Cs2-case6 eval=FALSE
###################################################
Z=c(1,1,1,1,2)
U.groups=c(1,2)
Q.groups=c(1,2)


###################################################
### chunk number 30: Cs2-case7 eval=FALSE
###################################################
Z=c(1,1,2,2,3)
U.groups=c(1,2,3)
Q.groups=c(1,2,3)



Created on Oct 21, 2009 at 11:49:23 AM by eli

No comments