# Estimating correlated parameters with multi-level model

I would like to estimate a multi level model in Stata or R (using lmer) where the first level coefficients are the same for all observations, but the coefficients within observation are correlated.

An example would look something like this:

$$Y_i = \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + ... + \varepsilon_{0i}$$ $$\beta_1=\gamma_1 z_1 + \gamma_2 z_2 + \varepsilon_{1}$$ $$\beta_2=\gamma_1 z_1 + \gamma_3 z_3 + \varepsilon_{2}$$ $$\beta_3=\gamma_2 z_2 + \gamma_3 z_3 + \varepsilon_{3}$$ and so on, with equations for each beta.

Clearly, I'd make a distributional assumption for the $\varepsilon$'s... like $\varepsilon \sim N(0,\sigma^2)$

The x variables vary by observation, but the z variables do not vary between observations. Thus, the parameters $\gamma$ and $\beta$ are also the same for all observations.

This differs from most hierarchical models I have seen in that parameters are related within an observation, rather than depending on observation-level characteristics.

As a specific application, consider a model where the dependent variable $Y$ is a student's test scores. The x variables are measures of performance in previous classes, and the $z$ variables are characteristics of those classes. Students have taken the same set of classes, but there are few students in each class, so I'd like to pool estimation of the coefficients $\beta$. Because the classes have similar characteristics, there may be far fewer $\gamma$ parameters than $\beta$ parameters, and pooling estimates to those lower level class characteristics may yield more precise estimates of $\beta$ than estimation without the 2nd level model.

At the same time, I'd like estimates of the $\beta$ parameters, so substituting in and estimating y as a function of $\gamma$ and x only gets me half way there.

What is the best way to estimate this type of model? I typically program in R, Stata and Python.

• why not to start from the hypothesis that second level is described by $\beta_i = c_j + \varepsilon_i$ and to test if all $c_i$ are equal? Is your restrictive assumption supported by data or theory? – Dmitrij Celov Jun 9 '11 at 17:18
• @Dmitrij, I haven't tested it, but in my case I am confident the $\c_i$ are not all equal (or approximately equal). For instance, some classes are very similar to the outcome test, and other classes are not similar. Those that are similar are sure to be better predictors. – DanB Jun 9 '11 at 21:04

Have you tried to use Bugs or Jags, calling one of them from R? The model you seem to be estimating is a simple varying slope model, with predictors at the second level.

Be $i = 1, ...n$ students and $k = 1, ... K$ classes. Assuming your data is in the form student-class (i.e. repeated measures), then your model is:

$y_{i} \sim N(\beta_{[k]}*x_{1,i} + \delta_{[k]}*x_{2,i} +..., \sigma^{2})$

$\beta_{[k]} \sim N(\gamma_{1}*Z_{1,k}, \sigma_{\beta1}^{2})$ ...

This model is quite easy to estimate using Bugs or jags and you can call them with function rjags or bugs. They're in package R2jags and here is a simple example o fitting a multilevel model (with winBugs) on R.

• I was under the impression it would be easier with lmer, but this looks promising and I will try it. Thanks! – DanB Jun 10 '11 at 15:44
• I guess there is a way to do this with lmer, I just don't know. – Manoel Galdino Jun 10 '11 at 18:30

How about just writing out the likelihood function and maximizing?

• I've considered this, and I will do it if there's no better alternative. There are about 60 $\beta$'s and they are formed from about 30 different $\gamma$'s. So it will be long/complicated likelihood function to write out. I'd prefer to take advantage of an existing framework to simplify the code. – DanB Jun 9 '11 at 21:08

How is this advantageous over a normal varying coefficient model such as:

fit<-lmer(score~1+vector of class_attributes+vector of student attributes
+(1+vector of class attributes+vector of student attributes)
+(1+vector of student attributes|class)
+(1+vector of class attributes|student))


?

In this example, there is an overall intercept and attribute effect, but each class has a different coefficient possible which can be viewed by typing ranef(fit)

Section 3.2 of the Bates book on lme4 seems exactly analogous to your situation.

https://r-forge.r-project.org/scm/viewvc.php/*checkout*/www/lMMwR/lrgprt.pdf?revision=656&root=lme4&pathrev=656


Update (I updated the line of code above):

I also ran these lines to try to simulate your situation, but without any student attributes

library(lme4)
n<-100 #class size
pool<-200 #student pool size
class=c(rep(1,n), rep(2,n), rep(3,n))
min_in_class=c(rep(45,n), rep(60,n), rep(90,n))
min_hw=c(rep(90,n), rep(60,n), rep(60,n))
student_id=c(sample(1:pool,n), sample(1:pool,n), sample(1:pool,n))
performance=55+10*class +.1*min_in_class    +.2*min_hw+ -.001*min_in_class*min_hw   +rnorm(3*n, 0,10)
df<-data.frame(class=as.factor(class), min_in_class, min_hw, student_id=as.factor(student_id), performance)
library(reshape2)
melted<-melt(df, id.vars=c('student_id', 'class'))
casted<-dcast(melted, student_id~class+variable)
casted$score<-rowMeans(casted[,c(4,7,10)],na.rm=T)+rnorm(nrow(casted),0,5) df$score<-casted$score[match(df$student_id, casted\$student_id)]


fit<-lmer(score~1+min_in_class+min_hw+(1|class)+(1+min_in_class+min_hw|student_id), data=df)

• You've got unmatched parenthesis in your example fit<- which make it unclear what you have in mind. – Aaron left Stack Overflow Jun 9 '11 at 20:58