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? Commented Jun 9, 2011 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
Commented Jun 9, 2011 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.

I'd rewrite your model as:

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
Commented Jun 10, 2011 at 15:44
• I guess there is a way to do this with lmer, I just don't know. Commented Jun 10, 2011 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
Commented Jun 9, 2011 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)]


I thought what you needed trying to do was this:

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


I ran it with various class sizes and pools and didn't get the results I was expecting; but perhaps with more than a few classes, things will look better.

• You've got unmatched parenthesis in your example fit<- which make it unclear what you have in mind. Commented Jun 9, 2011 at 20:58
• My situation differs from 3.2 in Bates because my first-level characteristics depend on multiple non-nested second-level characteristics. He considers classrooms within schools. I'd like to add characteristics of the class such as subject-matter and how long ago the class occurred. The code you suggest may be promising, but the unmatched parentheses are also making it hard for me to understand.
– DanB
Commented Jun 9, 2011 at 21:30
• I edited the code to add the necessary parentheses at the end. I also added some info after playing with this a bit. Commented Jun 9, 2011 at 22:19
• Do you mean for the last term to be (1+min_in_class+min_hw|student_id) or for it to be (1+min_in_class+min_hw|class)?
– DanB
Commented Jun 10, 2011 at 13:54
• In the little example I ran out, |class would give some trouble with identification (each class has no internal variation on the class level attributes) ; I was having some trouble with unexpected results once I tried the thing on the small data set I simulated, so I am not sure anymore if the formulation with lmer I was trying is helpful or not to you. Commented Jun 13, 2011 at 16:23