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

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

# Read in the data
# The data are logged already
head(sealData) #look at the first 6 lines
names(sealData) #look at the column names (site names)

#rename columns
years = sealData[,1] #first column is years
sealData=sealData[,2:dim(sealData)[2]] #columns 2:10 are the data
# make plots of the time series (year is in first column)
par(mfrow=c(3,3))
for(i in 1:9) {
plot(years, sealData[,i], xlab="Year", ylab="", main=names(sealData)[i])
}

# The following is an example to show you how to build and run a model
# for one of the hypotheses
# Hypothesis 5 build a simple panmictic model - 1 subpop, 1 R, 1 Q, 1 U
###############################################################
varcov.Q="diagonal"
varcov.R="diagonal"
Q.groups=c(1)
R.groups=rep(1,9) #repeat 1, 9 times
U.groups=c(1)
whichPop=rep(1,9) #all observation time series are measuring the same subpop

kem = KalmanEM(sealData, varcov.Q=varcov.Q, varcov.R=varcov.R, whichPop=whichPop,U.groups=U.groups, Q.groups=Q.groups,R.groups=R.groups)

# the parameter estimates
kem\$U
diag(kem\$Q) #since 0s are on the off-diagonals, we only need to see the diagonal
diag(kem\$R)
#num parameters, loglike, AIC (or use kem\$AICc if you wish)
c(kem\$K, kem\$loglike, kem\$AIC)

# plot the residuals
Xpred = t(kem\$states)
Xobs = sealData
par(mfrow=c(3,3))
for(i in 1:9) {
plot(years, Xpred[,whichPop[i]] - Xobs[,i], ylab="Predicted-Observed Data", main=names(sealData)[i] )
}

Created on Oct 21, 2009 at 11:53:27 AM by eli