Dear professors and collegues
I need to perform a analysis of dates from a nested experimental design.
From
"Bioestatical Analysis" of Zar
"Bimetry of Sokal" & Rohlf
"Design and Analysis of Experiments" of Montgomery
I have:
Sum (mean(x)_i - mean(x)_T)2 / (a-1) -> var(epsilon) + n sigma2_B + n b
(sum alfa_i)2 / (a-1)
Sum (mean(x)_ij - mean(x)_i)2 / (ba-a) -> var(epsilon) + n sigma2_B
Sum (x_ijl - mean(x)_ij)2 / (abn-ab) -> var(epsilon)
Dates for the execution of a example:
#> Mesures <- as.data.frame(matrix(rnorm(45*1, mean=10, sd=1), ncol=1))
#> colnames(Mesures) <- "VR"
#> Mesures$trat <- as.factor(rep(1:5,rep(9,5)))
#> Mesures$patient <- as.factor(rep(1:15,rep(3,15)))
The Analysis of Variance could be:
#> AnovaModel.1 <-aov(VR ~ trat + trat/patient, data=Mesures)
#> summary(AnovaModel.1)
Df Sum Sq Mean Sq F value Pr(>F)
trat 4 1.580 0.3950 0.751 0.565
trat:patient 10 4.811 0.4811 0.915 0.533
Residuals 30 15.778 0.5259
But this methods causes pseudoreplication because the F is the ratio of
"MS of trat" and "MS of res". However the variability of
"trat" is at
least as big as the variability of "patient".
The correct analysis suggested by ?aov is:
#> AnovaModel.2 <- aov(VR ~ trat + Error(patient),data=Mesures)
#> summary(AnovaModel.2)
Error: patient
Df Sum Sq Mean Sq F value Pr(>F)
trat 4 1.580 0.3950 0.821 0.541
Residuals 10 4.811 0.4811
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 30 15.78 0.5259
Althougth the Anova Model object it isn't a object of class "aov"
but a
object of class"aovlist":
#> attr(AnovaModel.2,"class")
[1] "aovlist" "listof"
This class isn't suitable for the glht function neither TukeyHSD function.
I can extract a object of class "aov" through:
#> ls.str(pat="AnovaModel.2")
AnovaModel.2 : List of 3
$ (Intercept):List of 9
$ patient :List of 9
$ Within :List of 6
#> AnovaModel.3<-AnovaModel.2$patient
#> attr(AnovaModel.3,"class")
[1] "aov" "lm"
This model is suitable for confidence intervals:
#> confint(AnovaModel.3)
2.5 % 97.5 %
trat2 -0.9107959 0.5462655
trat3 -0.9848351 0.4722263
trat4 -1.2997235 0.1573379
trat5 -1.0623118 0.3947496
But it is useless for multiple comparisons:
#> Pairs <- glht(AnovaModel.3, linfct = mcp(Tratamiento =
"Tukey"))
Error in model.matrix.default(model, data = structure(list(method = lm),
.Names = "method"), :
model structure is invalid
Error in factor_contrasts(model) :
no 'model.matrix' method for 'model' found!
#> TukeyHSD(AnovaModel.3, "Tratamiento")
Error in model.tables.aov(x, "means") :
this fit does not inherit from "lm"
Another way of analysis is through lme4 package:
#> AnovaModel.3 <- lmer(VR ~ trat+(1|trat/patient), data=Mesures)
But the anova test isn't the correct ratio above shown:
#> anova(AnovaModel.3)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
trat 4 0.43112 0.10778 0.2094
Other options for lmer:
#> anova(AnovaModel.4)
Analysis of Variance Table
Df Sum Sq Mean Sq F value
trat 4 1.5801 0.39501 0.7675
Neither is the F ratio searched and above shown:
#> summary(AnovaModel.2)
Error: patient
Df Sum Sq Mean Sq F value Pr(>F)
trat 4 1.580 0.3950 0.821 0.541
Residuals 10 4.811 0.4811
I can also take the mean of treatment, but if a measures of some patient
fails, and the design becomes unbalanced, the results would be very
difficult.
Can anyone give me the right model and multiple comparisons?
Thanks very much for the help and please excuse for my bad english.
--
_____---^---_____
Univ. de Extremadura
Dept. Matemáticas.
Despacho B29
Tf: + 34 924 289 300
Ext. 86823
[[alternative HTML version deleted]]
José Trujillo Carmona
2012-Feb-06 13:24 UTC
[R] multiple comparisons in nested design (error correction)
The last model is: AnovaModel.4 <- lmer(VR ~ trat+(1|patient:trat), data=Mesures) Sorry. El 06/02/12 13:43, José Trujillo Carmona escribió:> Dear professors and collegues > > I need to perform a analysis of dates from a nested experimental design. > > From > > "Bioestatical Analysis" of Zar > "Bimetry of Sokal" & Rohlf > "Design and Analysis of Experiments" of Montgomery > > I have: > > Sum (mean(x)_i - mean(x)_T)2 / (a-1) -> var(epsilon) + n sigma2_B + n > b (sum alfa_i)2 / (a-1) > Sum (mean(x)_ij - mean(x)_i)2 / (ba-a) -> var(epsilon) + n sigma2_B > Sum (x_ijl - mean(x)_ij)2 / (abn-ab) -> var(epsilon) > > Dates for the execution of a example: > > #> Mesures <- as.data.frame(matrix(rnorm(45*1, mean=10, sd=1), ncol=1)) > #> colnames(Mesures) <- "VR" > #> Mesures$trat <- as.factor(rep(1:5,rep(9,5))) > #> Mesures$patient <- as.factor(rep(1:15,rep(3,15))) > > The Analysis of Variance could be: > > #> AnovaModel.1 <-aov(VR ~ trat + trat/patient, data=Mesures) > #> summary(AnovaModel.1) > Df Sum Sq Mean Sq F value Pr(>F) > trat 4 1.580 0.3950 0.751 0.565 > trat:patient 10 4.811 0.4811 0.915 0.533 > Residuals 30 15.778 0.5259 > > But this methods causes pseudoreplication because the F is the ratio > of "MS of trat" and "MS of res". However the variability of "trat" is > at least as big as the variability of "patient". > > The correct analysis suggested by ?aov is: > > #> AnovaModel.2 <- aov(VR ~ trat + Error(patient),data=Mesures) > #> summary(AnovaModel.2) > > Error: patient > Df Sum Sq Mean Sq F value Pr(>F) > trat 4 1.580 0.3950 0.821 0.541 > Residuals 10 4.811 0.4811 > > Error: Within > Df Sum Sq Mean Sq F value Pr(>F) > Residuals 30 15.78 0.5259 > > Althougth the Anova Model object it isn't a object of class "aov" but > a object of class"aovlist": > > #> attr(AnovaModel.2,"class") > [1] "aovlist" "listof" > > This class isn't suitable for the glht function neither TukeyHSD function. > > I can extract a object of class "aov" through: > > #> ls.str(pat="AnovaModel.2") > AnovaModel.2 : List of 3 > $ (Intercept):List of 9 > $ patient :List of 9 > $ Within :List of 6 > #> AnovaModel.3<-AnovaModel.2$patient > #> attr(AnovaModel.3,"class") > [1] "aov" "lm" > > This model is suitable for confidence intervals: > > #> confint(AnovaModel.3) > 2.5 % 97.5 % > trat2 -0.9107959 0.5462655 > trat3 -0.9848351 0.4722263 > trat4 -1.2997235 0.1573379 > trat5 -1.0623118 0.3947496 > > But it is useless for multiple comparisons: > > #> Pairs <- glht(AnovaModel.3, linfct = mcp(Tratamiento = "Tukey")) > Error in model.matrix.default(model, data = structure(list(method = > lm), .Names = "method"), : > model structure is invalid > Error in factor_contrasts(model) : > no 'model.matrix' method for 'model' found! > > #> TukeyHSD(AnovaModel.3, "Tratamiento") > Error in model.tables.aov(x, "means") : > this fit does not inherit from "lm" > > Another way of analysis is through lme4 package: > > #> AnovaModel.3 <- lmer(VR ~ trat+(1|trat/patient), data=Mesures) > > But the anova test isn't the correct ratio above shown: > > #> anova(AnovaModel.3) > Analysis of Variance Table > Df Sum Sq Mean Sq F value > trat 4 0.43112 0.10778 0.2094 > > Other options for lmer: > > #> anova(AnovaModel.4) > Analysis of Variance Table > Df Sum Sq Mean Sq F value > trat 4 1.5801 0.39501 0.7675 > > Neither is the F ratio searched and above shown: > > #> summary(AnovaModel.2) > > Error: patient > Df Sum Sq Mean Sq F value Pr(>F) > trat 4 1.580 0.3950 0.821 0.541 > Residuals 10 4.811 0.4811 > > I can also take the mean of treatment, but if a measures of some > patient fails, and the design becomes unbalanced, the results would be > very difficult. > > Can anyone give me the right model and multiple comparisons? > > Thanks very much for the help and please excuse for my bad english. > > > > > -- > _____---^---_____ > > Univ. de Extremadura > Dept. Matemáticas. > Despacho B29 > Tf: + 34 924 289 300 > Ext. 86823 >-- _____---^---_____ Univ. de Extremadura Dept. Matemáticas. Despacho B29 Tf: + 34 924 289 300 Ext. 86823 [[alternative HTML version deleted]]