3
$\begingroup$

Does it matter that my time series are taken at random points? All the examples I've seen have nice time series, like every day or whatever. I have what is otherwise a pretty simple model: one fixed treatment condition and then a random intercept and slope for each subject, but the observations are essentially random, like s1 has observations at 1,3,7, and 14 days, whereas s2 has them at, say day 5, 6, and 15 (not even the same number of observations!) I was thinking of either binning (early/late), or considering the missing days as missing data. Can someone advise, and/or point me to a relevant example?

P.s., here's the model I'm using:

lmer(KPS ~ Day + Tx + (Day | Tx) + (Day | SID), data = d)
$\endgroup$
  • $\begingroup$ Does this tutorial help? $\endgroup$ – Vincent Guillemot Oct 14 at 21:22
  • $\begingroup$ Is the random slope with respect to days or to some other predictor? $\endgroup$ – EdM Oct 14 at 21:37
  • $\begingroup$ VG: Yes, thanks! (If I want to fill in N/A that looks very useful...although I'm assuming a linear model, so I would have though that this was implied in the technique). EdM: The random slope is wrt subject. (If I understand what you're asking.) I'm going to post dummy data. (Can't post my real data.) $\endgroup$ – jackisquizzical Oct 14 at 23:30
  • $\begingroup$ FYI, It doesn't look like your lmer model includes a random intercept. $\endgroup$ – JTH Oct 15 at 0:27
6
$\begingroup$

Short answer: this should be no problem. I'm pretty sure you'll have somewhat lower power than if you had a perfectly balanced design, but there is no fundamental difficulty.

Here's an example where I subsample the (complete/balanced) sleepstudy data set and show that it works fine (and the results don't change very much).

The only specific issue that I can think of is that if you are fitting a model with autoregressive structure (e.g. in lme with the correlation argument), you'll have to switch from an "AR1" specification to a "continuous AR" or "Ornstein-Uhlenbeck" or "exponential correlation" specification (in lme from corAR1 to corCAR: in glmmTMB from ar() to ou()).

This analysis of carbon exchange in tundra ecosystems shows an example of mixed models applied to an unevenly sampled time series.

subsample data set

library(lme4)
set.seed(101)
table(sleepstudy$Subject)
ss <- do.call("rbind",
           lapply(split(sleepstudy,sleepstudy$Subject),
                     function(x) {
                         x[sort(sample(1:10,
                          size=rbinom(1,size=10,prob=0.7),
                                       replace=FALSE)),]
                     }))
lapply(split(ss,ss$Subject), function(x) x$Days)

fit unbalanced/subsampled data set

lmer(Reaction~Days+(Days|Subject), data=ss)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (Days | Subject)
##    Data: ss
## REML criterion at convergence: 1130.586
## Random effects:
##  Groups   Name        Std.Dev. Corr
##  Subject  (Intercept) 22.634
##           Days         6.303   0.15
##  Residual             21.628
## Number of obs: 119, groups:  Subject, 18
## Fixed Effects:
## (Intercept)         Days
##      250.30        10.44

fit full data set

lmer(Reaction~Days+(Days|Subject), data=sleepstudy)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (Days | Subject)
##    Data: sleepstudy
## REML criterion at convergence: 1743.628
## Random effects:
##  Groups   Name        Std.Dev. Corr
##  Subject  (Intercept) 24.737       
##           Days         5.923   0.07
##  Residual             25.592       
## Number of obs: 180, groups:  Subject, 18
## Fixed Effects:
## (Intercept)         Days  
##      251.41        10.47  
$\endgroup$
  • $\begingroup$ I was intrigued by your comment that the power will be lower for an unbalanced design. I run a short simulation in my answer below that seems to suggest that the power is practically the same. Do you have any further insights on when there will be a difference? $\endgroup$ – Dimitris Rizopoulos Oct 20 at 19:22
  • $\begingroup$ I was guessing (I shouldn't have said "pretty sure"). e.g. "Researchers may also favor equal sample allocation because of a well-known fact that a balanced design produces the maximum statistical power for a given total sample size" (Liu 2003, jstor.org/stable/3701370). A question about balance and power might make a good follow-on CV question ... $\endgroup$ – Ben Bolker Oct 21 at 0:39
4
$\begingroup$

I was intrigued by the comment of @BenBolker that the power will be lower in an unbalanced design compared to a perfectly balanced one. The following simulation study seems to suggest that the power is the same:

simulate_mixed <- function (design = c("balanced", "unbalanced")) {
    design <- match.arg(design)
    n <- 100 # number of subjects
    K <- 8 # number of measurements per subject
    t_max <- 15 # maximum follow-up time

    # we constuct a data frame with the design: 
    DF <- data.frame(id = rep(seq_len(n), each = K),
                     sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))

    DF$time <- if (design == "unbalanced") {
        runif(n * K, 0, t_max)
    } else {
        rep(seq(0, t_max, length.out = K), n)
    }

    X <- model.matrix(~ sex * time, data = DF)
    Z <- model.matrix(~ time, data = DF)
    betas <- c(-2.13, 0.5, 1, -0.5) # fixed effects coefficients
    D11 <- 2 # variance of random intercepts
    D22 <- 1 # variance of random slopes
    D12 <- 0.8 # covariance random intercepts random slopes
    D <- matrix(c(D11, D12, D12, D22), 2, 2)

    # we simulate random effects
    b <- MASS::mvrnorm(n, rep(0, ncol(Z)), D)
    # linear predictor
    eta_y <- drop(X %*% betas + rowSums(Z * b[DF$id, ]))
    # we simulate normal longitudinal data
    DF$y <- rnorm(n * K, mean = eta_y, sd = 1.5)
    DF
}

run_simulation <- function (design, M = 2000L) {
    library("lmerTest")
    opt <- options(warn = (-1))
    on.exit(options(opt))
    p_values <- numeric(M)
    for (i in seq_len(M)) {
        set.seed(i + 2019)
        data_i <- simulate_mixed(design = design)
        fm <- lmer(y ~ sex * time + (time | id), data = data_i)
        p_values[i] <- anova(fm)$`Pr(>F)`[3L]
    }
    mean(p_values < 0.05)
}

#####################################################################
#####################################################################

run_simulation("balanced")
#> [1] 0.691
run_simulation("unbalanced")
#> [1] 0.69
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.