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

From: Douglas Bates <bates_at_stat.wisc.edu>
Date: Sat, 08 Mar 2008 08:07:29 -0600

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. :-)

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
>
> Family: poisson(log link)
>
> 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
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



R-help_at_r-project.org 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. Received on Sat 08 Mar 2008 - 14:13:43 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sat 08 Mar 2008 - 17:30:20 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive