Solved – Data visualization of average and standard deviation over a small time series

data visualizationmany-categoriesrtime series

I am trying to find the best way to visualize the following data:

I have values for 3 different times/dates, each time/date has the same 20 species. For each species I have the average height and standard deviation (that was obtained from n observations). I don't have the observations from where the average and standard deviation was obtained.

This is how the data looks like. Species are numerated from 1 to 20.

t1 = cbind(seq(1,20,1),runif(20,0.5,5.0),runif(20,0,1.0))
t2 = cbind(seq(1,20,1),runif(20,0.9,6.0),runif(20,0.2,1.5))
t3 = cbind(seq(1,20,1),runif(20,0.1,8.5),runif(20,0.4,1.0))

Can anyone give any advice on what would be the best plot/graphic to use? I want to show the evolution of the average and standard deviation in time, however I don't want to do a line plot, and for a heat map I have a small dataset.

Best Answer

Why not a line plot? A line plot seems pretty fitting if you'd like to show general trends in mean and SD provided individual specie is not the focal point.

Anyway, here is an alternative, which strictly speaking is still a line plot but time is not one of the axes. It is also good for discussing individual species. x-axis is your variable 1 (I guessed that's mean), y-axis is the variable 2 (SD).

enter image description here

And yes, heat map maybe too much, but I think it's worth a try to show 3D kernel density if you can provide the actual data as well. Codes are taken from this thread.

enter image description here

Another possibility is to make three plots:

t1 SD vs. mean, use black symbol
t1 SD vs. mean, but this time the symbol turned to light grey; overlay t2 SD vs. mean, use black symbol. Connect each species pair with an arrow.
t2 SD vs. mean, but this time the symbol turned to light grey; overlay t3 SD vs. mean, use symbol. Connect each species pair with an arrow.

Related Solutions

Time Series Analysis – Fixed and Random Factor Analysis in R

For mixed-effects modells of time series I usually use the nlme package, because it offers facilities to model auto-correlation structures. If I don't need to consider auto-correlation, I prefer the lme4 package, which offers more flexibility for specifying random effects and is also usually faster.

I recommend reading Zuur et al. 2009 (ISBN 978-0-387-87457-9) as an introduction to mixed effects modelling with R. The book contains a lot of nice and illustrative examples using the nlme package.

From your question and comments the structure of your model should be as follows:

fixed effects: intercept, plant type (make sure that this is a factor variable and not a numeric in R!), time and the interaction between both
random effects: random intercept grouped by plant nested within box, possibly also a random slope vs. time with the same nesting structure
correlation structure: some kind of auto-regressive moving average correlation structure with time as a co-variate and the same grouping structure as the random effect

Thus, a full model could look like this: fit1 <- lme(height ~ type * time, random= ~ 1|box/plant, correlation=corARMA(0.2, form=~time|box/plant, p=1, q=0), data=mydata). You should test by comparing models using the anova function if including a random slope improves the model. You should also test, which auto-correlation structure is most appropriate (although Zuur et al. advise against spending too much effort on finding the optimal auto-correlation structure). Of course, you also need to study residual plots. Possibly you might need to specify a variance structure or transform the dependent if the model suffers from heteroskedasticity. Judging from the plot the relationship between height and time is pretty linear, but you could also try to transform time.

Potential problems: You have an extremely small number of time points, which could result in problems when trying to fit the auto-correlation parameter(s). It could even be impossible to fit them. Also, your number of individuals and boxes is very small.

Solved – Average and standard deviation of timestamps (time wraps around at midnight)

Let's use the simplification you suggest: only use the data from positive readings and disregard the value of the reading, so we are left with a single set of circular data. You can use the circular dispersion, as whuber suggested, possibly multiplied by some constant to determine how much of the data should be seen as an outlier. A good text that is slightly easier to understand than the Wikipedia page would be Statistical Analysis of Circular data, by N.I. Fisher (1995).

I'll give some more straightforward formulas than the Wiki page, and give some sample code.

The dispersion can be calculated as (due to Fisher (p.32-34)):

Denote the data by $\boldsymbol\theta = \{\theta_1, \dots, \theta_n\}.$ An estimate of the mean $\bar\theta$ can be calculated with
- $S=\sum_{i=1}^{n}\sin(\theta_i)$,
- $C=\sum_{i=1}^{n}\cos(\theta_i)$,
- $\hat\mu = \text{atan2}(S, C)$. (See http://en.wikipedia.org/wiki/Atan2)
Calculate $\bar{R} = \frac{\sqrt{S^2 + C^2}}{n}$.
Calculate the dispersion as suggested by whuber. I'm not sure why, but Wikipedia's definition seems to differ slightly from Fisher. I will use Fisher's:
- $\hat\delta = \frac{1 - \left[ (1/n) \sum_{i=1}^{n} \cos 2 (\theta_i - \hat\mu) \right]}{2\bar{R}^2}.$
Then, choose some constant $c$ (1 is probably fine, but you may fine-tune). Then, the interval is given by
- $ \left[\hat\mu - c \hat\delta, \hat\mu + c \hat\delta \right]$.

I know you want to avoid R, but just to show how to calculate this in code, here is some basic R code anyway, which also generates a plot:

n  <- 200
th <- runif(n, 0.5 * pi, 1.5 * pi)

plot(cos(th), sin(th), xlim=c(-1, 1), ylim = c(-1, 1))

S <- sum(sin(th))
C <- sum(cos(th))

mu_hat <- atan2(S, C)

R_bar  <- sqrt(S^2 + C^2) / n

delta_hat <- (1 - sum(cos(2 * (th-mu_hat)))/n) / (2 * R_bar^2)

constant  <- 0.8

CI <- mu_hat + c(-1, 1) * constant * delta_hat

lines(x = c(0, cos(CI[1])), 
      y = c(0, sin(CI[1])), col="green")

lines(x = c(0, cos(CI[2])), 
      y = c(0, sin(CI[2])), col="blue")

Cut-off values for c=.8 with a uniform half-circle from a sample of 200.

In a final note, it may still be better to use the additional information provided by the value of the reading, and not only the sign, because it may provide better estimates. However, the simplification of only using the sign makes the problem much more manageable. If anyone has a good solution that incorporates the readings, I would love to know!

Best Answer

Related Solutions

Time Series Analysis – Fixed and Random Factor Analysis in R

Solved – Average and standard deviation of timestamps (time wraps around at midnight)

Related Question