Unidad 6 experim

Paula Esquivias Ruiz-Dana

17 de April, 2024

Contents
0.1 Ej Problema 8.5 libro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
0.1.1 Ajustemos el modelo de dos factores con interacción . . . . . . . . . . . . . . . . . . . . . . . 6
0.1.2 Comprobación de las suposiciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
0.2 Actividad 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
0.3 gráficas interacciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
0.4 Estimación de las medias y los efectos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
#teoria

conc <- c ( 16.5 , 14.5 , 39.1 , 32.0 , 18.4 , 11.0 , 26.2 , 23.8 , 12.7 , 10.8 , 21.3 , 28.8 , 14.0 , 14.3
hormone <- rep ( c ( "Treatment" , "Treatment" , "No Treatment" , "No Treatment" ), 5 )
gender <- rep ( c ( "Female" , "Male" ), 10 )
df <- data.frame (conc,hormone,gender)
interaction.plot (hormone,gender,conc)

1
 gender
30

Female
Male
mean of conc

25
20
15

No Treatment Treatment

hormone

interaction.plot (gender, hormone,conc)

2
 hormone
30

No Treatment
Treatment
mean of conc

25
20
15

Female Male

gender
No
hay interacción

y <- c ( 4 , 5 , 6 , 5 , 7 , 9 , 8 , 12 , 10 , 12 , 11 , 9 , 6 , 6 , 4 , 4 , 13 , 15 , 12 , 12 , 12 , 13 ,
temper <- rep ( rep ( c ( "frío" , "templada" , "caliente" ), each= 4 ), 2 )
deter <- rep ( c ( "A" , "B" ), each= 12 )
df <- data.frame (y,temper,deter)
interaction.plot (deter, temper,y)

3
 temper
12

templada
caliente
frío
10
mean of y

8
6

A B

deter

interaction.plot (temper, deter,y)

4
 deter
12

B
A
10
mean of y

8
6

caliente frío templada

temper
Sí hay
interacción.

0.1 Ej Problema 8.5 libro

y <- c ( 9 , 43 , 60 , 77 , 13 , 48 , 65 , 70 , 12 , 57 , 70 , 91 , 15 , 66 , 75 , 97 , 13 , 58 , 78 , 108
shape <- as.factor ( rep ( 1 : 4 , 6 ))
treat <- as.factor ( rep ( 1 : 2 , each = 12 ))
df <- data.frame (y, shape, treat)
str (df)

## ’data.frame’: 24 obs. of 3 variables:

## $ y : num 9 43 60 77 13 48 65 70 12 57 ...
## $ shape: Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
## $ treat: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...

head(df)

## y shape treat
## 1 9 1 1
## 2 43 2 1
## 3 60 3 1
## 4 77 4 1
## 5 13 1 1
## 6 48 2 1

5
0.1.1 Ajustemos el modelo de dos factores con interacción

modelo <- aov (y ~ shape + treat + shape : treat, data = df)

Fíjate con el uso del operador : para designar el término de la interacción dentro de la fórmula. También es posible
usar el operador * lo cual incluye tanto el término de interacción como los términos principales. Así la anterior
instrucción es equivalente a:

modelo <- aov (y ~ shape * treat, data = df)

0.1.2 Comprobación de las suposiciones

Para la comprobación de la homocedasticidad, podemos hacer un grafico de residuos frente a valores ajustados, y
opcionalmente, un test de Levene para testar si la varianza en las 8 condiciones experimentales es la misma. Para
aplicar el test de Levene (o alternativamente, el de Bartlett) debemos crear una variable nueva mediante la función
interaction con la que crear una variable factor con los cruces de los niveles de los dos factores.

#normalidad
plot(modelo)

Residuals vs Fitted

12
10

23
5
Residuals

0
−5
−10

20 40 60 80 100

Fitted values
aov(y ~ shape * treat)

6
 Q−Q Residuals

12
2
Standardized residuals

23
1
0
−1

−2 −1 0 1 2

Theoretical Quantiles
aov(y ~ shape * treat)

Scale−Location
1.5

12

823
Standardized residuals

1.0
0.5
0.0

20 40 60 80 100

Fitted values
aov(y ~ shape * treat)

7
 Constant Leverage:
Residuals vs Factor Levels

12
2

23
Standardized residuals

1
0
−1

8
−2

shape :
1 2 3 4

Factor Level Combinations

library(car)

## Loading required package: carData

condition <- with (df, interaction (shape, treat))

leveneTest (y ~ condition, data = df)

## Levene’s Test for Homogeneity of Variance (center = median)

## Df F value Pr(>F)
## group 7 0.45 0.856
## 16

Levene no resulta significativo, aceptaremos mantener la homocedasticidad. Podemos comprobar la normalidad,

mediante un qqplot y el test de Shapiro.

shapiro.test ( residuals (modelo))

##
## Shapiro-Wilk normality test
##
## data: residuals(modelo)
## W = 0.96243, p-value = 0.4892

Aceptamos la hipótesis de normalidad

8
anova(modelo)

## Analysis of Variance Table

##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## shape 3 19407.5 6469.2 143.4932 8.934e-12 ***
## treat 1 1305.4 1305.4 28.9547 6.122e-05 ***
## shape:treat 3 237.5 79.2 1.7557 0.1961
## Residuals 16 721.3 45.1
## ---
## Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1

Hay efecto de la forma y del tratamiento pero no hay interacción.

interaction.plot (shape,treat,y)
100

treat
2
80

1
mean of y

60
40
20

1 2 3 4

shape

interaction.plot (treat,shape,y)

9
 100

shape

4
80

3
2
mean of y

1
60
40
20

1 2

treat

Las medias por niveles y por conciciones experimentales son:

model.tables (modelo, type = "mean" )

## Tables of means
## Grand mean
##
## 58.625
##
## shape
## shape
## 1 2 3 4
## 13.67 57.50 73.00 90.33
##
## treat
## treat
## 1 2
## 51.25 66.00
##
## shape:treat
## treat
## shape 1 2
## 1 11.33 16.00
## 2 49.33 65.67
## 3 65.00 81.00
## 4 79.33 101.33

Y para obtener las estimación de los efectos principales e interacciones, haremos

10
model.tables (modelo, type = "effects" )

## Tables of effects
##
## shape
## shape
## 1 2 3 4
## -44.96 -1.13 14.37 31.71
##
## treat
## treat
## 1 2
## -7.375 7.375
##
## shape:treat
## treat
## shape 1 2
## 1 5.042 -5.042
## 2 -0.792 0.792
## 3 -0.625 0.625
## 4 -3.625 3.625

Dado que los factores principales han resultado significativos, deberemos realizar comparaciones por parejas para
determinar que condición experimental es óptima. Por ejemplo, mediante el test de Tukey

library (agricolae)
HSD.test (modelo, "shape" , console = TRUE )

##
## Study: modelo ~ "shape"
##
## HSD Test for y
##
## Mean Square Error: 45.08333
##
## shape, means
##
## y std r se Min Max Q25 Q50 Q75
## 1 13.66667 3.669696 6 2.741147 9 20 12.25 13.0 14.50
## 2 57.50000 11.077003 6 2.741147 43 73 50.25 57.5 64.00
## 3 73.00000 10.583005 6 2.741147 60 90 66.25 72.5 77.25
## 4 90.33333 14.306176 6 2.741147 70 108 80.50 94.0 98.50
##
## Alpha: 0.05 ; DF Error: 16
## Critical Value of Studentized Range: 4.046093
##
## Minimun Significant Difference: 11.09094
##
## Treatments with the same letter are not significantly different.
##
## y groups
## 4 90.33333 a
## 3 73.00000 b
## 2 57.50000 c
## 1 13.66667 d

11
HSD.test (modelo, "treat" , console = TRUE )

##
## Study: modelo ~ "treat"
##
## HSD Test for y
##
## Mean Square Error: 45.08333
##
## treat, means
##
## y std r se Min Max Q25 Q50 Q75
## 1 51.25 27.15653 12 1.938284 9 91 35.5 58.5 70.00
## 2 66.00 33.37664 12 1.938284 13 108 48.5 74.0 91.75
##
## Alpha: 0.05 ; DF Error: 16
## Critical Value of Studentized Range: 2.997999
##
## Minimun Significant Difference: 5.810973
##
## Treatments with the same letter are not significantly different.
##
## y groups
## 2 66.00 a
## 1 51.25 b

En este ejemplo, dado que el término de interacción no es significativo, para determinar la condición experimental
óptima (mayor producción de resina) es suficiente examinar en qué combinación de nivel del factor A y de nivel
de factor B tenemos máxima respuesta. A la vista de las comparaciones múltiples, concluimos que la condición
óptima corresponde a la forma rectangular y al tratamiento sin ácido. Siendo la respuesta media en dicho condición
experimental de 101.33 gramos de resina por agujero. Notar que si la interacción resultara significativa, deberíamos
hacer comparaciones múltiples de la interacción para determinar la condición experimental óptima.
Si la interacción hubieses resultado significativa no podemos interpretar los efectos principales del tratamiento y de
la forma por separado. En éste hipotético caso, para determinar qué combinación o combinaciones son las óptimas
deberíamos trabajar con el factor formado por la combinación de niveles de los dos factores. Las instrucciones
serían.

library (agricolae)
modelo1F <- aov (y ~ condition)
HSD.test (modelo1F, "condition" , console = TRUE )

##
## Study: modelo1F ~ "condition"
##
## HSD Test for y
##
## Mean Square Error: 45.08333
##
## condition, means
##
## y std r se Min Max Q25 Q50 Q75
## 1.1 11.33333 2.081666 3 3.876568 9 13 10.5 12 12.5
## 1.2 16.00000 3.605551 3 3.876568 13 20 14.0 15 17.5
## 2.1 49.33333 7.094599 3 3.876568 43 57 45.5 48 52.5

12
## 2.2 65.66667 7.505553 3 3.876568 58 73 62.0 66 69.5
## 3.1 65.00000 5.000000 3 3.876568 60 70 62.5 65 67.5
## 3.2 81.00000 7.937254 3 3.876568 75 90 76.5 78 84.0
## 4.1 79.33333 10.692677 3 3.876568 70 91 73.5 77 84.0
## 4.2 101.33333 5.859465 3 3.876568 97 108 98.0 99 103.5
##
## Alpha: 0.05 ; DF Error: 16
## Critical Value of Studentized Range: 4.89622
##
## Minimun Significant Difference: 18.98053
##
## Treatments with the same letter are not significantly different.
##
## y groups
## 4.2 101.33333 a
## 3.2 81.00000 b
## 4.1 79.33333 b
## 2.2 65.66667 bc
## 3.1 65.00000 bc
## 2.1 49.33333 c
## 1.2 16.00000 d
## 1.1 11.33333 d

0.2 Actividad 2

rate_delamination <- c (.83 , .68 , .18 , .25 , .78 , .90 , .16 ,.20 , .86 , .72 , .30 , .10 , .67 , .81 ,
firing_profile_time <- as.factor ( rep ( c ( "8 hrs" , "8 hrs" , "13 hrs" , "13 hrs" ), 4 ))
furnace_airflow <- as.factor ( rep ( c ( "Low" , "High" , "Low" , "High" ), 4 ))
laser <- as.factor ( c ( rep ( "Old" , 8 ), rep ( "New" , 8 )))
datos <- data.frame (rate_delamination,firing_profile_time,furnace_airflow,laser)

head(datos)

## rate_delamination firing_profile_time furnace_airflow laser

## 1 0.83 8 hrs Low Old
## 2 0.68 8 hrs High Old
## 3 0.18 13 hrs Low Old
## 4 0.25 13 hrs High Old
## 5 0.78 8 hrs Low Old
## 6 0.90 8 hrs High Old

str(datos)

## ’data.frame’: 16 obs. of 4 variables:

## $ rate_delamination : num 0.83 0.68 0.18 0.25 0.78 0.9 0.16 0.2 0.86 0.72 ...
## $ firing_profile_time: Factor w/ 2 levels "13 hrs","8 hrs": 2 2 1 1 2 2 1 1 2 2 ...
## $ furnace_airflow : Factor w/ 2 levels "High","Low": 2 1 2 1 2 1 2 1 2 1 ...
## $ laser : Factor w/ 2 levels "New","Old": 2 2 2 2 2 2 2 2 1 1 ...

library (gplots)

##
## Attaching package: ’gplots’

13
## The following object is masked from ’package:stats’:
##
## lowess

par(mfrow=c(1,3))
plotmeans (datos$rate_delamination ~ datos$firing_profile_time)
plotmeans (datos$rate_delamination ~ datos$furnace_airflow)
plotmeans (datos$rate_delamination ~ datos$laser)
0.8

0.7

0.7
0.7

0.6

0.6
datos$rate_delamination

datos$rate_delamination

datos$rate_delamination
0.6

0.5

0.5
0.5
0.4

0.4

0.4
0.3

0.3

0.3
0.2

0.2

0.2

n=8 n=8 n=8 n=8 n=8 n=8

13 hrs 8 hrs High Low New Old

datos$firing_profile_time datos$furnace_airflow datos$laser

Parece que firing_profile tendrá efecto

0.3 gráficas interacciones

with (datos, interaction.plot (firing_profile_time,furnace_airflow,datos$rate_delamination))

14
 0.8
mean of datos$rate_delamination

furnace_airflow
0.7

Low
High
0.6
0.5
0.4
0.3
0.2

13 hrs 8 hrs

firing_profile_time

with (datos, interaction.plot (firing_profile_time,laser,datos$rate_delamination))

15
 0.8
mean of datos$rate_delamination

laser
0.7

Old
New
0.6
0.5
0.4
0.3
0.2

13 hrs 8 hrs

firing_profile_time

with (datos, interaction.plot (furnace_airflow,laser,datos$rate_delamination))

16
 mean of datos$rate_delamination

laser
0.50

New
Old
0.48
0.46
0.44

High Low

furnace_airflow
##
Diseño experimental
El diseño experimental de tres factores A, B y C, donde el factor A tiene a niveles, el factor B tiene b niveles, el
factor C tiene c niveles, y hay n unidades experimentales asignadas a cada combinación factor-nivel. El modelo es:
yij = µ + αi + βj + γk + αβij + αγik + βγjk + αβγijk + ϵijkl
donde i=1,. . . , a, j=1,. . . , b, k=1, . . . , b, l=1,. . . ,n yijk denota la respuesta l en el i-ésimo nivel de A y j-ésimo
nivel de B y k-ésimonivel de C, µ la media general, αi efecto del nivel i del factor principal A, βj efecto del nivel j
del factor principal B, γk efecto del nivel k del factor principal C, αβij la interacción del nivel i del factor A con el
nivel j del factor B, αγik la interacción del nivel i del factor A con el nivel k del factor C , βγjk la interacción del
nivel j del factor B con el nivel k del factor C , αβγijk la interacción del nivel i del factor A con el nivel j del factor
B con el nivel k del factor C , ϵijkl el error es independiente y normalmente distribuido con media cero y varianza
α 2 . Hay un total de N = n × a × b × c unidades experimentales (observaciones).

model.3F <- aov (rate_delamination ~ firing_profile_time * furnace_airflow * laser, data= df)

summary(model.3F)

## Df Sum Sq Mean Sq F value Pr(>F)

## firing_profile_time 1 1.3748 1.3748 210.489 4.99e-07
## furnace_airflow 1 0.0028 0.0028 0.422 0.534
## laser 1 0.0014 0.0014 0.215 0.655
## firing_profile_time:furnace_airflow 1 0.0014 0.0014 0.215 0.655
## firing_profile_time:laser 1 0.0008 0.0008 0.116 0.742
## furnace_airflow:laser 1 0.0086 0.0086 1.310 0.285
## firing_profile_time:furnace_airflow:laser 1 0.0116 0.0116 1.769 0.220
## Residuals 8 0.0523 0.0065
##
## firing_profile_time ***
## furnace_airflow

17
## laser
## firing_profile_time:furnace_airflow
## firing_profile_time:laser
## furnace_airflow:laser
## firing_profile_time:furnace_airflow:laser
## Residuals
## ---
## Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1

Resulta un único factor principal significativo: firing_profile_time . Podemos interpretar fácilmente esta tabla
dado que ninguna interacción es significativa.

0.4 Estimación de las medias y los efectos

model.tables (model.3F, type= "means" )

## Tables of means
## Grand mean
##
## 0.488125
##
## firing_profile_time
## firing_profile_time
## 13 hrs 8 hrs
## 0.1950 0.7813
##
## furnace_airflow
## furnace_airflow
## High Low
## 0.4750 0.5013
##
## laser
## laser
## New Old
## 0.4788 0.4975
##
## firing_profile_time:furnace_airflow
## furnace_airflow
## firing_profile_time High Low
## 13 hrs 0.1725 0.2175
## 8 hrs 0.7775 0.7850
##
## firing_profile_time:laser
## laser
## firing_profile_time New Old
## 13 hrs 0.1925 0.1975
## 8 hrs 0.7650 0.7975
##
## furnace_airflow:laser
## laser
## furnace_airflow New Old
## High 0.4425 0.5075
## Low 0.5150 0.4875
##

18
## firing_profile_time:furnace_airflow:laser
## , , laser = New
##
## furnace_airflow
## firing_profile_time High Low
## 13 hrs 0.120 0.265
## 8 hrs 0.765 0.765
##
## , , laser = Old
##
## furnace_airflow
## firing_profile_time High Low
## 13 hrs 0.225 0.170
## 8 hrs 0.790 0.805

model.tables (model.3F, type= "effects" )

## Tables of effects
##
## firing_profile_time
## firing_profile_time
## 13 hrs 8 hrs
## -0.29313 0.29313
##
## furnace_airflow
## furnace_airflow
## High Low
## -0.013125 0.013125
##
## laser
## laser
## New Old
## -0.009375 0.009375
##
## firing_profile_time:furnace_airflow
## furnace_airflow
## firing_profile_time High Low
## 13 hrs -0.009375 0.009375
## 8 hrs 0.009375 -0.009375
##
## firing_profile_time:laser
## laser
## firing_profile_time New Old
## 13 hrs 0.006875 -0.006875
## 8 hrs -0.006875 0.006875
##
## furnace_airflow:laser
## laser
## furnace_airflow New Old
## High -0.023125 0.023125
## Low 0.023125 -0.023125
##
## firing_profile_time:furnace_airflow:laser
## , , laser = New
##
## furnace_airflow

19
## firing_profile_time High Low
## 13 hrs -0.026875 0.026875
## 8 hrs 0.026875 -0.026875
##
## , , laser = Old
##
## furnace_airflow
## firing_profile_time High Low
## 13 hrs 0.026875 -0.026875
## 8 hrs -0.026875 0.026875

20