Estimating correlated parameters with multi-level model

Question

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.

Answer 1

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.

Answer 2

How about just writing out the likelihood function and maximizing?

Answer 3

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.