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]]
______________________________________________
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.
