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

