Multivariate growth model or something else?

[jd_c]

I have some longitudinal data with the following variables, id, time, outcome, stress, and multiple measurements on stress and outcome over time for each id. The question of interest would be, does the change over time (ie growth, trajectory, slope, whatever you want to term it) in stress predict the change over time (growth, trajectory) in the outcome. When I think about this question it almost sounds like, “do the id level slopes in stress predict the id level slopes in outcome?”

I guess I am daft, but I am having a really hard time trying to figure out how to model this question!! (I did try to simulate some data, but I don’t think I ended up simulating quite what I wanted.)

I think what I might want is something called a “multivariate growth model”? But I have had trouble seeing exactly how to formulate it.

Here is what my ideas were:

bf_stress <- bf(stress ~ 0 + (1 + time|p|id))
bf_outcome <- bf(outcome ~ 0 + (1 + time|p|id))
m1 <- brm(bf_outcome + bf_stress + set_rescor(FALSE), data=data)

This would give me the correlation between the varying slopes for stress and the varying slopes for outcome but that isn’t exactly what I was looking for.
Another idea was:

bf_stress <- bf(stress ~ time + (time|p|id))
bf_outcome <- bf(outcome ~ time*stress + (time|p|id))
m2 <- brm(bf_stress + bf_outcome + set_rescor(FALSE), data=data)

This includes stress in the outcome model, but I don’t see that it is really looking at the association of the growth in stress over time with the growth in outcome over time…

I have also thought about running a growth model on stress and one on outcome and then extracting the point estimates and sd’s for the slopes for each and then running a measurement error model using the slopes as data, like this:

m.s <- brm(stress ~ 0 + (1 + time|id), data=data)
m.o <- brm(outcome ~ 0 + (1 + time|id), data=data)
#insert code to extract means and sd's from posterior samples and create new dataframe that is length(id) long
m.sl <- brm(slopes_outcome | mi(sd_outcome) ~ 1 + me(slopes_stress, sd_stress), data=data2)

But that just seems rather hackish and also weird…

Is there a way to answer my question in a single model? I’d be grateful if someone could help! My experience is limited in these sort of longitudinal growth models.

What it seems like I want is for the growth model for stress to be the predictor for the growth model for outcome… but how?

Thanks!

[martinmodrak]

Hi,
I am not sure “growth” models would necessarily be a good match to your problem. Is outcome and stress a continuous variable and only growing? If no, you may want to look at time series literature. Maybe this could be understood as a vector auto-regressive model? (though I don’t think those are supported in brms) Or a dynamic state-space model? (also not supported by brms, but the ctsem package seems to have decent support, although I’ve never used it myself). If outcome is binary, then maybe the joint models supported by stan_jm in rstanarm could also be helpful…

An important consideration is whether time is regularly spaced (which would simplify things a lot). If the time intervals were constant, you could for example use stress in one or two previous measurements as a predictor for the current measurement - which has some limitations, but could be made to run in brms without any problems.

Best of luck with your model!

[jd_c]

@martinmodrak

If the time intervals were constant, you could for example use stress in one or two previous measurements as a predictor for the current measurement - which has some limitations, but could be made to run in brms without any problems.

Thanks for the suggestions! Could you explain a little bit more about the above? So use lagged stress?

I tried to simulate some data, but not sure I simulated quite what I wanted… I would be looking to find beta2.

id <- rep(1:100, each=4)
z <- rnorm(100, 0, 1)
z_s <- rnorm(100, 0, 2)
alpha <- 20
beta <- 5
mu <- (alpha + z[id]) + (beta + z_s[id])*time
sd <- 1
stress <- rnorm(400, mu, sd)

sl_stress <- (beta + z_s[id])
alpha2 <- 1
beta2 <- 3
mu2 <- alpha2 + beta2*sl_stress
sd2 = 1
sl_out <- rnorm(100, mu2, sd2)

alpha3 <- 5
mu3 <- alpha3 + sl_out[id]*time
sd3 <- 1
outcome <- rnorm(400, mu3, sd3)

data <- cbind(id, time, stress, outcome)
data <- data.frame(data)

If you plot it using a 3D plot in plotly, you can see the slope (I don’t have the plotting code or the model code that I tried on this work computer). Usually, if I can simulate some data then I can recover the parameters, because hey, I programmed the generative model…but somehow I never could. Maybe I’m daft.

[simonbrauer]

You might find the (econometric) fixed effects model helpful, as well as its “within-between” growth curve formulation. In the case of a single X predictor, you decompose it into a respondent-level mean (which estimates the between-person variance) and the observation-level deviation from the respondent mean (which estimates the within-person variance). That second term is what indicates the how the change in X is associated with change in Y.

For a multivariate outcome? I haven’t seen this used, but it seems like it could be adapted. This is a shot in the dark, but maybe you can decompose it into three levels (for respondent i and item j at time t)

1. Respondent-level mean of all items across all time points (\bar{X}_i = \frac{\sum_{j=1}^J \sum_{t=1}^T X_{ijt}}{J\times T})
2. Item-level mean deviation from respondent-level mean (\delta_{ij} = \frac{ \sum_{t=1}^T X_{ijt}}{T} - \bar{X}_i)
3. Item x time-level deviation from respondent mean + item-level mean (\delta_{ijt} = X_{ijt} - (\delta_{ij}+- \bar{X}_i))

\bar{X}_i is time-constant and would correspond to how the mean of all X across all time points predicts the mean of Y across all time points. Each respondents’ \delta_{ij} 's are also time-constant and would correspond to how the average difference between item j and the item average is associated with the mean of Y (I think…). And $\delta_{ijt} would correspond to how changes in item j are associated with changes in Y.

[martinmodrak]

@jd_c

Thanks for the suggestions! Could you explain a little bit more about the above? So use lagged stress?

Yes, exactly, I was thinking about using lagged stress as a predictor (potentially you could also use me to handle the measurement error in stress).

I admit that I can’t follow completely what you intended in the simulated code. But I think that unless you go for a state-space model or a vector autoregressive model, you cannot model outcome and stress at the same time (you can put them in the same model, but they would be basically independent components of the model).

So the way to simulate the data for outcome would be just:

id <- rep(1:100, each=4)
time <- rep(1:4, times = 100)
z <- rnorm(100, 0, 1)
z_s <- rnorm(100, 0, 2)

stress_lagged <- rnorm(400, 0, 3) # Could use any distribution or even the actual observed data, because the stress is not modelled directly

alpha <- 20
beta <- 5
beta2 <- 3

mu <- (alpha + z[id]) + (beta + z_s[id])*time + beta2 * stress_lagged
sd <- 1
outcome <- rnorm(400, mu, sd)

You could obviously simulate stress analogously, but as I said, it will actually be an independent submodel.

Would that make sense?

[jd_c]

Thanks for the response. I will have to think about this more, because I don’t readily see how it will get me the effect of the individual level slopes of stress over time on the slopes of outcome over time.

[jd_c]

@martinmodrak

I admit that I can’t follow completely what you intended in the simulated code.

Well, I may have just done a poor job of simulating what I am looking for. Basically, I want to know if the change over time in stress, is predictive of the change in outcome over time. So maybe individuals who’s stress levels increase over time have an outcome that increases over time, and those who’s stress levels decrease have an outcome that decreases over time, irrespective of the actual levels of stress. In other words it is the trend in stress over time that matters for the trend in the outcome over time, not the actual values of stress. Does that make sense? I’m trying to model a scenario like that.
From my understanding, when using the lag, I would just be including predictors based on the previous values in time, but I am not sure that is actually finding the association I am describing in the above paragraph.
I don’t know enough about the dynamic state-space models or the vector auto-regressive models to know if they might do what I describe above.
Thanks for the help!

[martinmodrak]

@jd_c

So maybe individuals who’s stress levels increase over time have an outcome that increases over time, and those who’s stress levels decrease have an outcome that decreases over time, irrespective of the actual levels of stress . In other words it is the trend in stress over time that matters for the trend in the outcome over time, not the actual values of stress. Does that make sense?

Yes that does make sense. I think the simplest way to do this would be to use the change in stress as a predictor either for outcome or for change in outcome. I.e. you compute the diffs from previous time point (provided that the time intervals are the same) and use those in a simple linear regression (where you’ll have one less datapoint per participant as 4 timepoints only give you 3 differences). That’s IMHO the only approach you could do in brms without big hacks.

@jd_c

I don’t know enough about the dynamic state-space models or the vector auto-regressive models to know if they might do what I describe above.

I am also not super knowledgeable about them, but I think they wouldn’t. What you describe sounds like it could be handled by some sort of ODE model, but once again, that would greatly ramp-up the complexity of implementation.

[jd_c]

@martinmodrak

you compute the diffs from previous time point (provided that the time intervals are the same) and use those in a simple linear regression

Ohhh that’s so simple! Why didn’t I think of that?! That does sound like one way to do it. Thanks @martinmodrak ! (also thanks to @simonbrauer for thinking about my problem).

Along those same lines, @martinmodrak what do you think of one of my original ideas where I use this model:
brm(stress ~ 0 + (1 + time|id)
to find slopes and intercepts for each person for stress over time, and then use those slopes and intercepts (and their sd’s from the posterior samples) in a measurement error model to predict the outcome or predict the change in outcome computed either as differences in values from timepoints or as slopes taken from this model:
brm(outcome ~ 0 + (1 + time|id)
??
I think that might have a similar concept as what you suggest but maybe use all of the data in a single predictor for the slope, instead of several predictors for the differences at the different timepoints. Any thoughts on that?

[martinmodrak]

@jd_c

Along those same lines, @martinmodrak what do you think of one of my original ideas where I use this model:
brm(stress ~ 0 + (1 + time|id)
to find slopes and intercepts for each person for stress over time, and then use those slopes and intercepts (and their sd’s from the posterior samples) in a measurement error model to predict the outcome or predict the change in outcome computed either as differences in values from timepoints or as slopes taken from this model:
brm(outcome ~ 0 + (1 + time|id)
??
I think that might have a similar concept as what you suggest but maybe use all of the data in a single predictor for the slope, instead of several predictors for the differences at the different timepoints. Any thoughts on that?

In broad terms, I would advise against this. If I understand you correctly, what you describe is a heuristic approach to implement something like a bigger joint model. I.e.:

1. There is some true, unobserved level of stress (stress_true)
2. The observed stress is a noisy observation of stress_true
3. stress_true is predicted by a linear model
4. stress_true and some form of “trend in stress_true at given time” is then used as a predictor for outcome

If that would describe your assumed model, I would encourage you to use Stan to implement this model (it can’t AFAIK really be implemented in brms without a lot of hacks, which might in the end be more difficult than just using Stan). If I understand you correctly, what you describe is hacky way to implement this in two steps. I would expect it to face similar challenges as e.g. at Using posteriors as new priors The hacky way might give similar results to the full model, but it also might not and there is no easy way to tell - the two step procedure throws a lot of the nice guarantees Stan gives you out of the window.

Additionally, I would note that it is not clear which of the two models (using observed differences as predictors vs. the joint model described above) uses “all of the data” or uses the data “better” in some sense. The joint model IMHO has the potential to be a bit more efficient, but it is also making stronger assumptions about the data by enforcing a linear model for stress_true whereas just using the differences avoids those assumptions (but it makes the additional assumption that stress is interval scaled and thus difference is a meaningful quantity). So if the linear assumption for stress_true is correct, the joint model will “use all the data” but if it is incorrect, the joint model could artificially smooth-out some features of the observed stess levels while using differences as predictors will “use all the data” because it avoids this smoothing.

Finally, if you are going to use differences in stress as predictor, you may also want to use outcome at previous time as predictor (but you usally don’t want to use difference in outcome as response as discussed e.g. Statistical Thinking - Statistical Errors in the Medical Literature)

Does that make sense?

[Charles_Driver]

You probably also want to consider how much you care about the interpretation of the parameters. If it’s not important (ie you just want predictions) then you’re fine, if it is important then a latent variable approach (VAR / state space etc) would be more appropriate.

[jd_c]

@martinmodrak

If I understand you correctly, what you describe is a heuristic approach to implement something like a bigger joint model

That is probably true.

@martinmodrak

There is some true, unobserved level of stress ( stress_true )

Well, I am not sure my idea was really proposing this. The goal was really just to get the linear slopes and use them as predictors in another model. The measurement error part was just to include the uncertainty in the other model. Certainly it would make the assumption that stress changed in a linear fashion over time for each individual.
Everything you say makes sense. I didn’t think the measurement error proposal was sort of a hack for the latent variable model you describe, but maybe I am wrong! It certainly is likely a hack for some sort of joint model. Not having implemented either before, I wouldn’t know.

[jd_c]

Do you know of examples of how these are applied for the problem that I have posed?

[simonbrauer]

@jd_c

The question of interest would be, does the change over time (ie growth, trajectory, slope, whatever you want to term it) in stress predict the change over time (growth, trajectory) in the outcome

@jd_c

individual level slopes of stress over time on the slopes of outcome over time

You are correct that this model wouldn’t use the slope of stress to predict the slope of the outcome. It doesn’t model stress at all. If you are interested in that specifically, then your first model seems appropriate. Instead, it models how the change in stress is associated with the change in the outcome, which sounds like what you are looking for.

The model is mathematically connected to @martinmodrak 's suggestion of using change in stress to predict a change in the outcome. Paul Allison discusses this in his concise and highly recommended book Fixed Effects Regression Models (2009). The limitation to that model is that you can’t include any time-constant predictors of the intercept (such as race and gender). The multilevel model/within-between approach (Allison previously called it the “hybrid” model) allows you to estimate both. I believe that the lagged-variable approach would estimate a slightly different quantity.

Just to illustrate this, I’ve setup some contrived data below for a single predictor (i.e. one measurement of stress). Note the following

1. X and Y are positively correlated
2. The slopes of X and Y are (essentially) uncorrelated
3. The means of X and Y are positively correlated
4. X_dev and Y_dev are negatively correlated (where X_dev = X - respondent mean(X))

So a model predicting Y with X would indicate that they are positively related. A model correlating the slopes of X and Y would indicate that the slopes are unrelated. And a model predicting the change in Y with the change in X would predict that they are negatively associated. To the extent that you can decompose these relationships into different predictors, the better you can understand their unique contributions.

# Setup data
resp <- 100
obsper <- 4
obs <- resp * obsper

alpha <- 2            # Y intercept (fixed effect)
beta_mu <- 2          # Effect of respondent mean of X on Y
beta_dev <- -0.25     # Effect of change in X on Y
beta_y <- 0           # Slope of Y over time (fixed effect)
beta_x <- 0           # Slope of X over time (fixed effect)

x_dev <- c(-1, 1, 1, -1)        # Constrain the X slopes to be 0 across time (on average)

test_df <- data.frame(I = rep(1:resp, each = obsper),                         # Respondent indicator
                      wave = rep(1:obsper, times = resp),                     # Wave/time
                      delta = rep(rnorm(resp), each = obsper),                # Random intercept
                      X_mu = rep(3 + rnorm(resp, sd = 3), each = obsper),     # Respondent mean of X
                      resid = rnorm(obs, sd = 0.25))                          # Residual


test_df$X_dev <- test_df$wave * beta_x + rep(x_dev, times = resp) * rep(rnorm(resp), each = obsper)   # Deviation of X from respondent mean
test_df$X <- test_df$X_mu + test_df$X_dev                                                             # Total value of X
test_df$Y <- alpha +                                                                                  # Total value of Y
             test_df$delta +
             beta_mu * test_df$X_mu +
             beta_dev * test_df$X_dev +
             beta_y * test_df$wave +
             test_df$resid                                                                            
test_df$Y_mu <- rep(aggregate(test_df$Y, by = list(test_df$I), FUN = mean)$x, each = obsper)          # Respondent mean of Y
test_df$Y_dev <- test_df$Y - test_df$Y_mu                                                             # Deviation of Y from respondent mean

# Examine (naive) slope for each individual
y_slopes <- rep(NA, times = resp)
x_slopes <- rep(NA, times = resp)
for(i in 1:resp){
  y_mod <- lm(Y ~wave , data = test_df[test_df$I == i,])
  x_mod <- lm(X ~wave , data = test_df[test_df$I == i,])
  
  y_slopes[i] <- coef(summary(y_mod))[2,1]
  x_slopes[i] <- coef(summary(x_mod))[2,1]
}
test_df$y_slope <- rep(y_slopes, each = obsper)
test_df$x_slope <- rep(x_slopes, each = obsper)

# Plot
plot(test_df$X, test_df$Y)             # X predicting Y, positively associated
plot(test_df$X_mu, test_df$Y_mu)       # Mean of X predicting mean of Y, positively associated
plot(test_df$X_dev, test_df$Y_dev)     # Change in X predicting change in Y, negatively associated
plot(test_df$x_slope, test_df$y_slope) # Slope of X predicting slope of Y, no relationship

# Estimate
wb_results <- lmer(Y ~ X_mu + X_dev + wave + (1 + wave | I), data = test_df)
summary(wb_results

[jd_c]

@simonbrauer nice!! Ok, I see what you are suggesting now, thanks to the great simulation code that made it easy to follow and see what you were describing. It doesn’t specifically include the slope over time for stress, but you’re right, this could be a good way to go about it because it does relate change in stress to the outcome.

If I could mark two solutions, I would mark this as well.

One question about the code - when you write rep(rnorm(resp), each=obsper) in the line that calculates the deviation of X from respondent mean test_df$X_dev, this makes all of the deviations of the same magnitude for each person (which subsequently gives repeat values for the value X for each wave). This wouldn’t be the case in actual data, but I suppose for the model it wouldn’t matter if each X_dev was different for every observation within an individual (that would be a more likely scenario)?

Also, if you didn’t want to assume that the effect of change in X on Y was the same for all individuals, could you add (1 + X_dev + wave | I) in the model?

Thank you for the thorough explanation and code.
Always amazed how helpful this forum is.

[simonbrauer]

Glad to be helpful (as so many here are).

The data are a bit contrived in that the goal was to make the slope of X to be 0 across time. So the code ends up being a bit more confusing than necessary (I had to re-review some of it myself to figure out what I was doing!).

delta corresponds to the random intercept of Y, not the deviation from respondent mean for X. So that is why it is the same for each observation within each respondent. x_dev <- c(-1, 1, 1, -1) is the average change across time for respondents (i.e. the deviations from respondent means). rep(x_dev, times = resp) * rep(rnorm(resp), each = obsper) is meant to produce a random scaling of X_dev while keeping the pattern (and thus the slope) the same.

@jd_c

Also, if you didn’t want to assume that the effect of change in X on Y was the same for all individuals, could you add (1 + X_dev + wave | I) in the model?

That is correct. I believe it should work.