R：对LME4或BRMS对象的随机截距进行建模

Question

是否有某种方法可以直接（共同）建模使用lme4's lmer（）或brms估计的随机截距？例如，在以下代码中，我符合分层模型，提取随机截距，然后对它们进行建模。

这种两步方法的一个缺点是，我忽略了拦截中的错误。使用强大的协方差矩阵，加权最小二乘等方面很容易固定。

对于上下文：我对此很感兴趣，因为我正在估计一个项目响应模型，其中每个随机拦截是每个时间点的能力，我想预测这些能力。我将在更复杂的贝叶斯模型中做所有这一切。

library(lme4)
library(tibble)
set.seed(123)

# Simulate longitudinal data
N <- 100
time <- 2

# Time-varying data
df <- tibble(person = rep(1:N, time),
           x = rnorm(N*time),       
           y = 2 + x*runif(N*time)) 

# Fit hierarchical model
mod <- lmer(y ~ -1 + (1 | person), df)

# Time-invariant data (constant within person)
df_person <- data.frame(ints = data.frame(ranef(mod))$condval,
                        sex = rbinom(N, size = 1, prob = 0.5))

# Model intercepts as function of time-invariant feature
summary(lm(ints ~ 1 + sex, df_person))

Answer 1

我不知道lme4或brms，但可以直接在Stan中完成。这是复制模型的一种方法，但所有内容都共同估计；参见这个示例（（（在python中，不是r）和 。

模型

我们用副人的随机截距对观察到的值建模（一个对于每个人 j ）和一个比例参数。

我们对副人的随机截距对人的性别的函数进行建模，并且系数，加上第二个刻度参数。

stan代码

这是上述模型的stan代码。它使用非中心的参数化： 具有标准的正常先验，我们转换为<<< img src =“ https://latex.codecogs.com/svg.image?%5calpha_j” alt =“ alpha”>使用和。

data {
  int<lower=0> N_obs; // number of observations
  int<lower=2> N_person; // number of persons
  int<lower=1,upper=N_person> person[N_obs]; // person associated with each observation
  vector<lower=0,upper=1>[N_person] sex; // sex of each person
  vector[N_obs] y; // observed outcomes
}

parameters {
  real<lower=0> sigma; // observation-level variation
  vector[N_person] alpha_person_raw; // random intercepts for persons
  real mu_alpha; // mean of random intercepts for persons
  real<lower=0> sigma_alpha; // scale of random intercepts for persons
  real beta; // coefficients for sex
}

transformed parameters {
  vector[N_person] alpha_person = (alpha_person_raw * sigma_alpha) +
                                  mu_alpha +
                                  (sex * beta);
}

model {
  sigma ~ exponential(1);
  alpha_person_raw ~ std_normal();
  mu_alpha ~ std_normal();
  sigma_alpha ~ exponential(1);
  beta ~ std_normal();
  y ~ normal(alpha_person[person], sigma);
}

我使用与原始示例相同的规则重新创建了数据集

，但是以稍微不同的格式，Stan可以使用略有不同的格式。

library(tidyverse)
person.df = data.frame(person = 1:N, sex = rbinom(N, size = 1, prob = 0.5))
obs.df = data.frame(person = rep(1:N, time),
                    x = rnorm(N * time)) %>%
  mutate(y = 2 + (x * runif(N * time)))

library(rstan)
stan.data = list(
  N_obs = nrow(obs.df),
  N_person = nrow(person.df),
  person = obs.df$person,
  sex = person.df$sex,
  y = obs.df$y
)
stan.model = stan("two_level_model.stan", data = stan.data, chains = 3)

首先比较两种方法

，我使用新数据集重新构建了两阶段模型。

first.level.m = lmer(y ~ -1 + (1 | person), obs.df)
second.level.m = lm(intercept ~ 1 + sex,
                    person.df %>%
                      mutate(intercept = ranef(first.level.m)$person[["(Intercept)"]]))

现在，让我们比较估计的参数值。两种方法都同意，性的效果包括0，随机拦截的平均值约为2（自然，这些估计值受x尚未作为预测器的事实的影响。在模型中。）它们在观察级别估计了大约相同的误差项，但是联合模型估计随机截距模型的误差项大大较小。我不确定lme4的模型拟合过程如何与指数级相关，请先发。

library(tidybayes)
bind_rows(
  spread_draws(stan.model, mu_alpha, beta, sigma_alpha, sigma) %>%
    ungroup() %>%
    dplyr::select(.draw, mu_alpha, beta, sigma_alpha, sigma) %>%
    pivot_longer(cols = -.draw, names_to = "parameter") %>%
    group_by(parameter) %>%
    summarise(lower.95 = quantile(value, 0.025),
              lower.50 = quantile(value, 0.25),
              est = median(value),
              upper.50 = quantile(value, 0.75),
              upper.95 = quantile(value, 0.975)) %>%
    ungroup() %>%
    mutate(parameter = case_when(parameter == "mu_alpha" ~ "mu[alpha]",
                                 parameter == "beta" ~ "beta",
                                 parameter == "sigma_alpha" ~ "sigma[alpha]",
                                 parameter == "sigma" ~ "sigma"),
           model = "joint"),
  summary(second.level.m)$coefficients %>%
    data.frame() %>%
    rownames_to_column("parameter") %>%
    mutate(lower.95 = Estimate + (Std..Error * qt(0.025, second.level.m$df.residual)),
           lower.50 = Estimate + (Std..Error * qt(0.25, second.level.m$df.residual)),
           est = Estimate,
           upper.50 = Estimate + (Std..Error * qt(0.75, second.level.m$df.residual)),
           upper.95 = Estimate + (Std..Error * qt(0.975, second.level.m$df.residual)),
           parameter = case_when(parameter == "(Intercept)" ~ "mu[alpha]",
                                 parameter == "sex" ~ "beta"),
           model = "two-stage") %>%
    dplyr::select(model, parameter, est, matches("lower|upper")),
  data.frame(est = summary(second.level.m)$sigma,
             parameter = "sigma[alpha]",
             model = "two-stage"),
  data.frame(est = summary(first.level.m)$sigma,
             parameter = "sigma",
             model = "two-stage")
) %>%
  mutate(parameter = fct_relevel(parameter, "beta", "mu[alpha]", "sigma[alpha]",
                                 "sigma")) %>%
  ggplot(aes(x = parameter, color = model)) +
  geom_linerange(aes(ymin = lower.95, ymax = upper.95), size = 1,
                 position = position_dodge(width = 0.3)) +
  geom_linerange(aes(ymin = lower.50, ymax = upper.50), size = 2,
                 position = position_dodge(width = 0.3)) +
  geom_point(aes(y = est), size = 3, position = position_dodge(width = 0.3)) +
  scale_x_discrete(labels = scales::parse_format()) +
  labs(x = "Parameter", y = "Estimated parameter value", color = "Model") +
  coord_flip() +
  theme_bw()