From: Manuel Morales <Manuel.A.Morales_at_williams.edu>

Date: Sat, 08 Mar 2008 11:00:38 -0500

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
**> >
**> > 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.
*

-- http://mutualism.williams.eduReceived on Sat 08 Mar 2008 - 16:42:18 GMT______________________________________________ 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.

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 Mon 10 Mar 2008 - 16:30:22 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.
*