10
$\begingroup$

I've been asked to look at a dataset of repeated measures taken during exercise in participants under two conditions control and treatment. Measures were recorded at each minute and at maximum exercise tolerance.

Participants reached different max exercise durations and so have different number of samples. Most range from 3-15 minutes. One outlier went 30+ minutes.

Individual data look like this enter image description here

They are interested in the effect of condition on the outcome variable ic.

What would be an appropriate approach here? Is a linear mixed effects model suited to this? GEE? Something else I know even less about?

An LMEM just feels strange to me. Is it estimating the effects of condition based on assuming that all participants could reach >30 min, and partially pooling an estimate based almost entirely on the outlier id = 16 to fill those "missing" data? enter image description here

model_full <- lmerTest::lmer(ic ~ (I(time) + I(time^2)) * condition + (time | id), data = df)

> summary(model_full)
...
Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 id       (Intercept) 0.34428  0.5868        
          I(time)     0.00163  0.0403   -0.40
 Residual             0.10893  0.3300        
Number of obs: 203, groups:  id, 16

Fixed effects:
                              Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)                   3.65e+00   1.59e-01  1.90e+01   22.94  2.7e-15 ***
I(time)                       4.28e-03   1.67e-02  2.66e+01    0.26   0.8000    
I(time^2)                     3.02e-04   5.81e-04  1.82e+02    0.52   0.6037    
conditiontreatment           -3.89e-01   8.90e-02  1.69e+02   -4.37  2.2e-05 ***
I(time):conditiontreatment    6.80e-02   2.07e-02  1.74e+02    3.28   0.0013 ** 
I(time^2):conditiontreatment -2.34e-03   8.91e-04  1.75e+02   -2.63   0.0093 **

We're not trying to predict individual responses per se, just describe within-subject differences related to condition. Any advice would be appreciated.

Improve this question
$\endgroup$
8
  • 1
    $\begingroup$ Each participant did the exercise just once under the control condition and once under the treatment condition? Have you considered summarizing each trial into a single number of most scientific interest? If you can do that, then you could do a much simple analysis eg. func(ic) ~ condition + (1 | participant). $\endgroup$
    – dipetkov
    Commented Oct 28, 2023 at 23:13
  • $\begingroup$ Thanks. Correct. Yeah, what they did previously was compare RM-ANOVA at complete case time points (only 0-3 min in this case) and paired t-test at peak. But that leaves out a lot of information. So they are interested to describe the trend over time. Would it make any sense to plot individual simple linear models and test slope or quadratic coefficient by condition? $\endgroup$
    – Jem Arnold
    Commented Oct 28, 2023 at 23:32
  • $\begingroup$ In a way the information is limited by the fact that there is one replicate per participant:condition. There is a lot of variability between subjects (so I would ask: is there a trend?) and judging from the plots, the quadratic is not a necessarily a great fit. Note also that you add random intercept and slope per participant but the quadratic effect is fixed. Was there evidence for difference between the conditions "at the peak"? $\endgroup$
    – dipetkov
    Commented Oct 28, 2023 at 23:38
  • $\begingroup$ I was told that this parameter is expected to increase during exercise and fall back toward resting values near maximum. The quadratic function seemed to describe individual data well enough, but ya I'm not committed to that. Comparing linear to quadratic models by AIC & anova() Chisq suggested quadratic was "better". Including I(time^2) in random effects was -0.99 correlated to the linear time term and gave singular fit issues. Paired t-tests were not sig for highest common time (greatest measure shared by both trials). Nor at peak (last measure of each trial). $\endgroup$
    – Jem Arnold
    Commented Oct 29, 2023 at 3:10
  • 1
    $\begingroup$ You can look at the fitted against the measured data points, by participant in time, to get a sense of the fit. The expected behavior doesn't seem to be the observed behavior incl. for the person who exercises the longest. Another point is that "model 2 is a better fit" is not the same as "model 2 is a good fit". $\endgroup$
    – dipetkov
    Commented Oct 29, 2023 at 6:53

3 Answers 3

Reset to default
12
$\begingroup$

A few points:

  • Random effects provide a flexible framework for modeling serial correlations. This is achieved by specifying nonlinear functions of time in the random-effects part of the model. An illustration of this is given in this shiny app (Chapter 3, Section 3.3); also check Section 3.3 in these slides.
  • The shape of the fitted lines you obtained is affected by using polynomials, which do not have a local nature. It would be better to use natural cubic splines. You should also suitably set the boundary knots; for a relevant discussion, check Section 2.4 in my course notes.

EDIT: A comparison between a continuous AR1 structure and random effects

set.seed(1234)
n <- 300 # number of subjects
K <- 6 # number of measurements per subject

# We construct a data frame with the design: 
DF <- data.frame(id = rep(seq_len(n), each = K),
                 time = rep(c(0, 0.5, 1, 1.2, 5, 7.5), n),
                 sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))

# Design matrix for the average longitudinal outcome
X <- model.matrix(~ sex * time, data = DF)

betas <- c(20.1, -0.5, 0.24, -0.05) # fixed effects coefficients
sigma <- 1.2 # error variance
# Continuous AR1 correlation matrix
d <- data.frame(id = rep(1, 6), time = c(0, 0.5, 1, 1.2, 5, 7.5))
cs <- nlme::corCAR1(form = ~ time | id, value = 0.8)
cs <- nlme::Initialize(cs, data = d)
# Covariance matrix 
V <- sigma * nlme::corMatrix(cs)

# linear predictor
eta_y <- matrix(c(X %*% betas), K)
# We simulate normal longitudinal data
DF$y <- c(apply(eta_y, 2, MASS::mvrnorm, n = 1, Sigma = V))


library("nlme")
library("splines")
# ΅We compare the true model with a mixed model that captures the serial 
# correlations using splines
fm <- gls(y ~ sex * time, correlation = corCAR1(form = ~ time | id), data = DF)
gm <- lme(y ~ sex * time, data = DF, random = ~ ns(time, 3) | id,
          control = lmeControl(opt = 'optim'))

intervals(fm, which = "coef")
intervals(gm, which = "fixed")

anova(fm, gm, test = FALSE)
Improve this answer
$\endgroup$
5
  • $\begingroup$ I worry a bit about convergence, Bayesian posterior sampling problems, and effective number of parameters being estimated with that approach. When group level contrasts are the main interest it’s worth less trouble accounting for individual variation. $\endgroup$
    – Frank Harrell
    Commented Oct 30, 2023 at 19:27
  • $\begingroup$ @FrankHarrell if you have a sufficient number of repeated measurements per individual (e.g., 4-5 measurements on average), natural cubic splines with two or three degrees of freedom are pretty stable in my experience. Also, in some datasets, you may select to make the random effects covariance matrix diagonal. This still models the correlations more flexibly than only assuming random intercepts and linear random slopes with an unstructured covariance matrix. E.g., see the example on slides 229-230 in the slides mentioned above. $\endgroup$
    – Dimitris Rizopoulos
    Commented Oct 30, 2023 at 20:26
  • 2
    $\begingroup$ Hi Dimitris - those resources are fantastic by the way. I have a remaining worry. Suppose that the time effects truly linear for every subject. I've seen a paper that derived the covariance matrix induced by a random slopes and intercepts model. This would be applicable to the discussion because the nonlinear spline terms would be negligible. The induced correlation structure is not reasonable. So I still worry that modeling covariance structure with random effects is not the best way. Also there is quite a computational burden and higher effective d.f. with random effects. $\endgroup$
    – Frank Harrell
    Commented Oct 30, 2023 at 21:35
  • 2
    $\begingroup$ Hi @FrankHarrell I've edited my response above to include a simulated example in which longitudinal data are generated by a continuous AR1 structure. I compare the true model with a random effects model using splines. The results for the regression coefficients are very similar, with a log-likelihood of similar magnitude. The fact that we spent more parameters to model the correlations in the random-effects model does not influence the standard errors for the fixed effects. $\endgroup$
    – Dimitris Rizopoulos
    Commented Oct 31, 2023 at 13:07
  • $\begingroup$ Extremely useful. We need to keep monitoring the SE of fixed effects for various dataset when random effects are used. We also need to routinely draw variograms so we can compare correlations from raw data with correlation structure induced by models. $\endgroup$
    – Frank Harrell
    Commented Oct 31, 2023 at 13:42
9
$\begingroup$

I believe that you are fine to use the approach you outlined in your question. There is more between id variance at time==0 (the random intercept) than there is in the rate of change (random slope for time) in the random part of the model. That suggests you might consider eliminating the time random slope, however, I think you are fine to keep it. If you wanted to, you could utilize OLS with indicators for each of the ids in the data. This is sometimes called a fixed effect or no-pooling model. I simulated some data to try to match what I think your data looks like and ran these models. They all converge on the same fixed effect parameter estimates. I'll put my code at the end of the answer, but please note that I did it in Stata as I know how to run simulations better in it than I do in R.

clear *
version 16
set seed 23401
** level 2:
set obs 16
gen id = _n 
matrix C = (1, -.4 \ -.4, 1)
corr2data u_i u_1i, sds(.59 .04) corr(C)
** level 1
expand int(rnormal(6,4))        // observations/id dran from normal distrib
expand 30 if id==16             // for outlier
sort id
bysort id: gen time = _n
expand 2
gen cond = 0
bysort id time: replace cond = 1 if cond[_n-1]==0
replace time = time-1
gen time2 = time*time
*interactions for condition*time and condition*time2
gen condXtime= cond*time
gen condXtime2= cond*time2
gen e_ij = rnormal(0,.33)

* outcome
gen ic = 3.65 + .00428*time + .0003*time2 + (-.389)*cond +  ///
    .068*condXtime + (-.00234)*condXtime2 + u_i + time*u_1i + e_ij

* data structure for a single individual
list id time cond ic if id==11, noobs sep(6)
/*
+-----------------------------+
| id   time   cond         ic |
|-----------------------------|
| 11      0      0   3.512198 |
| 11      0      1   4.249703 |
| 11      1      0   4.185723 |
| 11      1      1   4.209756 |
| 11      2      0   4.182796 |
| 11      2      1   4.213877 |
+-----------------------------+
*/

* Models
mixed ic c.time##c.time##i.cond || id: time, cov(un) reml stddev
est sto full_mixed 

mixed ic c.time##c.time##i.cond || id: , reml stddev
est sto reduced_mixed 

regress ic c.time##c.time##i.cond i.id
est sto reg 

est table full_mixed reduced_mixed reg,  b(%7.3f) se(%7.3f)

/*
--------------------------------------------
    Variable | full_~d   reduc~d     reg    
-------------+------------------------------
ic           |
        time |  -0.004    -0.008     -0.007     // coeff
             |   0.016     0.013      0.013     // SE
             |
      c.time#|
      c.time |   0.001     0.003      0.003     // coeff
             |   0.001     0.001      0.001     // SE
             |
        cond |
          1  |  -0.411    -0.411      -0.411    // coeff  
             |   0.067     0.070       0.070    // SE 
             |
 cond#c.time |
          1  |   0.081     0.081       0.081   // coeff
             |   0.017     0.018       0.018   // SE
             |
 cond#c.time#|
      c.time |
          1  |  -0.003    -0.003      -0.03    // coeff
             |   0.001     0.001       0.01    // SE
             |
       _cons |   3.648     3.643       2.810   // coeff   
             |   0.143     0.143       0.093   // SE
             |
  id(dummies)|   --        --          included
-----------------------------------------------
*/
Improve this answer
$\endgroup$
9
$\begingroup$

I have found that for data similar to yours the fit of the model to the actual correlation patterns is crucial. For example if you use random effects and the induced correlation pattern is not what the random effects model assumes, the random effects can go haywire and inference for the other parameters (especially standard errors) can be distorted. Random effects models also have no way to taking absorbing states into account. This would be relevant to you if exercise had to be terminated for some bad reason such as pain, but may also be relevant with regard to exhaustion.

Generalized least squares and Markov models are two ways to model serial correlation flexibly. If you find that results are not interpretable because of censoring (exhaustion/reached exercise tolerance) then a Markov multi-state model with absorbing states may be applicable. The previous link shows how to use Markov ordinal longitudinal modeling for an ordinal response variable. This ordinal response can include hundreds of levels of a continuous response, with additional discrete levels at the high end to capture events that may or may not be terminating/absorbing.

See this for a related discussion. There is no reason to automatically use random effects for longitudinal data.

Improve this answer
$\endgroup$
12
  • $\begingroup$ In this dataset there are n = 16 participants. How likely is it that the Markov model can be fitted in the first place? The examples in Ch 22 have 100x patients. $\endgroup$
    – dipetkov
    Commented Oct 29, 2023 at 13:45
  • $\begingroup$ The effective sample size for the Markov model is greater than the effective sample size of the methods you are entertaining, because it’s more likely to fit the correlation structure. The only place the Markov model will have an inadequate effective sample size would be if you modeled exhaustion as the most extreme Y and you wanted to make a separate assessment concerning exhaustion without borrowing of any information about pre-exhaustion. $\endgroup$
    – Frank Harrell
    Commented Oct 29, 2023 at 14:02
  • $\begingroup$ The method I'm entertaining (well, proposing to the OP to consider) is to devise summaries that describe ic from each exercise run (so one number per run). I asked the question about $n$ precisely because estimating correlations (or even more challenging, the exercise curves) seems a bit of a non-starter to me. $\endgroup$
    – dipetkov
    Commented Oct 29, 2023 at 14:05
  • $\begingroup$ This proposal is inspired by how bioequivalence trials are analyzed. $\endgroup$
    – dipetkov
    Commented Oct 29, 2023 at 14:14
  • 1
    $\begingroup$ Inclusion of corAR1 will likely make the model fit a lot better. The question then is whether it’s worth all the trouble and approximations to add random effects to that. $\endgroup$
    – Frank Harrell
    Commented Oct 29, 2023 at 18:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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