# How to check for a relationship between water bacterial levels

I have a set of data: http://pastebin.com/KHLKD8XB

which is based on TVC (total viable counts - i.e., numbers of bacteria) from water going into a machine (BM), water taken from the machine (IM) and water taken from a device which is washed in that machine (Scope).

Whilst the counts are performed on a 200ml sample, results are reported as cfu/100ml, hence you appear to find 1/2 a bacteria in some samples.

library(ggplot2)
testscope$date<-as.Date(testscope$date)
ggplot(data=testscope,aes(x=date, y=tvc, colour=class)) +
geom_step()  + facet_grid(class ~ ., scales="free")  +
labs(title="Machine TVC versus Scope TVC") + ylab("cfu/100ml")


Plotted using ggplot2:

My question is: How can I check whether TVC readings from the Scope are related to TVC readings from the machine ?

By eye, I can't discern a pattern between the machine level (IM) and the Scope levels.

Note that although I have the BM data, that is not something we are particularly interested in at present - it is IM versus Scope that we need to understand.

EDIT: Note that the dataset uses natural log, not base 10. My mistake.

• Hi @PJP, welcome to the site! Forgive my ignorance, but if TVC is a count variable, why does it not always come in integers? Count data are usually modelled by a Poisson regression or a Negative binomial regression if your data show overdispersion. – COOLSerdash May 28 '13 at 9:52
• Hi @COOLSerdash. It's because we report in standardised TVC/100ml but we filter larger quantities (I think around 400ml-500ml). So 10 bugs in 400ml =2.5/100ml. – PJP May 28 '13 at 9:55
• That's important information that should be in the question I think. – Glen_b -Reinstate Monica May 28 '13 at 12:09

Granger causality tests if a time series can be better predicted with the additional information from another time series compared to only the time series itself. The problem is that under certain conditions the test does not work, for example if the time series are integrated.

Case 1: Time series are not integrated (easy case) The time series above do not look integrated to me, but be sure to carry out the tests for integration (KPSS, ADF, ...). If they are not integrated, and there is no structural break you can just carry out the standard Granger-causality test. For example you could use the vars package for that. So, set up are VAR model first and then use the causality() function.

Case 2: At least one time series is integrated (more difficult): In this case you could use the Toda-Yamamoto procedure. Another (more demanding) method is the vector error correction model (VECM). HTH