重复测量线性回归（R语言）：16名受试者实验数据分析请求

Got it, let's walk through how to run repeated measures linear regression for your dataset in R step by step. Since you have multiple measurements per subject, we need to account for within-subject correlation—standard linear regression won't cut it here because it assumes all observations are independent. Mixed-effects models are the right tool for this job.

First, let's make sure your data is properly structured in long format (which you mentioned it already is—great!). Each row represents one bloodlevelB measurement, with baseline variables (Age, BMI, bloodlevelA) repeated for each subject's rows. Here's how to create a reproducible example matching your dataset structure:

# Load essential packages
library(lme4)       # For mixed-effects models
library(lmerTest)   # To get p-values for fixed effects
library(tidyverse)  # For data manipulation

# Set seed for reproducibility
set.seed(123)

# Create baseline data for 16 subjects
baseline <- tibble(
  subject = paste0("subj_", 1:16),
  Age = sample(40:75, 16, replace = TRUE),
  BMI = sample(22:36, 16, replace = TRUE),
  bloodlevelA = sample(5:13, 16, replace = TRUE)
)

# Generate repeated measures for days 1-5
long_data <- expand.grid(subject = baseline$subject, day = 1:5) %>%
  left_join(baseline, by = "subject") %>%
  # Simulate realistic bloodlevelB values (with subject-specific variation)
  mutate(bloodlevelB = 8 + 0.6*day - 0.12*Age + 0.25*BMI + 0.3*bloodlevelA + rnorm(n(), 0, 1.2))

We'll use a linear mixed-effects model (LMM) to account for both fixed effects (your predictors: day, Age, BMI, bloodlevelA) and random effects (subject-specific variation, since each subject might have a unique baseline bloodlevelB).

Basic Model (Random Intercept Only) #

This model assumes each subject has their own starting bloodlevelB value, but the effect of day is the same across all subjects:

# Fit the model
model_int <- lmer(bloodlevelB ~ day + Age + BMI + bloodlevelA + (1 | subject),
                  data = long_data)

# View model summary
summary(model_int)

Model with Random Slope for Day #

If you suspect the effect of day on bloodlevelB varies by subject, add a random slope for day:

# Fit model with random intercept + random slope for day
model_slope <- lmer(bloodlevelB ~ day + Age + BMI + bloodlevelA + (1 + day | subject),
                    data = long_data)

# Compare the two models to see if the random slope improves fit
anova(model_int, model_slope)

If the p-value from the ANOVA is < 0.05, the random slope model is a better fit.

It's critical to check if your model meets the assumptions of normality and homoscedasticity of residuals:

# Extract residuals and fitted values
resids <- residuals(model_int)
fitted_vals <- fitted(model_int)

# QQ plot to check residual normality
qqnorm(resids)
qqline(resids, col = "red")

# Residual vs fitted values plot to check homoscedasticity
plot(fitted_vals, resids, xlab = "Fitted Values", ylab = "Residuals")
abline(h = 0, lty = 2, col = "red")

If residuals are not normally distributed, consider transforming bloodlevelB (e.g., log transformation) or using a generalized linear mixed model (GLMM) if the data is count/non-normal.

From the summary() output:

• Fixed Effects: Look at the Estimate column to see the direction of each predictor's effect. The Pr(>|t|) column (from lmerTest) tells you if the effect is statistically significant.
• Random Effects: The VarCorr section shows how much variation in bloodlevelB is due to subject differences.
• Model Fit: Use the r.squaredGLMM() function from the MuMIn package to calculate the marginal (fixed effects only) and conditional (fixed + random effects) R-squared:

  library(MuMIn)
  r.squaredGLMM(model_int)
  

If you're focused on testing main effects and interactions (e.g., does day interact with Age?), you can use repeated measures ANOVA with the afex package. Note that ANOVA works best with balanced data (no missing measurements) and categorical predictors, so mixed models are still preferred for your continuous baseline variables:

library(afex)
aov_model <- aov_car(bloodlevelB ~ day * Age * BMI * bloodlevelA + Error(subject/day),
                     data = long_data)
summary(aov_model)

内容的提问来源于stack exchange，提问作者MLyall