# Re: [R] analysing mixed effects/poisson/correlated data

From: Manuel Morales <Manuel.A.Morales_at_williams.edu>
Date: Sat, 08 Mar 2008 11:00:38 -0500

On Sat, 2008-03-08 at 08:07 -0600, Douglas Bates wrote:
> On Sat, Mar 8, 2008 at 2:57 AM, Alexandra Bremner
> <alexandra.bremner_at_uwa.edu.au> wrote:
> > I am attempting to model data with the following variables:
>
> > timepoint - n=48, monthly over 4 years
> > hospital - n=3
> > opsn1 - no of outcomes
> > total.patients
> > skillmixpc - skill mix percentage
> > nurse.hours.per.day
>
> > Aims
> > To determine if skillmix affects rate (i.e. no.of.outcomes/total.patients).
> > To determine if nurse.hours.per.day affects rate.
> > To determine if rates vary between hospitals.
>
> > There is first order autoregression in the data. I have attempted to use the lmer function (and lmer2) with correlation structure and weights:
>
> > test1 <-lmer(opsn1~timepoint+as.factor(hospital)+ skillmixpc + nursehrsperpatday +(timepoint|hospital) +offset(log(totalpats)),family=poisson, data=opsn.totals)
> > test2 <-lmer(opsn1~timepoint+as.factor(hospital)+ skillmixpc + nursehrsperpatday +(timepoint|hospital)+offset(log(totalpats)),family=poisson, data=opsn.totals, correlation=corAR1(form=~1|hospital))
> > test3 <-lmer(opsn1~timepoint+as.factor(hospital)+ skillmixpc + nursehrsperpatday +(timepoint|hospital)+offset(log(totalpats)),family=poisson, data=opsn.totals, correlation=corAR1(form=~1|hospital),weights=varIdent(form=~1|hospital))
>
> You are mixing arguments for lme or nlme into a call to lmer. Because
> the weigths argument doesn't have the form required by lmer you get an
> error message. The effect of the correlation argument is more subtle
> - because lmer has ... in the argument list your correlation
> specification is absorbed without an error message but it has no
> effect.
>
> The lmer documentation doesn't say that you can use the forms of the
> correlation and weights arguments from the lme function, although you
> are not the first person to decide that it should. :-)

The documentation for weights in lmer references lm. It looks to me like the weights argument for lm requires a vector of weights a priori - does that mean lmer cannot estimate heteroscedasticity like lme can?

> There are theoretical problems with trying to specify a separate
> correlation argument in a generalized linear mixed model. In a GLMM
> the conditional variance of the response (i.e. the variance of Y given
> a value of B, the random effects) depends on the conditional mean so
> it is more difficult to decide what would be the effect if a
> correlation structure or a non-constant weighting structure were
> overlaid on it.
>
> >
> >
> > Test1 & test2 give the same output (below). Does this mean that the correlation structure is not being used?
> >
> > Test3 produces the following error message (I notice there are others who have had problems with weights).
> >
> > Error in model.frame(formula, rownames, variables, varnames, extras, extranames, :
> >
> > variable lengths differ (found for '(weights)')
> >
> >
> >
> > > summary(test1)
> >
> > Generalized linear mixed model fit using Laplace
> >
> > Formula: opsn1 ~ timepoint + as.factor(hospital) + skillmixpc + nursehrsperpatday + (timepoint | hospital) + offset(log(totalpats))
> >
> > Data: opsn.totals
> >
> >
> > AIC BIC logLik deviance
> >
> > 196.2 223.0 -89.12 178.2
> >
> > Random effects:
> >
> > Groups Name Variance Std.Dev. Corr
> >
> > hospital (Intercept) 3.9993e-03 6.3240e-02
> >
> > timepoint 5.0000e-10 2.2361e-05 0.000
> >
> > number of obs: 144, groups: hospital, 3
> >
> >
> >
> > Estimated scale (compare to 1 ) 1.057574
> >
> >
> >
> > Fixed effects:
> >
> > Estimate Std. Error z value Pr(>|z|)
> >
> > (Intercept) -2.784857 1.437283 -1.9376 0.0527 .
> >
> > timepoint -0.002806 0.002709 -1.0358 0.3003
> >
> > as.factor(hospital)2 -0.030277 0.120896 -0.2504 0.8022
> >
> > as.factor(hospital)3 -0.349763 0.171950 -2.0341 0.0419 *
> >
> > skillmixpc -0.026856 0.015671 -1.7137 0.0866 .
> >
> > nursehrsperpatday -0.025376 0.043556 -0.5826 0.5602
> >
> > ---
> >
> > Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
> >
> >
> >
> > Correlation of Fixed Effects:
> >
> > (Intr) timpnt as.()2 as.()3 skllmx
> >
> > timepoint -0.606
> >
> > as.fctr(h)2 -0.382 0.132
> >
> > as.fctr(h)3 -0.734 0.453 0.526
> >
> > skillmixpc -0.996 0.596 0.335 0.715
> >
> > nrshrsprptd 0.001 -0.297 0.204 -0.130 -0.056
> >
> >
> >
> > Can anyone suggest another way that I can do this analysis please? Or a work around for this method?
> >
> > Any help will be gratefully received as I have to model 48 different opsn outcomes.
> >
> >
> >
> > Alex
> >
> >
> >
> >
> > ______________________________________________
> > R-help_at_r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > and provide commented, minimal, self-contained, reproducible code.
> >
>
> ______________________________________________
> R-help_at_r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> and provide commented, minimal, self-contained, reproducible code.

```--
http://mutualism.williams.edu

______________________________________________
R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help