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