# Disentangling and estimating between-day and between-time variability in a longitudinal study

In a longitudinal study, the hand grip strength of subjects is measured. Each subject is measured at 7 times ($T_1-T_7$):

• $T_1-T_6$ are measurements at different days and different but fixed hours (at 7, 10, 13, 16, 20, 21 hours).
• The starting time was randomized. For example, one subject's $T_1$ could be 13 o'clock whereas another subject would start at 16 o'clock. All participants continued from their starting time according to the sequence. Example: A subject who started at 10 would continue with measurements at 13, 16, 20, 21 and 7 hours.
• $T_7$ is measured at a different day but at the same time as $T_6$. Which time is measured twice is randomized (e.g. for one subject $T_6=T_7=10$ hours, for another $T_6=T_7=21$ hours etc.). This is a direct consequence of randomizing the starting time described above.

At each time, one measurement per subject is taken so that there are 7 data points for each participant.

As I understand it, the measurements at $T_6$ and $T_7$ capture the between-day variability whereas the measurements at $T_1-T_6$ capture the between-day and the between-time variability.

Question: Is it possible to disentangle the between-day and the between-time variability with these data and if so, how?

I thought of using a mixed model or a repeated measures ANOVA but I'm unsure how to set up the model to estimate the two sources of variability separately (if it is possible at all). Thanks for your time.

Edit: Implementing and fitting the model suggested by @MartinModrák using brms:

library(brms)

prior <- c(prior_string("normal(0, 1000)", class = "sd", group = "ID")
, prior_string("normal(0, 100)", class = "sd", group = "Timepoint"))

mod <- brm(legisom ~ (1|Timepoint) + (1|ID)
, data = dat
, prior = prior
, chains = 5
, iter = 50000
, warmup = 10000
, thin = 1
, seed = 142857
, control = list(adapt_delta = 0.99))

summary(mod)


The results are as follows:

Group-Level Effects:
~ID (Number of levels: 18)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)   837.46    154.85   594.11  1196.58      34814 1.00

~Timepoint2 (Number of levels: 6)
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)    55.08     40.91     2.19   152.99     107288 1.00

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept  3971.88    198.37  3580.54  4365.21      25146 1.00

Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma   418.01     29.30   365.28   479.93     136039 1.00


From the output, we can infer the estimates: $\hat{\tau}=837.46, \hat{\kappa}=55.08, \hat{\sigma}=418.01, \hat{\gamma}=3971.88$.

• Are the hours at $T_1 - T_6$ randomized per subject? Or is only $T_6$ randomized? – Martin Modrák Mar 15 '18 at 12:35
• @MartinModrák I have updated the question with the according information. Shortly: The starting time was randomized and as a consequence, the last measurement time which was repeated was also random. – COOLSerdash Mar 15 '18 at 13:45
• One minor point on your model description: do you have a good justification for the priors you set? Note that the estimate for between-person variability (sd(Intercept) of ID) is quite close to the boundary of your prior (N(0,400) means 95% of mass is below 800). This indicates that the estimate might change noticeably if you put a different prior so unless you have good reasons to expect between-person variation below 800, you should try a fit with wider prior to check the robustness of your results. – Martin Modrák Mar 19 '18 at 15:43
• @MartinModrák You're right, chosing a wider prior leads to slightly different results, albeit not too different. No, I don't have a reason for strongly informative priors. – COOLSerdash Mar 19 '18 at 16:49

TLDR: I think you can disentangle in principle, but in practice you would likely need a lot of data to achieve any reasonable precision/reliability. Mixed models are IMHO the way to go.

To me the bigger problem is not disentangling between-day and between-time variability but distinguishing between-day variability and other residual variability not accounted for. If you don't mind mixing the latter two, there is some hope. In the following I will use between-day variability to mean the total variability left after accounting for between-person and between-time variability.

I think mixed models are the way to go. The key is to think generatively - how do you expect your data to have been created? One way to think about the probelm you describe is that there is a per-person true strength $\alpha_i$ which is modified by a Gaussian noise with sd of $\sigma$ independently for each day. Then there is a separate effect for each time point $\beta_t$ which I will assume is the same for all participants (it doesn't have to). $Y_{i,j}$ is the j-th measurement of subject i. Then

$Y_{i,j} \sim N(\alpha_i + \beta_{time(j)} + \gamma, \sigma)$

$\alpha_i \sim N(0,\tau)$

$\beta_t \sim N(0,\kappa)$

Here $\gamma$ is the overall intercept and $\tau$ is the between-person variability, $\kappa$ the between time variability and $\sigma$ combines the measurement error, between-day variability. All of those are treated as model parameters. Since measurement error of your device is likely to be known and quite small and you are actually measuring the variable of interest (not a proxy), you can compute between-day variability simply by subtracting the known measurement error from $\sigma$.

However, your data are only weakly informative about $\kappa$ and $\sigma$ individually (you only have one data point per subject to spot the difference). It might be prudent to use a fully Bayesian approach to quantify this uncertainty, as methods such as lme4 may not be reliable. The above model should be expressible in both INLA and brms which both should handle the uncertainty well (the former is computationally more efficient, while the latter is more flexible and IMHO better documented).

You can also improve your ability to estimate $\beta_t$ by assuming that they are a continuous function of time (e.g. that $\beta_7$ is correlated with $\beta_{10}$ but less bound to $\beta_{21}$). Using splines or Gaussian process or autoregressive models (there is support for some of those in the aforementioned packages). With this assumption, you could probably even work with person-specific between-time variability.

EDIT: It will also help if you have some form of horseshoe prior on the per-time and per-person coefficients e.g.:

$\tau \sim HalfNormal(0,1)$

$\kappa \sim HalfNormal(0,1)$

The priors should obviously be modified to reflect the scale of the values you work with ($HalfNormal(0,1)$ should be OK if your data are somewhere between -5 and +5). Alternatively, you could assume a covariance strucutre on those parameters (once again, this is supported in INLA and brms out of the box, not sure about lme4).

• (+1) Thanks a lot for this comprehensive answer. Am I correct in assuming that the following model corresponds to your model (in lme4 notation): legisom ~ Timepoint + (1|ID)? Where Timepoint is a categorical variable. – COOLSerdash Mar 16 '18 at 11:45
• I have little experience with lme4, but I think the formula would look like legisom ~ Timepoint + (1|ID) + 1  (thats what it would look like in brms and maybe the +1 is implicit) and - very importantly, there is the horseshoe prior (or similar - corresponds to the distribution defined for $\alpha_i$ and $\beta_t$ above) on both the Timepoint and ID intercepts that both ties them together (increasing your power to determine them) and helps you determine the between-person and between-time variability. – Martin Modrák Mar 16 '18 at 12:41
• I also added some notes on priors for the hyperparameters. – Martin Modrák Mar 16 '18 at 13:00
• Thanks that's very helpful information. I managed to fit the model. One question: Because Timepoint is categorical, I think there is a separate prior on each beta. I'm not sure how to extract or see $\tau$ in this. – COOLSerdash Mar 16 '18 at 19:58
• Just to be clear: you are able to get the between-person variability $\tau$, but cannot get the between-time variability $\kappa$? It seems that when working with categorical population-level effects, you cannot easily bind together the individual coefficients. But I think you can have legisom ~ (1|Timepoint) + (1|ID) and then you'll be able to extract both $\kappa$ and $\tau$ (at least this is what is the default in brms, not sure about lme4). – Martin Modrák Mar 19 '18 at 6:14

This could be analyzed via a linear mixed model similar to repeated measures ANOVA. If I understand your design correctly, the layout below shows the order of treatment for days and times for 6 subjects. Below that layout is an example showing data entry for analysis. I admit that I don't understand repeated time from day 6 and 7, so not sure how to address that.

Design layout

Days
1  2  3  4  5  6  7
Subject
1         7 10 13 16 20 21 21
2        10 13 16 20 21  7  7
3        13 16 20 21  7 10 10     (hours within subjects)
4        16 20 21  7 10 13 13
5        20 21  7 10 13 16 16
6        21  7 10 13 16 20 20
etc

Data entry
sub  grip day hour
1    x     1   7
1    x     2  10
1    x     1  13
1    x     4  16
1    x     5  20
1    x     6  21
1    x     7  21
2    x     1  10
2    x     2  13
2    x     1  16
2    x     4  20
2    x     5  21
2    x     6   7
2    x     7   7
etc.


Is it possible to disentangle the between-day and the between-time variability with these data and if so, how?

Yes it is in principle. Mixed models mentioned in the previous answers and comments are based on a principle that the variance observed in a variable affected by compounded variances is the sum of individual sources of variance plus their covariance. Variability is the square root of variance, so total variability is the square root of the sum of the squared variabilities plus covariances.

As I understand it, the measurements at $T_6$ and $T_7$ capture the between-day variability whereas the measurements at $T_1−T_6$ capture the between-day and the between-time variability.

That is a useful first order approximation to start with, but the answer is more subtle than that. It isn’t clear if the rest interval is consistent (i.e. is it every day for a week?). If the interval is one day comparing $T_6$ and $T_7$ will also include the effect of a 24 hour recovery interval (which may or may not be enough depending on your application, the target population and the rigor of the test), while $T_1−T_6$ captures the effect of 27, 28 or 25 hour recovery. Hidden factors could be complicating the picture further, but without explicitly identifying them (daily schedules, weekly schedules, one off activities, medication regimes etc ) you would just lump them all into a heterogenous variability measure (which is after all what propagation of error is all about). This is not necessarily a bad thing, but is important to remember when working out what your can infer from your results.

A further complication is that intuitively one would expect a high degree of covariance between the two variabilities. It would not be surprising if was discovered that subjects with high within-day variability in grip would also exhibit a high between-day grip variability, with tendency the other way for those that exhibit low variability in grip (note I am talking about variability, not average grip characteristics). Your data collection would make it difficult to get a comprehensive estimate of the covariance structure between the sources of variability.

Your experimental design under characterises between-day variability compared to within-day variability, making your estimate of between-day variability and also any covariance more speculative and less well defined. Using this imprecise estimate of between-day variability and covariance to correct you observed variability in order to obtain an estimate of within-day variability will necessarily under-characterise the true within-day variability. You will still get useful estimates that can guide future work, but it is important to understand this principle and not to over infer about either within-day or between-day variabilities based on your results.

If you are doing further work I would recommend a better mix up of the timings so that you get a better balance of repetitions for within and between days. Also, the non randomised nature of the hour would lead to a risk of bias, depending on whether people with the worst grip characteristics have fully recovered from the test previous (not knowing your setup I have no idea of the risk, but randomising would avoid bias). This would allow a more precise estimate of between-day variability and also the covariance and therefore make derivation of within-day variability more accurate.