0
$\begingroup$

How can a linear mixed-effects modeling approach using the lmer() function be applied to investigate the relationship between a single-measured independent variable (covariate) and a dependent variable with multiple measurements. What considerations should be taken into account when interpreting the results?

$\endgroup$

2 Answers 2

2
$\begingroup$

Yes, you can use a multilevel model in a situation like this. However, you can only include a random intercept, not a random slope, because your predictor would only have one value per cluster. In lmer

model<-lmer(y ~ (1|cluster) + covariate, data = data)

This model will give you the estimate of (between-person differences in) covariate predicting y while controlling for the non-independence of your y observations.

An example would be a study where participants' mood is measured, say, 5 times a day for 10 days (i.e. 50 times). We also measure their trait Anxiety (only once per participant, because we consider it a trait). We then use a multilevel model with a random intercept of participant for predicting mood from trait anxiety because the mood observations are non-independent. The estimate for trait Anxiety would tell you how strongly between-person differences in trait Anxiety are related to momentary mood scores. Say it's negative - we'd then conclude that on average, people with higher trait Anxiety experienced more negative mood during the follow-up 10 days than people with low trait Anxiety.

Edited to add: I realized I don't actually know if the results of the above would be any different from a single-level model with (e.g.) trait Anxiety predicting average mood (averaged over the 50 occasions). I ran some tests with actual data and the results were almost identical, but I'm not sure about the mathematical differences.

$\endgroup$
5
  • $\begingroup$ Thank you for the answer. I see that you have not included any fixed effects Is there any specific reason to do so? $\endgroup$ Aug 30 at 12:48
  • 1
    $\begingroup$ In the model code above, there is a fixed effect of "covariate" included. This is how fixed effects are specified in lmer code (in the multilevel model framework, "fixed effect" refers to a predictor). $\endgroup$
    – Sointu
    Aug 30 at 13:28
  • 1
    $\begingroup$ ...so the "covariate" in the code is just meant to refer to the name of your independent variable. $\endgroup$
    – Sointu
    Aug 30 at 13:36
  • $\begingroup$ I see it now. Missed it earlier. Thank you. $\endgroup$ Aug 30 at 14:17
  • $\begingroup$ @OmarRafique if you felt that Sointu addressed your questions, you should hit the checkmark next to the solution to accept the answer. $\endgroup$ Sep 4 at 14:02
1
$\begingroup$

Sointu's answer is already correct, but I wanted to contextualize his answer. Below is some code for a simulated model which shows a bivariate relationship with a clustered mean by year.

#### Load Libraries #### ####
library(tidyverse)
library(lmerTest)
library(lattice)

#### Simulate Data ####
set.seed(123)
x <- rnorm(n=10000)
year <- rep(1:10,1000)
y <- x + year + rnorm(n=10000)
df <- tibble(x,y,year)

#### Plot By Year Averages of Y ####
df %>% 
  ggplot(aes(x=factor(year),y=y))+
  geom_boxplot()+
  labs(x="Year",
       y="Y",
       title = "Average Value of Y by Year")+
  theme_classic()

#### Plot Bivariate Relationship by Year ####
df %>% 
  ggplot(aes(x,y))+
  geom_point(color="gray",
             alpha = .2)+
  geom_smooth(aes(x,y,group=year),
              se=F,
              color="steelblue",
              method = "lm")+
  theme_classic()+
  labs(x="X",
       y="Y",
       title = "Simulated By-Year Regressions")

#### Fit Model ####
fit <- lmer(y ~ x + (1|year))

#### Get Random Effects and Plot ####
ranef(fit)
ranef(fit)
dotplot(ranef(fit))

If we assume your time points are years for this simulated case, the reason why it works is because your $y$ value here exists across every time point, but the time points are technically clustered by year into these bins, where their means and variance are freely observed:

enter image description here

As such, you can simply enter this by-year factor variable in as a random intercept in a Gaussian mixed model, wherein:

$$ \begin{aligned} \operatorname{y}_{i} &\sim N \left(\alpha_{ij} + \beta_{1}(\operatorname{x}), \sigma^2 \right) \\ \alpha_{j} &\sim N \left(\mu_{\alpha_{j}}, \sigma^2_{\alpha_{j}} \right) \text{, for year j = 1,} \dots \text{,J} \end{aligned} $$

so that your dependent variable $y$ is estimated by a conditional mean $a$ which is predicted by $i$ individuals and $j$ years and has a fixed variance $\sigma^2$. What this means is that each year technically has its own conditional mean value of $y$ called $a_j$, or the random intercept.

An example is below, where the simulated data encompasses a single bivariate relationship that varies by year. You can clearly see that in general the trend is positive, but the baseline conditional mean for each year varies a lot (though the random slopes do not vary at all). Now visualizing our single vector of $x$ values and our single vector of $y$ values regressed on by each year, we get a relationship like this (in the example I have simulated here):

enter image description here

Here you can see in general the relationship between our two variables is universally positive, but the conditional mean/intercept varies a lot. Year 10, for example, has much higher values of y on average and thus the regression line is at the top of our scatterplot.

Mixed models allow you to directly observe these clusters either through the ranef command or through the use of caterpillar plots. As an example with this data:

> ranef(fit)
$year
   (Intercept)
1   -4.4466178
2   -3.4817729
3   -2.5182757
4   -1.4839000
5   -0.5300438
6    0.4787734
7    1.5254932
8    2.5065977
9    3.4858812
10   4.4638646 

We can see that the conditional mean of Year 10 is indeed the highest, where the opposite is true for Year 1, as

$$ a_1 = -4.4466178 \\ a_{10} = 4.4638646 $$

The positive values represent random intercepts which are greater than the fixed intercept, and the negative values are lower than the fixed intercept or the grand mean. A caterpillar plot will showcase this visually.

enter image description here

Now of course this is a highly unrealistic example, as its probably not likely that each year will have almost exactly the same slopes and variance around the mean, but this hopefully shows what the other answerer was explaining.

$\endgroup$

Your Answer

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

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