Dear all
I am trying to do multiple comparisons in a lme model using s-plus
6.2. Below is the code for the model, where income (4 levels),
leisure (4 levels) and hs (3 levels) are categorical fixed effect
predictor variables, and B some continuous response variable.
> m3<-lme(fixed = B ~ income + leisure + hs, data = d, random = ~ 1
| id, correlation = corSymm(), na.action = na.omit)
> summary(m3)
Linear mixed-effects model fit by REML
Data: d
AIC BIC logLik
1800.773 1848.027 -886.3867
Random effects:
Formula: ~ 1 | id
(Intercept) Residual
StdDev: 9.20643 12.15391
Correlation Structure: General
Formula: ~ 1 | id
Parameter estimate(s):
Correlation:
1 2
2 0.162
3 -0.496 0.277
Fixed effects: B ~ income + leisure + hs
Value Std.Error DF t-value p-value
(Intercept) 75.22513 3.683161 142 20.42407 <.0001
income2 -3.00479 3.706030 142 -0.81078 0.4188
income3 -0.34330 3.508361 142 -0.09785 0.9222
income4 4.76078 3.576659 142 1.33107 0.1853
leisure2 7.95498 3.038379 142 2.61817 0.0098
leisure3 -3.86922 3.284164 142 -1.17814 0.2407
leisure4 -0.85617 3.877200 142 -0.22082 0.8255
hs2 -19.38267 2.288641 142 -8.46907 <.0001
hs3 -49.44063 3.110348 142 -15.89553 <.0001
Number of Observations: 225
Number of Groups: 75
> anova(m3)
numDF denDF F-value p-value
(Intercept) 1 142 1761.043 <.0001
income 3 142 63.039 <.0001
leisure 3 142 22.158 <.0001
hs 2 142 150.489 <.0001
THE ANOVA OUTPUT INDICATES THAT INCOME SIGNIFICANTLY CONTRIBUTES TO
THE MODEL, AND I WOULD LIKE TO DO MULTIPLE COMPARISONS AMONG THE FOUR
LEVELS OF THE INCOME VARIABLE
> contrasts(d$income) TREATMENT CONTRASTS USED
2 3 4
1 0 0 0
2 1 0 0
3 0 1 0
4 0 0 1
> m3<-lme(fixed = B ~ income-1 + leisure + hs, data = d, random = ~
1 | id, correlation = corSymm(), na.action = na.omit)
IN ORDER TO COMPARE THE DIFFERENCE IN RESPONSE AMONG THE FOUR LEVELS
OF INCOME I HAVE CHANGED THE CODING TO THE "CELL MEANS"
PARAMETRIZATION AND ADDED THE TERM -1.
> summary(m3)
Linear mixed-effects model fit by REML
Data: d
AIC BIC logLik
1800.773 1848.027 -886.3867
Random effects:
Formula: ~ 1 | id
(Intercept) Residual
StdDev: 9.206424 12.15392
Correlation Structure: General
Formula: ~ 1 | id
Parameter estimate(s):
Correlation:
1 2
2 0.162
3 -0.496 0.277
Fixed effects: B ~ income - 1 + leisure + hs
Value Std.Error DF t-value p-value
income1 75.22513 3.683161 142 20.42407 <.0001
income2 72.22034 4.068986 142 17.74898 <.0001
income3 74.88183 3.683089 142 20.33126 <.0001
income4 79.98592 3.332672 142 24.00054 <.0001
leisure2 7.95498 3.038379 142 2.61817 0.0098
leisure3 -3.86922 3.284164 142 -1.17814 0.2407
leisure4 -0.85617 3.877200 142 -0.22082 0.8255
hs2 -19.38267 2.288642 142 -8.46907 <.0001
hs3 -49.44063 3.110348 142 -15.89553 <.0001
Number of Observations: 225
Number of Groups: 75
> fe<-fixef(m3)[1:4]
> vc<-m3$varFix[1:4,1:4]
> fe
income1 income2 income3 income4
75.22513 72.22034 74.88183 79.98592
> vc
income1 income2 income3 income4
income1 13.565678 8.193832 7.411114 5.939948
income2 8.193832 16.556648 6.739330 7.166021
income3 7.411114 6.739330 13.565143 8.680778
income4 5.939948 7.166021 8.680778 11.106703
>
> multicomp(fe,vc)
EXTRACT FIXED EFFECT ESTIMATES FOR INCOME AND THE CORRESPONDING VCOV
MATRIX, AND THE APPLY MULTICOMP
95 % simultaneous confidence intervals for specified
linear combinations, by the Sidak method
critical point: 2.636
response variable:
intervals excluding 0 are flagged by '****'
Estimate Std.Error Lower Bound Upper Bound
income1-income2 3.000 3.71 -6.76 12.80
income1-income3 0.343 3.51 -8.90 9.59
income1-income4 -4.760 3.58 -14.20 4.67
income2-income3 -2.660 4.08 -13.40 8.09
income2-income4 -7.770 3.65 -17.40 1.86
income3-income4 -5.100 2.70 -12.20 2.02
NONE OF THE DIFFERENCES IS NOW SIGNIFICANT. IS THE RIGHT WAY TO GO
AND THANKS A LOT FOR COMMENTS REGARDING THIS APPROACH.
PD Dr. med. Pedram Sendi, DSc
Health Economics Research Unit
Institute for Clinical Epidemiology
& Division of Infectious Diseases
Basel University Hospital
Hebelstrasse 10, 3rd Floor
CH-4031 Basel
Switzerland
Phone: +41 61 265 31 02
Fax: +41 61 265 31 09
E-Mail: psendi@uhbs.ch
www.bice.ch
|