Multilevel regression, multiple response (multiple outcomes) using lmer function in R?
===

I am trying to run a multilevel regression with multiple responses (multiple outcome measures) using lmer in R. Or at least that is what I think I want to be doing. My key resource is this article by Gelman ([1]), section 5 "Multiple Outcomes and Other Challenges".

I have data on workers from different countries (*Country*) and two different outcome measures: the effort they exerted (*Hours_Work*) and the amount they communicated with others (*Nr_Messages*). Right now, I have two linear models:

fit1 <- lm( Hours_Worked ~ Country -1, m)
fit2 <- lm( Nr_Messages ~ Country -1, m)

From what I read in the literature (and Gelman's work in [1] in particular) I think I am running into some kind of "multiple comparisons" issue, which could/should be eliminated by moving to a multilevel/hierarchical modeling approach. My ultimate goal is to plot some kind of two-dimensional coefficient plot of the different countries using the coefficient for *Hours_Worked* on the x-axis and the coefficient for *Nr_Messages* on the y-axis but somehow take the joint distribution of the two models into account. Does this make sense?

I realize that this help request is rather unspecific but I'm a little bit lost here. I'll quickly respond to any questions asking for clarification.

Chris

  [1]: http://stat.columbia.edu/~gelman/research/published/multiple2f.pdf