Hello R users,
I'm trying to take into acount the center effect in a clinical study
comparing 3 treatments (data from the book Applied Mixed Models in Medicine
Statistics by Helen Brown and Robin Prescott):
- dbp: diastolic blood pressure at 8 weeks
- dbp0: diastolic blood pressure at inclusion
- treat: three treatments A, B and C
- centre: centres of the study
First, I have fitted a marginal model with gls (nlme package), with a
compound symmetry structure:
require(nlme)
gls1 <- gls(dbp ~ factor(centre) + treat + dbp0, data = tension1,
correlation = corCompSymm(form = ~ 1 | centre))
summary(gls1)
anova(gls1)
Denom. DF: 224
numDF F-value p-value
(Intercept) 1 33349.45 <.0001
factor(centre) 28 2.38 0.0003
treat 2 1.79 0.1700
dbp0 1 2.61 0.1075
Second, I tried to fit a random effect model with lme:
lme1 <- lme(dbp ~ factor(centre) + treat + dbp0, random = ~ 1 | centre, data
= tension1)
anova(lme1)
numDF denDF F-value p-value
(Intercept) 1 224 3753.228 <.0001
factor(centre) 28 0 0.335 NaN
treat 2 224 1.786 0.1700
dbp0 1 224 2.612 0.1075
Message d'avis :
In pf(q, df1, df2, lower.tail, log.p) : production de NaN
I thought this two models equivalent, and I don't understand why the fixed
center effect produced NaN in lme model.
I also tried with lmer (from lme4 package):
require(lme4)
lmer1 <- lmer(dbp ~ factor(centre) + treat + dbp0 + (1 | centre), data =
tension1)
anova(lmer1)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
factor(centre) 28 1214.42 43.372 0.7026
treat 2 220.52 110.262 1.7862
dbp0 1 161.23 161.231 2.6118
F values for centre effect are very different in the 3 models.
It is surely a naive question, but does someone could to me explain what
happens here?
Thank you very much in advance.
Here the data :
tension1 <-
structure(list(patient = c(1L, 3L, 4L, 5L, 7L, 8L, 9L, 10L, 11L,
13L, 14L, 15L, 18L, 19L, 23L, 24L, 25L, 28L, 30L, 31L, 32L, 34L,
35L, 36L, 37L, 38L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 50L, 52L,
53L, 54L, 55L, 56L, 57L, 58L, 60L, 63L, 64L, 70L, 71L, 72L, 73L,
74L, 80L, 81L, 82L, 84L, 92L, 93L, 94L, 95L, 96L, 97L, 98L, 99L,
100L, 101L, 102L, 103L, 104L, 105L, 106L, 107L, 108L, 109L, 110L,
111L, 112L, 113L, 114L, 116L, 117L, 118L, 119L, 120L, 122L, 124L,
125L, 127L, 128L, 129L, 130L, 131L, 132L, 133L, 134L, 135L, 136L,
137L, 139L, 140L, 141L, 142L, 143L, 144L, 145L, 146L, 147L, 150L,
157L, 158L, 160L, 161L, 162L, 165L, 166L, 167L, 169L, 170L, 171L,
172L, 173L, 174L, 175L, 176L, 177L, 178L, 180L, 181L, 182L, 183L,
185L, 186L, 187L, 189L, 190L, 191L, 192L, 199L, 200L, 201L, 202L,
203L, 204L, 205L, 206L, 207L, 208L, 212L, 213L, 214L, 216L, 217L,
218L, 219L, 220L, 221L, 223L, 224L, 225L, 226L, 227L, 228L, 229L,
230L, 237L, 241L, 242L, 243L, 244L, 245L, 246L, 247L, 248L, 249L,
250L, 251L, 252L, 255L, 256L, 258L, 260L, 261L, 262L, 263L, 264L,
265L, 266L, 267L, 268L, 269L, 270L, 271L, 272L, 273L, 274L, 275L,
279L, 280L, 281L, 282L, 283L, 284L, 285L, 286L, 287L, 290L, 291L,
292L, 293L, 294L, 301L, 302L, 303L, 304L, 305L, 306L, 307L, 308L,
309L, 310L, 311L, 312L, 319L, 320L, 322L, 323L, 324L, 325L, 326L,
327L, 328L, 329L, 330L, 332L, 333L, 334L, 335L, 337L, 338L, 339L,
340L, 341L, 342L, 343L, 344L, 345L, 346L, 347L, 348L, 350L, 352L,
355L, 357L, 361L, 362L, 363L, 364L, 365L, 366L), centre = c(29L,
5L, 5L, 29L, 3L, 3L, 3L, 3L, 3L, 36L, 36L, 36L, 36L, 5L, 5L,
5L, 15L, 15L, 15L, 4L, 4L, 4L, 4L, 4L, 23L, 23L, 7L, 7L, 7L,
7L, 7L, 7L, 30L, 30L, 30L, 30L, 30L, 18L, 18L, 18L, 18L, 18L,
24L, 23L, 35L, 35L, 35L, 27L, 27L, 1L, 1L, 1L, 1L, 26L, 26L,
26L, 26L, 26L, 14L, 14L, 14L, 14L, 14L, 14L, 11L, 11L, 11L, 11L,
11L, 11L, 7L, 7L, 7L, 7L, 7L, 7L, 6L, 6L, 6L, 6L, 6L, 32L, 32L,
32L, 25L, 25L, 25L, 25L, 25L, 25L, 31L, 31L, 31L, 31L, 31L, 14L,
14L, 14L, 14L, 14L, 14L, 37L, 37L, 37L, 37L, 31L, 31L, 31L, 31L,
31L, 13L, 13L, 13L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 9L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 31L, 31L, 31L, 31L, 31L,
31L, 8L, 8L, 8L, 8L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L,
12L, 36L, 36L, 36L, 36L, 36L, 36L, 40L, 40L, 41L, 14L, 14L, 14L,
14L, 14L, 14L, 4L, 4L, 4L, 4L, 4L, 4L, 36L, 36L, 36L, 1L, 1L,
1L, 1L, 1L, 7L, 7L, 7L, 7L, 7L, 7L, 1L, 1L, 1L, 1L, 1L, 5L, 5L,
5L, 5L, 25L, 26L, 26L, 26L, 26L, 12L, 12L, 12L, 12L, 12L, 15L,
15L, 15L, 15L, 15L, 15L, 1L, 1L, 1L, 1L, 1L, 1L, 36L, 36L, 36L,
36L, 36L, 31L, 31L, 31L, 31L, 31L, 31L, 1L, 1L, 1L, 1L, 14L,
14L, 14L, 14L, 14L, 14L, 31L, 31L, 31L, 31L, 31L, 31L, 1L, 1L,
36L, 36L, 31L, 31L, 31L, 31L, 31L, 31L), treat = structure(c(3L,
2L, 1L, 1L, 1L, 2L, 2L, 1L, 3L, 2L, 1L, 3L, 1L, 2L, 1L, 3L, 3L,
2L, 2L, 2L, 3L, 2L, 1L, 3L, 3L, 3L, 2L, 1L, 1L, 3L, 3L, 2L, 2L,
3L, 3L, 2L, 1L, 2L, 1L, 3L, 3L, 2L, 3L, 1L, 3L, 2L, 1L, 2L, 3L,
2L, 3L, 1L, 3L, 3L, 1L, 1L, 2L, 3L, 2L, 1L, 1L, 3L, 2L, 3L, 3L,
2L, 3L, 2L, 1L, 1L, 3L, 2L, 1L, 3L, 1L, 2L, 3L, 1L, 3L, 2L, 1L,
1L, 2L, 3L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 2L, 2L, 1L, 3L, 1L, 1L,
3L, 3L, 2L, 2L, 1L, 3L, 2L, 3L, 2L, 2L, 3L, 3L, 1L, 2L, 3L, 3L,
3L, 1L, 2L, 2L, 1L, 3L, 2L, 3L, 2L, 1L, 3L, 2L, 3L, 1L, 3L, 2L,
1L, 2L, 3L, 3L, 2L, 2L, 2L, 1L, 3L, 1L, 3L, 3L, 1L, 1L, 3L, 3L,
3L, 1L, 2L, 2L, 1L, 3L, 3L, 1L, 3L, 2L, 3L, 1L, 2L, 1L, 2L, 1L,
3L, 2L, 1L, 3L, 3L, 1L, 2L, 3L, 2L, 1L, 2L, 1L, 3L, 1L, 3L, 1L,
2L, 3L, 1L, 1L, 3L, 1L, 2L, 1L, 3L, 3L, 2L, 2L, 2L, 3L, 1L, 1L,
1L, 2L, 2L, 3L, 1L, 3L, 2L, 2L, 1L, 3L, 1L, 1L, 2L, 2L, 2L, 3L,
1L, 1L, 3L, 2L, 2L, 3L, 3L, 1L, 2L, 1L, 2L, 2L, 1L, 3L, 1L, 1L,
1L, 3L, 3L, 2L, 2L, 3L, 1L, 2L, 2L, 1L, 1L, 3L, 3L, 2L, 2L, 1L,
2L, 1L, 3L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 1L, 1L, 2L, 3L, 3L), .Label = c("A",
"B", "C"), class = "factor"), dbp = c(86L, 109L, 87L, 85L, 100L,
80L, 90L, 100L, 94L, 80L, 85L, 80L, 100L, 102L, 98L, 84L, 83L,
87L, 110L, 94L, 90L, 86L, 86L, 100L, 92L, 94L, 92L, 88L, 106L,
92L, 70L, 81L, 100L, 90L, 82L, 92L, 100L, 90L, 96L, 90L, 100L,
95L, 88L, 90L, 96L, 90L, 74L, 92L, 92L, 100L, 72L, 116L, 90L,
94L, 86L, 86L, 92L, 94L, 90L, 85L, 90L, 90L, 100L, 80L, 88L,
82L, 100L, 96L, 100L, 98L, 106L, 80L, 102L, 92L, 98L, 80L, 98L,
90L, 86L, 92L, 98L, 80L, 100L, 80L, 95L, 98L, 85L, 90L, 93L,
93L, 70L, 80L, 80L, 80L, 80L, 88L, 84L, 90L, 80L, 96L, 90L, 86L,
80L, 84L, 86L, 85L, 80L, 85L, 90L, 85L, 98L, 110L, 106L, 90L,
80L, 85L, 86L, 94L, 95L, 80L, 85L, 80L, 80L, 96L, 98L, 92L, 88L,
96L, 98L, 78L, 102L, 88L, 94L, 74L, 90L, 80L, 90L, 80L, 85L,
85L, 105L, 100L, 95L, 100L, 80L, 70L, 96L, 84L, 88L, 105L, 100L,
85L, 105L, 85L, 92L, 105L, 100L, 85L, 80L, 98L, 96L, 90L, 88L,
90L, 90L, 80L, 90L, 95L, 100L, 82L, 84L, 82L, 84L, 78L, 90L,
70L, 88L, 98L, 78L, 98L, 88L, 82L, 76L, 86L, 100L, 78L, 90L,
86L, 92L, 92L, 92L, 80L, 82L, 104L, 80L, 86L, 84L, 92L, 100L,
96L, 96L, 92L, 88L, 85L, 110L, 100L, 85L, 102L, 79L, 94L, 89L,
89L, 84L, 90L, 80L, 102L, 90L, 102L, 104L, 90L, 90L, 102L, 90L,
90L, 90L, 90L, 80L, 90L, 90L, 90L, 92L, 92L, 108L, 90L, 80L,
80L, 90L, 95L, 80L, 90L, 80L, 90L, 70L, 85L, 90L, 80L, 92L, 98L,
82L, 85L, 80L, 85L, 76L, 80L, 84L, 80L), dbp0 = c(97L, 117L,
100L, 105L, 114L, 105L, 100L, 102L, 105L, 100L, 100L, 100L, 100L,
100L, 100L, 98L, 101L, 98L, 102L, 104L, 102L, 96L, 98L, 102L,
104L, 104L, 100L, 105L, 100L, 100L, 105L, 101L, 100L, 112L, 100L,
102L, 114L, 105L, 108L, 105L, 110L, 100L, 102L, 102L, 100L, 96L,
104L, 114L, 96L, 100L, 100L, 100L, 100L, 102L, 104L, 102L, 102L,
106L, 105L, 110L, 100L, 110L, 105L, 100L, 105L, 100L, 95L, 102L,
98L, 105L, 115L, 95L, 104L, 115L, 105L, 95L, 104L, 96L, 114L,
102L, 108L, 105L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L,
100L, 100L, 105L, 100L, 100L, 110L, 105L, 110L, 104L, 100L, 100L,
96L, 92L, 104L, 96L, 105L, 100L, 100L, 100L, 95L, 112L, 114L,
114L, 102L, 110L, 110L, 104L, 100L, 100L, 100L, 100L, 100L, 100L,
104L, 100L, 108L, 100L, 98L, 98L, 110L, 100L, 100L, 104L, 100L,
105L, 100L, 105L, 105L, 100L, 105L, 105L, 100L, 100L, 100L, 100L,
104L, 110L, 108L, 95L, 110L, 105L, 95L, 98L, 105L, 95L, 105L,
95L, 100L, 102L, 96L, 108L, 104L, 105L, 102L, 105L, 110L, 100L,
100L, 102L, 100L, 98L, 98L, 96L, 100L, 100L, 100L, 110L, 100L,
100L, 104L, 100L, 106L, 105L, 100L, 95L, 102L, 95L, 105L, 102L,
102L, 106L, 104L, 108L, 108L, 104L, 102L, 102L, 112L, 110L, 104L,
106L, 96L, 95L, 110L, 105L, 105L, 103L, 105L, 103L, 99L, 102L,
103L, 98L, 100L, 100L, 98L, 96L, 102L, 110L, 95L, 110L, 100L,
105L, 105L, 105L, 105L, 105L, 110L, 105L, 105L, 100L, 104L, 98L,
108L, 100L, 100L, 110L, 105L, 100L, 105L, 110L, 95L, 100L, 105L,
110L, 100L, 98L, 104L, 102L, 98L, 105L, 100L, 110L, 110L, 100L,
100L)), .Names = c("patient", "centre", "treat", "dbp", "dbp0"
), row.names = c(NA, 256L), class = "data.frame")
[[alternative HTML version deleted]]
______________________________________________
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.
