Recently I have been working on massive long-term group-group networks for both the Serengeti lions and Yellowstone wolves. We have tracked territory size, average pack/pride size, the number (and strength) of between pack/pride contacts every year from 1971-until today. Basically a series of time series in which we want to know which one is dependent on which. Not being particularly familiar with time series analysis I didn’t know where to start.

After doing a heap of reading I decided that vector auto-regression was the way to go. Vector autoregression (VAR) are stochastic process models that capture linear dependencies between multivariate time series. Mostly used for economic forecasting, the method seems pretty robust and quite straightforward to implement in the R package ‘vars’. However, finding out all of the steps/assumptions required to run the model was tricky so here is my adapted code to fill the gap:

—————————————————————–

rm(list = ls())

library(“vars”)

#——————————————————————-#

############import data#######################

#——————————————————————-#

data1 <- read.csv(“Data.csv”, head=T)

str(data1)

#——————————————————————-#

############detrend with regression#######################

#——————————————————————-#

m1 <- lm(model~0+Year, data=data1) #lm with no intercept

summary(m1)

m1resid <- residuals(m1)

…..

#make a datframe again

dataResid <- cbind(m1resid)

#——————————————————————-#

############Vector Autoregression#######################

#——————————————————————-#

#make a ts object – Freq here is how many obs per year.

ts.obj <- ts(dataResid,frequency=1, start = 1997, end = 2016); str(ts.obj)

#test for the most appropriate lag for your data (eg., does a 2 year time lag best predict the next years connectivity.

VARselect(ts.obj, lag.max=3, type=”const”)$selection

# ‘p’ below is the the lag factor to test.

varLag1 <- VAR(ts.obj, p=1, type=”const”) #p is is the lag factor

#testing normality (has to be ‘insignificant’ at alpha 0.05 to trust the results)

serial.test(varLag1, lags.pt=10, type=”PT.asymptotic”)

arch.test(varLag1) #test for heteroskedasticity. Error terms are fine if p>0.05

roots(varLag1) #have to be under 1 to trust model results.

#extensive list of summary results.

summary(varLag1)

#links nicely to the forcast package to predict the future

library(“forecast”)

fcstL1 <- forecast(varLag1)

plot(fcstL1, xlab=”Year”)