Tabla de análisis de varianza para ensayos cruzados 2x2

Question

Estoy tratando de reproducir la tabla anova para un 2x2 diseño cruzado. He utilizado los datos que se enumeran en las tablas 2.1 y 2.2 del "Diseño y Análisis de cross-over" pruebas " de libro, por Jones y Kenward. Me estoy dando todo el código para que pueda reproducir los cálculos.

g1AB <- read.table(textConnection("
Label Per1 Per2
7 121.905 116.667
8 218.5 200.5
9 235 217.143
13 250 196.429
14 186.19 185.5
15 231.563 221.842
17 443.25 420.5
21 198.421 207.692
22 270.5 213.158
28 360.476 384
35 229.75 188.25
36 159.091 221.905
37 255.882 253.571
38 279.048 267.619
41 160.556 163
44 172.105 182.381
58 267 313
66 230.75 211.111
71 271.19 257.619
76 276.25 222.105
79 398.75 404
80 67.778 70.278
81 195 223.158
82 325 306.667
86 368.077 362.5
89 228.947 227.895
90 236.667 220
")->con,header=T)
close(con)

g2BA <- read.table(textConnection("
Label Per1 Per2
3 138.333 138.571
10 225 256.25
11 392.857 381.429
16 190 233.333
18 191.429 228
23 226.19 267.143
24 201.905 193.5
26 134.286 128.947
27 238 248.5
29 159.5 140
30 232.75 276.563
32 172.308 170
33 266 305
39 171.333 186.333
43 194.737 191.429
47 200 222.619
51 146.667 183.81
52 208 241.667
55 208.75 218.81
59 271.429 225
68 143.81 188.5
70 104.444 135.238
74 145.238 152.857
77 215.385 240.476
78 306 288.333
83 160.526 150.476
84 353.81 369.048
85 293.889 308.095
99 371.190 404.762
")->con,header=T)
close(con)

n1 <- nrow(g1AB)
n2 <- nrow(g2BA)
p <- 2 # periods

Algunas cantidades que pueden ser útiles:

y11. <- sum(g1AB$Per1)
y12. <- sum(g1AB$Per2)
y21. <- sum(g2BA$Per1)
y22. <- sum(g2BA$Per2)
y1.. <- sum(c(g1AB$Per1,g1AB$Per2))
y2.. <- sum(c(g2BA$Per1,g2BA$Per2))
y... <- y1.. + y2..

y11.bar <- 1/n1 * y11.
y12.bar <- 1/n1 * y12.
y21.bar <- 1/n2 * y21.
y22.bar <- 1/n2 * y22.

y1..bar <- 1/(p*n1)*(y11.+y12.)
y2..bar <- 1/(p*n2)*(y21.+y22.)

y...bar <- 1/(p*(n1+n2)) * (y1.. + y2..)

Con el fin de realizar el análisis de la varianza de R, he creado el siguiente conjunto de datos

mydata1 <- data.frame(PEFR=g1AB$Per1,Subjects=g1AB$Label,Time=1,Groups="AB",Treatment=1)
mydata2 <- data.frame(PEFR=g2BA$Per1,Subjects=g2BA$Label,Time=1,Groups="BA",Treatment=2)
mydata3 <- data.frame(PEFR=g1AB$Per2,Subjects=g1AB$Label,Time=2,Groups="AB",Treatment=2)
mydata4 <- data.frame(PEFR=g2BA$Per2,Subjects=g2BA$Label,Time=2,Groups="BA",Treatment=1)
mydata <- rbind(mydata1,mydata2,mydata3,mydata4)
mydata$Subjects<-factor(mydata$Subjects)

(He leído que a) La correcta tabla anova se obtiene mediante el siguiente comando

res <- summary(PEFR.aov <- aov(PEFR~Groups+Time+Treatment+Error(Subjects), data=mydata))

Error: Subjects
          Df Sum Sq Mean Sq F value Pr(>F)
Groups     1  10573   10573   0.899  0.347
Residuals 54 634866   11757               

Error: Within
          Df Sum Sq Mean Sq F value  Pr(>F)   
Time       1    480   479.6   1.470 0.23061   
Treatment  1   3026  3026.1   9.276 0.00359 **
Residuals 54  17617   326.2 

1) ¿Puede por favor que me explique el uso del término de Error en la fórmula?

2) La tabla anova es la misma que la de la tabla 2.10 del libro, a excepción de la SS para la Época) que tiene un valor de 396.858 en lugar de 480. De acuerdo a la Tabla 2.8 esta cantidad se calcula mediante

> n1*n2/(2*(n1+n2))*(y11.bar-y12.bar+y21.bar-y22.bar)^2
[1] 396.8583

También, el Tratamiento de las SS se calcula por

n1*n2/(2*(n1+n2))*(y11.bar-y12.bar-y21.bar+y22.bar)^2
[1] 3026.12

Si la fórmula fue especificado como este

PEFR~Groups+Treatment+Time+Error(Subjects)

la SS por el Tiempo sería correcto, pero el Tratamiento de las SS tendría un valor de 3109. ¿Cómo puedo obtener 3026 para el Tratamiento y 397 de Tiempo, como yo lo veo en el libro? Lo que debe tener cuidado en el orden de las variables especificadas en la fórmula?

Gracias por su tiempo.

Answer 1

1) El término de Error, básicamente, dice el aov() función que, además de la aleatoriedad entre cada observación, algunas de las observaciones están relacionadas entre ellas por ser las observaciones sobre el mismo tema. Así que una componente aleatoria de la PERF valor de cada tema, antes de calcular la componente aleatoria de cada observación individual sobre el tema.

2) el problema con la suma de los cuadrados parece que está relacionado con la cuestión de Tipo III suma de cuadrados (el predeterminado en algunas otras estadísticas de los paquetes). Tipo III SS son los marginales de la suma de cuadrados explicada por una variable o grupo de variables después de todo la varianza explicada por las otras variables tomadas en cuenta. Como usted ha descubierto, R sólo informa de este para el último la variable introducida en el modelo. Aunque a veces controversial, hay buenas razones para esto y yo mismo prefiero el R - defecto que me gusta a mi las sumas de cuadrados en una tabla ANOVA para añadir el número de la derecha. Como usted ha descubierto, hay un montón de maneras para generar la marginal de la suma de los cuadrados de un grupo de variables si desea probar.

Para responder a su pregunta ¿qué se debe tener cuidado en el orden... usted debe poner las variables que están interesados en probar la última. De esa manera usted está implícitamente prueba dos modelos anidados, uno con todas las variables, incluyendo los que están en cuestión; otro sólo con las otras variables (por ejemplo molestia variables de confusión, o que son sólo positivos están en la mejor modelo y, por tanto, no necesita ser probado). Usted puede entonces llevar a cabo una prueba de hipótesis sobre la suma de los cuadrados de los reportados para la última variable/s.