1
$\begingroup$

I'm working on a longitudinal analysis using the nlme package in R, where I'm modeling the effect of treatment (tx) and menopausal status (menscat) on a biomarker level that tested at two visits: V00 and V30. Also, in my research question, I care about the change of biomarker at two visits, which may be affected by both treatment (tx) and menopausal status (menscat). The biomarker levels are normally distributed.

My model currently includes visit both as a fixed effect and a random effect (random slope), as shown below:

# Simulate the data
set.seed(123)
n <- 100  # Number of subjects
data <- data.frame(
 id = rep(1:n, each = 2),
 visit = factor(rep(c("V00", "V36"), n)),
 tx = factor(rep(c("Placebo", "Active"), each = n)),
 menscat = factor(rep(1:2, each = n)),
 biomarker = rnorm(2 * n, mean = 50, sd = 10)
)

# Model with visit as both fixed and random effects
model <- lme(
 fixed = biomarker ~ tx * visit + menscat * visit,
 random = ~ visit | id,
 data = data,
 method = "REML"
)

summary(model)

My questions are:

  1. Is it statistically appropriate to include visit as both a fixed effect and a random effect?
  2. What are the implications of this model structure? Specifically, how does having visit in both the fixed and random parts of the model affect the interpretation?
  3. Are there any potential issues (e.g., identifiability, multicollinearity) with this approach that I should be aware of?
$\endgroup$
1
  • 2
    $\begingroup$ It would be fine, if you had repeated first visits and repeated second visits for each id. Since you don't have that, you can't include visit in the random effects. The parameters are unidentifiable, which is why lme complains about "Singularity in backsolve at level 0, block 1". $\endgroup$
    – Roland
    Commented Aug 29 at 6:41

1 Answer 1

1
$\begingroup$

One issue in your code is that the way you have tx and menscat variables created, they end up perfectly correlated. I adjusted that bit of code, below, with a few other changes:

n_visits <- 2  
data <- data.frame(
  id = rep(1:n, each = n_visits),
  visit = as.numeric(rep(0:(n_visits-1), length.out = n * n_visits)),
  tx = (rep(0:1, each = n)),
  menscat = rep(rep(0:1, each = n_visits), length.out = n * n_visits),
  biomarker = rnorm(2 * n, mean = 50, sd = 10)
)

The other issue I see with your simulation is that you give no values for fixed effects of (coefficients associated with) treatment, visit, or menscat. Likewise, you do not provide means and standard deviations for the random effects (visit and id). So if the goal is to simulate some data that might look like what you expect/hope to find, you would need to add that information into the simulation.

In terms of your questions, Roland has answered #1 affirmatively. I'll focus on question 2. When you specify a random slope, you are allowing for both an average slope in the population (the fixed slope) and for each individual in your study to have a unique slope that deviates from the population (random slope). You can see how much the model predicts each individual to deviate from both the fixed intercept and slope by using ranef() on your model object in R. For question 3, issues arise when your data does not support a complex random effects specification. You will get warnings and often see random effect estimates that look weird (e.g., near zero variance estimates or a correlation near abs(1)).

$\endgroup$

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.