Linear Mixed Models – Using Group-Level Variables in LMM

lme4-nlmer

I am running an analysis on a national sample of 20,000, representative at the province level (34 provinces)

After checking for linearity and normality of my dependent variable I have run a preliminary OLS in order to see how the covariates perform in explaining the variation of the variable of interest.
I have selected the relevant independent variables following accreditee literature in my field of analysis, explored the covariance matrix in order to avoid problems of multicollinearity etc..

The result from the OLS is good in term of significance level of the coefficients, the sign and the magnitute of the latter fullfil my expectations and match the results find by other analysis.
However, the value of R^2 is quite low: only 0.09 . Thus, knowing that some variation could be explained by the differences between provinces I have first estimate the OLS adding first provincial dummies and after distric dummies (398 districts).

The R^2 improved much, reaching respectively the 36% and 41%.
However, what I would like to see are the underling cause of the regional differences: why do they perform differently?

Among the variables I have some take a unique value according to each observation's province. I cannot use them in the OLS while using the province dummies because there would be perfect collinearity.

In my view using a mixed linear model would help.

I have run a random intercept null.model in which only the dependent variable is regressed against an intercept. For the estimation I have used the command lmer from the {lme4}package in R.

The InterCorrelation Coeffient equals 0.30, suggesting that the 30% of the variation happens between groups, the values of the group-mean reliance are quite high too (not less than 0.9)- I repeat myself: the sample is province representative.

I finally run a set of random mixed intercept model with 2 levels:

where: i indicates the household and j indicates the province.

Y_i = beta0j + beta1 X_i + e_ij
Y_i = beta0j + beta1 X_i + beta2 Z_j + e_ij
Y_i = beta0j + beta1 X_i + beta2 W_j + e_ij

The first question is: am I allowed to include a variable which varies only at the provincial level (as Z and W do) given that their coefficient is not computed by taking in cosideration the groups?
As far as I have understood it would be a mistake to use those variables in the random part of the model in order to get random coefficients.

The second question: given that by running the command anova(), also in lme4, model 1 is statistically different from model 2, and model 2 is statistically not different from model 3, can I say that Z and W have the same power in explaing the variation in Y despite the fact that Z's coefficient is significant and W's is not significant?

Think as if Z and W were proxies for the same dimension. They are in fact statistically and conceptually high correlated.

Sorry but I cannot give more details on the actual problem I am woking on.

Thanks in advance.

Best Answer

am I allowed to include a variable which varies only at the provincial level (as Z and W do) given that their coefficient is not computed by taking in consideration the groups?

Yes. People include group-level covariates (i.e., continuous variables that vary among groups, not within groups) in multilevel/mixed models all the time. The first-order effect, however, is to partition the variance among groups into a component explained by the group-level covariate and a residual component (i.e., the among-group standard deviation in the new model, less than the among-group standard deviation in the original model). It doesn't generally change the predictive accuracy/goodness-of-fit at the individual-observation level very much, although it can improve the reliability of the assumption that the group-level effects are Normally distributed.

given that by running the command anova(), also in lme4, model 1 is statistically different from model 2, and model 2 is statistically not different from model 3

There's a conceptual/theoretical problem here, as models 2 and 3 are not nested; there is no way to reduce model 2 to model 3 but setting one of the parameters to a special value (typically $\{\beta_i\} \to 0$ for some subset of parameters), or vice versa. Thus the likelihood ratio test doesn't apply. There are tests for non-nested models such as Vuong's test, but I personally don't think it's generally a good idea to take null-hypothesis significance tests too seriously when building/selecting among models.
Informally, if model 2 and model 3 have negative log-likelihood values within about 2 units of each other (i.e. less than approximately $\chi^2_1(0.95)/2$), you can say that they're approximately reasonably good fits to the data/Z and W are approximately equal in their predictive capability. Similarly, people often take AIC differences of <2 to mean that models are "approximately equivalent" (in this case, corresponding to negative log-likelihoods within 1 unit).
you're right that the fact that

Z's coefficient is significant and W's is not significant

does not tell you that they have different explanatory power (e.g. see Gelman and Stern, The Difference Between “Significant” and “Not Significant” is not Itself Statistically Significant).

Related Solutions

Solved – Linear mixed effect models with two independent variables

In linear mixed effects mode, definitely you can include more than one independent variable. The general linear mixed model has the form

$y_i=X_i \beta + Z_i bi + \epsilon_i$

where $b_i\sim N(0,D)$, $\epsilon_i\sim N(0,I_{n_i} \sigma^2)$, $X_i$ and $Z_i$ are the known design matrices for fixed($\beta's$) and random effects regression coefficients($bi's$) respectively. $I_{n_i}$ is the $n_i$ dimensional identity matrix and $D$ is the variance co-variance matrix of $b_i$.

So you can include as many independent variable you need in the model and also the random component. But you should be concerned that as you increase the dimension of random effects, the estimation procedure becomes more complex as $b_i$'s are unobserved and so does their variance covariance matrix $D$.

The R codes are as follows

randominterceptfixedslope<-lmer(y ~ x1 + x2 + (1|state),data=data,method="ML") 
randominterceptrandomslope<-lmer(y ~ x1 + x2 + (1+x1+x2|state),data=data,method="ML")

Its not always the same set of covariate contribute both as fixed and random component. That means you can also think of the following model

randominterceptrandomslope<-lmer(y ~ x1 + x2 + (1+x1|state),data=data,method="ML")

randominterceptrandomslope<-lmer(y ~ x1 + x2 + (1+x2|state),data=data,method="ML")

You better use xyplot under lattice to see which covariate can be served in the random component.

Best Answer

Related Solutions

Solved – Linear mixed effect models with two independent variables

Related Question