# How to measure and visualise relationship between two sets of count data over time?

I have two sets of count data - one for different viral species, and one for different mosquito species. The count data spans 7 weeks, in 3 different locations.

I want to be able to show that changes in the mosquito species counts correspond to changes in the viral species counts. This is how far I've gotten in visualising that relationship:

While this clearly shows what kind of data I have, it is pretty clunky and I haven't actually measured the correlation between the two datasets.

I'm struggling to figure out what kind of statistical test I should use for this. I've looked into a different options - I thought perhaps something like:

cor(as.vector(mosquito_matrix), as.vector(virus_matrix))


But my matrices are of different dimensions and that only gives me one value. I was hoping for something that could show the correlation over time. Is there a statistical test suitable for this?

I am working in R and have been using the Phyloseq package to analyse my data so far.

Modelling the data: You have multivariate count data over time, so you are going to need to look up models that are appropriate for that kind of data. There are a number of different classes of models that have been applied to this kind of data. One commonly used class of models are multivariate extensions to standard count models, such as the negative-binomial GLM (e.g., Davis and Wu 2009; Motta et al 2012; Zhang et al 2017). Another class of models that has recently become popular for this kind of data are copula models (see e.g., Song et al 2008; Nicoloulopolous and Karlis 2009; Madsen and Fang 2010).

I recommend you read some of this literature to get an idea of the available models for this kind of data. So far as I know, there are computational packages in R to implement most of these model forms, but I will leave it to you to check the details. Personally, I have always found that the negative-binomial GLM models most count data well, so I would start by fitting a multivariate negative-binomial GLM. If you are willing to ignore possible auto-correlation over time, one option would be to drop the time variable from the model, and treat the outcomes at different times as exchangeable data points. This should allow you to estimate the statistical association between the counts of mosquito species and virus species.

Visualising the data (and model output): I think you have made quite a good attempt at visualising your data. You are trying to visualise the relationship between three variables at a time (species, time, count), and this can be difficult, so kudos on your attempt. I agree with you that the present plots are quite complicated, and not easy to interpret, so I will offer some suggestions that can simplify them. Probably the most obvious way to simplify this visualisation would be to replace the bar-charts for each time point with a $$5 \times 7$$ heatmap showing the log-counts of each combination of virus and mosquito species. (I am assuming that none of your counts are zero, so a log-scale would not drop any of the values.) This would give you one heatmap showing each time period, and you could then see the change in the pattern in the heatmap over time. Another way to further simplify the visualisation would be to change your time-series plot to an animation. If you made both of these changes then you would end up with a single animation of a $$5 \times 7$$ heatmap that is changing over time. (If you can post your data I might have time to knock up a plot to show you what I have in mind.)

That is how I would suggest you visualise the raw data. However, the relationships between the species counts will probably not be all that clear in the plot of the raw data, and for this latter task you will need to make plots of the outputs of your model. The usual thing to do in these cases is to make barplots of the regression coefficients, with error bars showing confidence intervals. Since the coefficients in the model directly estimate the relationships in the species counts, this will give the viewer/reader a clear visualisation of the estimated relationships in the data.

• Thank you, this is really helpful. For the visualisation, I think the heatmap is a good idea, but I do have some zero counts - what does this mean for the log-scale? Also, does having zero counts effect what models I should use? If it helps, here is my data: Mosquito matrix: figshare.com/articles/mosquito_matrix/10229129 Virus matrix: figshare.com/articles/virus_matrix/10229453 Nov 3, 2019 at 7:33
• Zero counts shouldn't really be a problem. For the visualisation it just means they'll be off the scale and show up as blank (which is fine anyway). Standard models for count data can handle zero counts, so they are certainly not a problem for the model.
– Ben
Nov 3, 2019 at 10:43

I'm sorry that this probably isn't very encouraging, but if you've only got 7 datapoints and you want to perform a statistical test on the relationships between 5 independent variables and 7 dependent variables, your results probably won't be significant. You might have better luck simply checking the correlation between the total mosquito count and the total viral species count.

If all you're interested in is some sensible summary statistics or plots, though, you could try simply evaluating cor(mosquito_matrix[i,], virus_matrix[j,]) for each mosquito species i and virus species j. Again, the correlations will probably be spurious, since you've so few samples and 35 possible combinations i,j.

When you say that you're interested in "correlation over time" I assume you mean checking if, say, an increase in the mosquito population takes a week or two before the viral population increases in response. If so, you'll be interested in a cross correlation plot, which you can get from the ccf function.

I hope this helps!