[R] analysing mixed effects/poisson/correlated data

[Manuel Morales]

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