Experimentos con factores aleatorios en R.
Orlando Moscote Flórez
2024-09-26
Más aplicaciones estadísticas en:
https://rpubs.com/orlandoan
Los concetpos sobre este tipo de experimentos pueden encontarse en el libro Diseño y análisis de experimentos de Douglas Montgomery.
Los ejemplos utilizados son tomados del libro citado arriba.
1.Modelo con efectos aleatorios.
1.1. Datos.
Página 514 del libro de Montgomery.
resis<-c(98,97,99,96,91,90,93,92,96,95,97,95,95,96,99,98)
telar<-factor(rep(1:4,each=4))
datos<-data.frame(telar,resis)
datos## telar resis
## 1 1 98
## 2 1 97
## 3 1 99
## 4 1 96
## 5 2 91
## 6 2 90
## 7 2 93
## 8 2 92
## 9 3 96
## 10 3 95
## 11 3 97
## 12 3 95
## 13 4 95
## 14 4 96
## 15 4 99
## 16 4 98boxplot(resis ~ telar)
1.2. Análisis de varianza.
ana<-aov(resis~telar)
summary(ana)## Df Sum Sq Mean Sq F value Pr(>F)
## telar 3 89.19 29.729 15.68 0.000188 ***
## Residuals 12 22.75 1.896
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11.3.Componentes de la varianza.
library(lme4)## Cargando paquete requerido: Matrixalea<-lmer(resis ~(1|telar), data =datos)
summary(alea)## Linear mixed model fit by REML ['lmerMod']
## Formula: resis ~ (1 | telar)
## Data: datos
##
## REML criterion at convergence: 63.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.38018 -0.57260 -0.04342 0.82574 1.52491
##
## Random effects:
## Groups Name Variance Std.Dev.
## telar (Intercept) 6.958 2.638
## Residual 1.896 1.377
## Number of obs: 16, groups: telar, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 95.438 1.363 70.01plot(alea)
2. Diseño factorial de dos factores aleatorios.
2.1. Datos.
Página 519 del libro de Montgomery.
medida<-c(21,24,20,27,19,23,22,19,24,25,21,18,23,24,29,26,20,19,25,19,20,23,21,
27,18,21,21,17,23,23,20,19,25,24,30,26,20,21,26,19,20,24,19,28,19,24,22,18,25,
26,20,17,25,23,30,25,19,19,25,18,20,24,21,26,18,21,24,20,23,25,20,19,25,25,28,
26,20,19,24,17,19,23,20,27,18,23,22,19,24,24,21,18,25,24,31,25,20,21,25,19,21,
24,22,28,21,22,20,18,24,25,20,19,25,25,30,27,20,23,25,17)
pieza<-factor(c(rep(1:20,6)))
operador<-factor(rep(c(1,2,3),c(40,40,40)))
datos<-data.frame(pieza,operador,medida)
datos## pieza operador medida
## 1 1 1 21
## 2 2 1 24
## 3 3 1 20
## 4 4 1 27
## 5 5 1 19
## 6 6 1 23
## 7 7 1 22
## 8 8 1 19
## 9 9 1 24
## 10 10 1 25
## 11 11 1 21
## 12 12 1 18
## 13 13 1 23
## 14 14 1 24
## 15 15 1 29
## 16 16 1 26
## 17 17 1 20
## 18 18 1 19
## 19 19 1 25
## 20 20 1 19
## 21 1 1 20
## 22 2 1 23
## 23 3 1 21
## 24 4 1 27
## 25 5 1 18
## 26 6 1 21
## 27 7 1 21
## 28 8 1 17
## 29 9 1 23
## 30 10 1 23
## 31 11 1 20
## 32 12 1 19
## 33 13 1 25
## 34 14 1 24
## 35 15 1 30
## 36 16 1 26
## 37 17 1 20
## 38 18 1 21
## 39 19 1 26
## 40 20 1 19
## 41 1 2 20
## 42 2 2 24
## 43 3 2 19
## 44 4 2 28
## 45 5 2 19
## 46 6 2 24
## 47 7 2 22
## 48 8 2 18
## 49 9 2 25
## 50 10 2 26
## 51 11 2 20
## 52 12 2 17
## 53 13 2 25
## 54 14 2 23
## 55 15 2 30
## 56 16 2 25
## 57 17 2 19
## 58 18 2 19
## 59 19 2 25
## 60 20 2 18
## 61 1 2 20
## 62 2 2 24
## 63 3 2 21
## 64 4 2 26
## 65 5 2 18
## 66 6 2 21
## 67 7 2 24
## 68 8 2 20
## 69 9 2 23
## 70 10 2 25
## 71 11 2 20
## 72 12 2 19
## 73 13 2 25
## 74 14 2 25
## 75 15 2 28
## 76 16 2 26
## 77 17 2 20
## 78 18 2 19
## 79 19 2 24
## 80 20 2 17
## 81 1 3 19
## 82 2 3 23
## 83 3 3 20
## 84 4 3 27
## 85 5 3 18
## 86 6 3 23
## 87 7 3 22
## 88 8 3 19
## 89 9 3 24
## 90 10 3 24
## 91 11 3 21
## 92 12 3 18
## 93 13 3 25
## 94 14 3 24
## 95 15 3 31
## 96 16 3 25
## 97 17 3 20
## 98 18 3 21
## 99 19 3 25
## 100 20 3 19
## 101 1 3 21
## 102 2 3 24
## 103 3 3 22
## 104 4 3 28
## 105 5 3 21
## 106 6 3 22
## 107 7 3 20
## 108 8 3 18
## 109 9 3 24
## 110 10 3 25
## 111 11 3 20
## 112 12 3 19
## 113 13 3 25
## 114 14 3 25
## 115 15 3 30
## 116 16 3 27
## 117 17 3 20
## 118 18 3 23
## 119 19 3 25
## 120 20 3 172.2. Análisis de varianza.
ana2<-aov(medida~pieza+operador+pieza:operador)
summary(ana2)## Df Sum Sq Mean Sq F value Pr(>F)
## pieza 19 1185.4 62.39 62.915 <2e-16 ***
## operador 2 2.6 1.31 1.319 0.275
## pieza:operador 38 27.0 0.71 0.718 0.861
## Residuals 60 59.5 0.99
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.3. Componentes de la varianza.
alea2<-lmer(medida ~(1|pieza)+(1|operador)+(1|pieza:operador), data =datos)## boundary (singular) fit: see help('isSingular')summary(alea2)## Linear mixed model fit by REML ['lmerMod']
## Formula: medida ~ (1 | pieza) + (1 | operador) + (1 | pieza:operador)
## Data: datos
##
## REML criterion at convergence: 409.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0313 -0.6595 0.1270 0.5374 2.7345
##
## Random effects:
## Groups Name Variance Std.Dev.
## pieza:operador (Intercept) 0.00000 0.0000
## pieza (Intercept) 10.25130 3.2018
## operador (Intercept) 0.01063 0.1031
## Residual 0.88316 0.9398
## Number of obs: 120, groups: pieza:operador, 60; pieza, 20; operador, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 22.3917 0.7235 30.95
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')3. Modelo reducido.
En el modelo anterior, la interacción resultó no ser significativa. Entonces se corre el modelo sin la interacción.
3.1. Análisis de varianza.
ana3<-aov(medida~pieza+operador)
summary(ana3)## Df Sum Sq Mean Sq F value Pr(>F)
## pieza 19 1185.4 62.39 70.645 <2e-16 ***
## operador 2 2.6 1.31 1.481 0.232
## Residuals 98 86.5 0.88
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.2. componentes de la varianza.
alea3<-lmer(medida ~(1|pieza)+(1|operador), data =datos)
summary(alea3)## Linear mixed model fit by REML ['lmerMod']
## Formula: medida ~ (1 | pieza) + (1 | operador)
## Data: datos
##
## REML criterion at convergence: 409.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0313 -0.6595 0.1270 0.5374 2.7345
##
## Random effects:
## Groups Name Variance Std.Dev.
## pieza (Intercept) 10.25127 3.2018
## operador (Intercept) 0.01063 0.1031
## Residual 0.88316 0.9398
## Number of obs: 120, groups: pieza, 20; operador, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 22.3917 0.7235 30.954. Modelo mixto con 2 factores.
4.1. Datos.
Se utilzan los datos del ejemplo anterior pero suponiendo que los operadores son el factor fijo.
4.2. Análisis de varianza.
ana4<-aov(medida~pieza+operador+pieza:operador)
summary(ana4)## Df Sum Sq Mean Sq F value Pr(>F)
## pieza 19 1185.4 62.39 62.915 <2e-16 ***
## operador 2 2.6 1.31 1.319 0.275
## pieza:operador 38 27.0 0.71 0.718 0.861
## Residuals 60 59.5 0.99
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 14.3. Componentes de la varianza.
alea4<-lmer(medida ~(1|pieza)+operador+(1|pieza:operador), data =datos)## boundary (singular) fit: see help('isSingular')summary(alea4)## Linear mixed model fit by REML ['lmerMod']
## Formula: medida ~ (1 | pieza) + operador + (1 | pieza:operador)
## Data: datos
##
## REML criterion at convergence: 409.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1809 -0.7183 0.1445 0.5253 2.5849
##
## Random effects:
## Groups Name Variance Std.Dev.
## pieza:operador (Intercept) 0.0000 0.0000
## pieza (Intercept) 10.2513 3.2018
## Residual 0.8832 0.9398
## Number of obs: 120, groups: pieza:operador, 60; pieza, 20
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 22.3000 0.7312 30.498
## operador2 -0.0250 0.2101 -0.119
## operador3 0.3000 0.2101 1.428
##
## Correlation of Fixed Effects:
## (Intr) oprdr2
## operador2 -0.144
## operador3 -0.144 0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')—|—|
O.M.F.
—|—|