Multilevel modeling with count variables

[dadrivr]

I am using a multilevel modeling approach to model change in a person's symptom score over time (i.e., longitudinal individual growth models).  I have been using the lme function in the multilevel package for the analyses, but my problem is that my outcome (symptoms) and one of my predictors (events) are count data, and are non-normal.  Do you have any suggestions for how to deal with them?  Are there poisson-regression or similar techniques that can be applied to multilevel modeling (lme, nlme, etc.)?  Some of my predictors are normal, however, so the technique should be able to accommodate normally distributed data, as well.  I have tried square root transformations, and it appears that the data become more normal, but at the expense of some of the consistency in my findings modeling the raw data.  Your input would be very helpful and greatly appreciated.  Thanks so much!

[dadrivr]

By the way, my concern with lmer and glmer is that they don't produce p-values, and the techniques used to approximate the p-values with those functions (pvals.fnc, HPDinterval, mcmcsamp, etc.) only apply to Gaussian distributions.  Given that I would likely be working with quasi-poisson distributions, is there a good alternative mixed effects modeling approach that would output significance values?  Thanks!

[Emmanuel Charpentier]

Your request might find better answers on the R-SIG-mixed-models
list ...

Anyway, some quick thoughts :

Le vendredi 26 mars 2010 à 15:20 -0800, dadrivr a écrit :
> By the way, my concern with lmer and glmer is that they don't produce
> p-values,

The argumentation of D. Bates is convincing ... A large thread about
these issues exist on this (R-help) list archive.

I won't get in the (real) debate about the (dubious) value of hypothesis
testing, and I'm aware that there are domains where journal editors
*require* p-values. Sigh... This insistence seems to be weakening,
however.

>            and the techniques used to approximate the p-values with those
> functions (pvals.fnc, HPDinterval, mcmcsamp, etc.) only apply to Gaussian
> distributions.

(Probably ?) true for pvals.fnc, false for HPDinterval and mcmcsamp. You
can always generate a "sufficient" number of estimates and assess a
confidence interval and a p-value from such samples. See any good text
on use of simulation in statistics...

The fly in the ointment is that current versions of lmer seem to have
problems with mcmcsamp(). Bootstrap comes to mind, but might be
non-trivial to do efficiently, especially in a hierarchical/multilevel
context.

> Given that I would likely be working with quasi-poisson
> distributions, is there a good alternative mixed effects modeling approach

glmer, Bayesian modeling through BUGS (whose "probabilities"
interpretation is (quite) different from p-values).

> that would output significance values?

No. You'd have to compute them yourself ... if you're able to derive the
(possibly asymptotic) distribution of your test statistic under the
hypothesis you aim to refute. This nut might be harder to crack than it
appears...

Now, to come back to your original question, you might try 1) to use
glm/glmer to model your dependent variable as a (quasi-)Poisson
variable, and 2) use log transformations of the independent cout
variables ; a more general solution is given by the Box-Cox
transformation family : see the relevant function in MASS, which also
offers the logtrans family. ISTR that John Fox's car package offers a
function aiming at finding the optimal simultaneous transformations of a
dependent variable and its predictors. Other packages I'm not aware of
might offer other solutions ; in particular, I'd recommend to peek at
Frank Harrell's Hmisc and rms packages documentation, whose wealth I did
not yet seriously assess...

Bayesian modeling with BUGS would also allow you to try to fit any model
you might wish to test (provided that it can be written as a directed
acyclic graph and that the  distributions of you variables are either
from a "standard" family available in BUGS or that you are able to
express the (log-)density of the non-standard distribution you wish to
use). But, again, no p-values in sight. Would you settle for Bayes
factors between two models ? or DIC comparisons ?

HTH,

                                        Emmanuel Charpentier

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[Corey Sparks]

have you tried using glmer?
If your dependent variable is poisson distributed, you can try something like
fit<-glmer(y~x+(1|group), family=poisson)

and if you have differential exposure, you can do

fit<-glmer(y~offset(log(exposure))+x+(1|group), family=poisson)

Is this what you are asking?
With regard to the t-statistics generated from lmer/glmer, you can get p-values by using dt(), or look at your confidence intervals for the parameters.
Does this help?

Corey

[Corey Sparks]

Whoops, sorry that's pt(), not dt()
Thanks Dennis!

[dadrivr]

Thanks everyone for the helpful ideas.  It appears that this will be more difficult than I thought.  I don't necessary have an inclination toward p-values, but many journals certainly do.  I would be willing to try to calculate the confidence intervals around the estimates, but I haven't gotten any functions to work when applied to glmer & (quasi-)poisson distributions.  Alternatively, it appears that I could look at using model comparison, Bayesian (BUGS), or other software (SAS/Stata).  Model comparison seems a little inadequate for reporting on fixed effects.  I'm not yet very familiar with Bayesian approaches, so I would prefer to avoid that direction -- at least for now.  As for other software, I may very well end up using SAS for confirmatory analysis (and potentially reporting), but for exploratory analysis, I'd prefer to stick with R to continue to learn about the ever-advancing approach to analyzing statistics.  I'm still a novice with R, but I recognize that its flexibility and user community will make it a strong statistics package now and into the future.  Thanks again guys!