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