On Sat, 2008-03-08 at 08:07 -0600, Douglas Bates wrote: > On Sat, Mar 8, 2008 at 2:57 AM, Alexandra Bremner > <[EMAIL PROTECTED]> 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@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@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
______________________________________________ R-help@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.
