[R-sig-ME] random slope models with several terms with the same grouping variable

[Hans Ekbrand]

Dear list,

I am trying to model the probability of poverty at the individual
level using predictors at the country level. The dependent variable is
measuring poverty is cdepidxs, and individuals are nested in
household/cluster/country. ("cluster" is a small area consisting of a
few hundred households).

I want to test the effect of three predictors variables at the country
level: loggdp, wbgi_gee, and killed.per.million. In a random intercept
version with these variables as fixed terms they all come out as
significant:

## If you want to load the data set example.data, use 
(load(url("http://dl.dropboxusercontent.com/u/99038959/data/example.data.RData")))

lmer(cdepidxs ~ (1|country) + (1|cluster) + (1|household) + loggdp + wbgi_gee + killed.per.million, data = example.data)

Now, I want to test if the effect of these three predictors vary by a dicotomous grouping variable rural (measuringing whether the cluster is a rural or urban), so I want a random slope for each the conditioned on the variable rural:

lmer(cdepidxs ~ (1|country) + (1|cluster) + (1|household) + (loggdp|rural) + (wbgi_gee|rural) + (killed.per.million|rural), data = example.data)

I ran this model, but I don't understand how to interpret the unique intercept I got for each of them. I thought I might missed to include a term for rural, so I tried a new version with a term for a random interept for rural included on its own:

lmer(cdepidxs ~ (1|country) + (1|cluster) + (1|household) + (loggdp|rural) + (wbgi_gee|rural) + (killed.per.million|rural) + (1|rural), data = example.data)

But the intercept of that term turned out to be almost zero. In fact, it looks like one of the predictors, wgbi_gee "got" all the intercept as the other intercepts are practically zero.

> summary(my.fit)
Linear mixed model fit by REML 
Formula: formula 
   Data: my.df 
     AIC     BIC   logLik deviance REMLdev
 4007625 4007812 -2003797  4007595 4007595
Random effects:
 Groups    Name               Variance   Std.Dev.   Corr  
 household (Intercept)        3.9593e-01 6.2923e-01       
 cluster   (Intercept)        3.3303e-01 5.7709e-01       
 country   (Intercept)        1.7662e-01 4.2026e-01       
 rural     (Intercept)        3.0174e-01 5.4931e-01       
           wbgi_gee           2.0499e-03 4.5276e-02 1.000 
 rural     (Intercept)        4.0165e-25 6.3376e-13       
 rural     (Intercept)        1.7462e-24 1.3214e-12       
           killed.per.million 6.9565e-02 2.6375e-01 0.000 
 rural     (Intercept)        1.4500e-24 1.2042e-12       
           loggdp             2.1524e-01 4.6393e-01 0.000 
 Residual                     2.4374e-01 4.9370e-01       
Number of obs: 1999809, groups: household, 658191; cluster, 39150; country, 69; rural, 2

Fixed effects:
            Estimate Std. Error t value
(Intercept)   4.7075     0.5519   8.529

> ranef(my.fit, whichel = "rural")
$rural
    (Intercept)   wbgi_gee   (Intercept)   (Intercept) killed.per.million   (Intercept)     loggdp
no    -2.011992 -0.1658346 -2.584028e-24 -1.122397e-23        0.001312586 -9.360747e-24 -0.2788045
yes    1.481114  0.1220780  2.583515e-24  1.128187e-23        0.006438184  9.257922e-24 -0.6008917

Should I specifically specify that I don't want a random intercept for each of the variables loggdp, wgbi_gee and killed.per.million? Like this?

lmer(cdepidxs ~ (1|country) + (1|cluster) + (1|household) + (0+loggdp|rural) + (0+wbgi_gee|rural) + (0+killed.per.million|rural) + (1|rural), data = example.data)