Skip to content

analysing mixed effects/poisson/correlated data

4 messages · Alexandra Bremner, Manuel Morales, Douglas Bates

Alexandra Bremner
# 
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))

 

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)')
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
Douglas Bates
# 
On Sat, Mar 8, 2008 at 2:57 AM, Alexandra Bremner
<alexandra.bremner at uwa.edu.au> wrote:
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.
Manuel Morales
#
On Sat, 2008-03-08 at 08:07 -0600, Douglas Bates wrote:
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?
1 day later
Douglas Bates
# 
On Sat, Mar 8, 2008 at 11:00 AM, Manuel Morales
<Manuel.A.Morales at williams.edu> wrote:

            

      
      
      
    

        
Yes.  lmer behaves the way that the documentation says it does.

As I indicated in my answer to Alexandra, it is not as easy as it may
seem to define what "modeling heteroscedasticity" should mean for
models like generalized linear mixed models.  Even for linear mixed
models incorporating correlation structures in addition to the
correlation in the marginal distribution of the response induced by
the random effects will frequently result in an overspecified model.

The purpose of lmer is to allow general specification of mixed-effects
models.  If someone wants to wrap the underlying representation in a
more general model specification they are free to do so.  It's open
source.