r - 指定不同组中具有相同水平的不同截距之间的相关性

Question

假设我有 2 个因子变量 foo和 bar两者都包含相同的级别 "a" , "b" , 和 "c" .有什么方法可以在 lme4(或任何其他包)中指定一个具有随机截距的模型 foo和 bar与相同级别的截距之间的相关性？换句话说，我认为"a"的效果在 foo应该与 "a" 相关在 bar (类似于 "b" 和 "c" )。形式上，这可能类似于:



每个级别k在 ["a", "b", "c"] .

这是一些估计 sigma^2_foo 的代码和 sigma^2_bar :

library(lme4)

levs <- c("a", "b", "c")
n <- 1000

df <- data.frame(y = rpois(n, 3.14),
                 foo = sample(levs, n, TRUE),
                 bar = sample(levs, n, TRUE))

mod <- glmer(y ~ (1 | foo) + (1 | bar), df, poisson)

> mod
Formula: y ~ (1 | foo) + (1 | bar)
Random effects:
 Groups Name        Std.Dev.
 foo    (Intercept) 0.009668
 bar    (Intercept) 0.006739

但当然错过了相关项 rho .是否可以将此相关结构添加到此模型中？

更新

希望对熟悉 Stan的人有所帮助，在 Stan 中，此随机效应模型的基本实现如下所示:

data {
    int<lower = 1> num_data;
    int<lower = 1> num_levels;

    int<lower = 0> y[num_data];

    int<lower = 1, upper = num_levels> foo_ix[num_data];
    int<lower = 1, upper = num_levels> bar_ix[num_data];
}

parameters {
    real alpha;

    vector[num_levels] alpha_foo;
    vector[num_levels] alpha_bar;

    real<lower = 0.0> sigma_foo;
    real<lower = 0.0> sigma_bar;

    real<lower = -1.0, upper = 1.0> rho;
}

transformed parameters {
    matrix[2, 2] Sigma;
    Sigma[1, 1] = square(sigma_foo);
    Sigma[2, 1] = rho * sigma_foo * sigma_bar;
    Sigma[1, 2] = rho * sigma_foo * sigma_bar;
    Sigma[2, 2] = square(sigma_bar);
}

model {
    for (i in 1:num_levels) {
        [alpha_foo[i], alpha_bar[i]] ~ multi_normal([0.0, 0.0], Sigma);
    }

    y ~ poisson_log(alpha + alpha_foo[foo_ix] + alpha_bar[bar_ix]);
}

Answer 1

您的模型没有固定效应，这就是您没有得到相关矩阵的原因。根据您的描述，您指的是foo之间的互动和 bar在某些层面上。要添加此类交互，您需要添加 foo:bar项作为固定效应的模型，如下:

mod <- glmer(y ~ (1 | foo) + (1 | bar) + foo:bar, df, poisson)
summary(mod)

这将提供以下输出:

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod]
 Family: poisson  ( log )
Formula: y ~ (1 | foo) + (1 | bar) + foo:bar
   Data: df

     AIC      BIC   logLik deviance df.resid 
  3962.1   4016.1  -1970.1   3940.1      989 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.8572 -0.6665 -0.0947  0.5406  3.8695 

Random effects:
 Groups Name        Variance Std.Dev.
 foo    (Intercept) 0        0       
 bar    (Intercept) 0        0       
Number of obs: 1000, groups:  foo, 3; bar, 3

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.07131    0.05882  18.212   <2e-16 ***
fooa:bara    0.16682    0.07692   2.169   0.0301 *  
foob:bara    0.04549    0.08039   0.566   0.5715    
fooc:bara   -0.08801    0.08464  -1.040   0.2984    
fooa:barb    0.08196    0.08370   0.979   0.3275    
foob:barb    0.05421    0.08006   0.677   0.4983    
fooc:barb    0.08886    0.07712   1.152   0.2492    
fooa:barc   -0.02109    0.07884  -0.268   0.7891    
foob:barc    0.12437    0.07720   1.611   0.1072    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
          (Intr) foo:br fob:br foc:br fo:brb fb:brb fc:brb fo:brc
fooa:bara -0.765                                                 
foob:bara -0.732  0.560                                          
fooc:bara -0.695  0.531  0.509                                   
fooa:barb -0.703  0.537  0.514  0.488                            
foob:barb -0.735  0.562  0.538  0.511  0.516                     
fooc:barb -0.763  0.583  0.558  0.530  0.536  0.560              
fooa:barc -0.746  0.571  0.546  0.519  0.524  0.548  0.569       
foob:barc -0.762  0.583  0.558  0.530  0.535  0.560  0.581  0.569

当然，正如您可能已经看到的，这里的交互发生在所有级别之间(不按您的意愿进行调节)。我希望你能得到我的回答，作为获得你想要的解决方案的第一步。我也会尝试为您解决问题，一旦发现有用的东西，我会更新我的答案。我的第一印象是您需要以一种可以控制同级拦截之间交互的方式修改数据框。

[更新]

您可以手动添加交互变量，如下所示:

df <- transform(df,foo_bar.inter=interaction(foo,bar, sep = ":"))

然后你可以只保留 a与 a , b与 b , 和 c与 c如下:

df$foo_bar.inter[df$foo != df$bar] <- NA

您可以尝试一下，如果您需要任何进一步的帮助，请告诉我。

祝你好运。