Specify correlation between varying intercepts with same level in different groups

Question

Say I have 2 factor variables foo and bar which both contain the same levels "a", "b", and "c". Is there any way to specify in lme4 (or any other package) a model with random intercepts for foo and bar with correlation between intercepts with the same level? In other words, I think the effect of "a" in foo should be correlated with "a" in bar (similar for "b" and "c"). Formally, this might look something like:

for each level k in ["a", "b", "c"].

Here is some code which estimates sigma^2_foo and 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

but of course misses the correlation term rho. Is it possible to add this correlation structure to this model?

UPDATE

In hope that it is helpful to people who are familiar with Stan, in Stan a basic implementation of this random effects model would look like this:

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

Your model does not have fixed effects and that's why you don't get the correlation matrix. According to your description, you are referring to an interaction between foo and bar at some levels. To add such interaction, you need to add the foo:bar term to the model as a fixed effect, as follows:

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

This will give the following output:

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

Of course, as you might have seen, the interaction here happens between all the levels together (not conditioned as you desire). I hope you can obtain my answer as an initial step to get your desired solution. I'll also try to figure it out for you and will update my answer once I find something useful. My first impression is that you need to modify your dataframe in a way that you can control the interactions between same-level intercepts.

---

[UPDATE]

You can add an interaction variable manually, something like the following:

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

and then you may keep only a with a, b with b, and c with c as follows:

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

You may give it a try and let me know if you need any further help.

Best of luck.