Generally it is best to send questions like this to the
R-SIG-Mixed-Models@R-project.org mailing list (see instructions at
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models).  Many of
those who read that list can reply and often much faster than I am
able to.

I have taken the liberty of cc:ing the list on this reply.



On Mon, Aug 26, 2013 at 4:57 PM, Lijun Chen <lchen@chapinhall.org> wrote:

>  Dear Professor Bates,
>
> My name is Lijun Chen, I am researcher with Chapin Hall at the University
> of Chicago.
>
> Recently I am learning to use R to analyze some administrative data on
> foster children in IL. I came across some documentation you wrote about
> using R for multilevel data analysis. That is why I am writing to seek some
> advice on using R for poisson models.
>
> One of the tasks in my analysis is to learn about how each county in IL
> stand in the number of children placed in public foster care in a year. For
> each county, we have the children placed in foster care, its total child
> population, and a bunch of county-level attributes from the census data.
>
> I decided to adopt the poisson model under the lme4 package in R to run a
> mixed effect poisson model. Here is the syntax I use:
> >
> fit2=lmer(fp~1+offset(log(childn))+(1|county),family=poisson,data=ildata,reml=TRUE)
> Here are the results:
> > fit2
> Generalized linear mixed model fit by the Laplace approximation
> Formula: fp ~ offset(log(childn)) + (1 | county)
>    Data: ildata
>    AIC   BIC logLik deviance
>  334.1 339.3   -165    330.1
> Random effects:
>  Groups Name        Variance Std.Dev.
>  county (Intercept) 0.53857  0.73388
> Number of obs: 101, groups: county, 101
>
> Fixed effects:
>             Estimate Std. Error z value Pr(>|z|)
> (Intercept) -6.25386    0.07988  -78.29   <2e-16 ***
> --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> However, I found that the R results (coefficients and variance) are
> different from the coefficients obtained with the HLM package (-6.2152
> unit-specific, 0.5324 variance). I wonder whether this is because my
> specification of the model in R is incorrect, or it is because estimation
> methods used are different in R and HLM. I tried out both Full Maximum
> Likelihood and REML in HLM, and both sets of results are different from the
> R results.
>
> To get the predicted result for each county, I need to get the random
> effect (residuals) for each county as well as their variances. However, the
> ranef command only produce the random effect without giving the confidence
> interval / variance. I would be glad if you can let me know how to produce
> the confidence interval of these random effects or predicted coefficients
> in R.
>  > ranef(fit2)
> $county
>       (Intercept)
> 4001  0.612613842
> 4003  0.066916182
> 4005  1.293875058
> My last question is whether the random residuals/coefficients generated in
> R are empirically bayesian estimates, as in HLM.
>
> Thanks a lot
>
> Lijun
>
>
>
>
>
> Lijun Chen
> Senior Research Specialist
> Chapin Hall at the University of Chicago
> 1313 E. 60th St., Chicago, IL 60637
> V: 773-256-5140 F: 773-256-5340
> http://www.chapinhall.org
>

	[[alternative HTML version deleted]]

