Diseño completamente aleatorio
Orlando Moscote Flórez
7/3/2022
1. Datos
Se desarrola el ejemplo del libro Diseño y anális de experimentos de Douglas Montgomery. Página 60. Se mide la resistencia a la tensión (resis) de una fibra sintética y el peso porcentual de algodón utilizado (peso).
Se trata de un experimento de un solo factor de efectos fijos y diseño balanceado.
Los datos se muestran a continuación.
peso<-factor(c(rep(15,5),rep(20,5),rep(25,5),rep(30,5),rep(35,5)))
resis<-c(7,7,15,11,9,12,17,12,18,18,14,18,18,19,19,19,25,22,19,23,7,10,11,15,11)
datos<-data.frame(peso,resis)2. Análisis de varianza.
anava<-aov(resis~peso)
summary(anava)## Df Sum Sq Mean Sq F value Pr(>F)
## peso 4 475.8 118.94 14.76 9.13e-06 ***
## Residuals 20 161.2 8.06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13. Estimación de parámetros.
3.1 media global.
media<-mean(resis);media## [1] 15.043.2 media de los tratamientos.
trata<- lsmeans(anava, "peso");trata## peso lsmean SE df lower.CL upper.CL
## 15 9.8 1.27 20 7.15 12.4
## 20 15.4 1.27 20 12.75 18.0
## 25 17.6 1.27 20 14.95 20.2
## 30 21.6 1.27 20 18.95 24.2
## 35 10.8 1.27 20 8.15 13.4
##
## Confidence level used: 0.953.3 Efectos.
ma<-tapply(resis, peso, mean)-mean(resis);ma##efectos## 15 20 25 30 35
## -5.24 0.36 2.56 6.56 -4.24regre<-lm(resis~peso)
eme <- emmeans(regre, "peso")
eff_size(eme, sigma = sigma(regre), edf = df.residual(regre))## contrast effect.size SE df lower.CL upper.CL
## 15 - 20 -1.973 0.705 20 -3.443 -0.5015
## 15 - 25 -2.747 0.767 20 -4.348 -1.1469
## 15 - 30 -4.156 0.912 20 -6.059 -2.2538
## 15 - 35 -0.352 0.635 20 -1.677 0.9721
## 20 - 25 -0.775 0.644 20 -2.119 0.5689
## 20 - 30 -2.184 0.721 20 -3.687 -0.6808
## 20 - 35 1.620 0.682 20 0.197 3.0437
## 25 - 30 -1.409 0.671 20 -2.808 -0.0102
## 25 - 35 2.395 0.737 20 0.857 3.9329
## 30 - 35 3.804 0.873 20 1.984 5.6248
##
## sigma used for effect sizes: 2.839
## Confidence level used: 0.953.4 Intervalo de confianza tratamiento 4.
n<-5
me<-tapply(resis,peso,mean);me## 15 20 25 30 35
## 9.8 15.4 17.6 21.6 10.8df<-df.residual(anava)
cme<-deviance(anava)/df
LI<-me[4]+qt(0.025,df)*sqrt(cme/n);LI## 30
## 18.95157LS<-me[4]+qt(0.975,df)*sqrt(cme/n);LS## 30
## 24.248434. Verificación de la adecuación del modelo.
4.1. Residuales sin estructura.
residuos<-resid(anava)
m<-mean(residuos)
d<-sd(residuos)
plot(residuos,pch=18,xlab="peso")
4.2. Histograma.
hist(residuos,prob=TRUE)
curve(dnorm(x,m,d),col="red",add=TRUE)
4.3. Gráfico cuantil-cuantil.
qqnorm(residuos)
qqline(residuos)
4.4. Residuos estandarizados.
dij<-residuos/sqrt(cme);dij## 1 2 3 4 5 6 7
## -0.9862579 -0.9862579 1.8316219 0.4226820 -0.2817880 -1.1975989 0.5635760
## 8 9 10 11 12 13 14
## -1.1975989 0.9158109 0.9158109 -1.2680459 0.1408940 0.1408940 0.4931290
## 15 16 17 18 19 20 21
## 0.4931290 -0.9158109 1.1975989 0.1408940 -0.9158109 0.4931290 -1.3384929
## 22 23 24 25
## -0.2817880 0.0704470 1.4793869 0.0704470mean(dij)## [1] 0var(dij)## [1] 0.83333334.5. residuales en secuecia en el tiempo.
tiempo<-c(15,19,25,12,6,8,14,1,11,3,18,13,20,7,9,22,5,2,24,10,17,21,4,16,23)
plot(tiempo,residuos,pch=19,xlab="tiempo")
abline(h=0)
4.6. Residuales contra los valores ajustados.
pre<-predict(anava);pre## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 9.8 9.8 9.8 9.8 9.8 15.4 15.4 15.4 15.4 15.4 17.6 17.6 17.6 17.6 17.6 21.6
## 17 18 19 20 21 22 23 24 25
## 21.6 21.6 21.6 21.6 10.8 10.8 10.8 10.8 10.8plot(pre,residuos,pch=19)
abline(h=0)
4.7. Normalidad de los residuos.
ad.test(residuos)##
## Anderson-Darling normality test
##
## data: residuos
## A = 0.51857, p-value = 0.1699cvm.test(residuos)##
## Cramer-von Mises normality test
##
## data: residuos
## W = 0.080455, p-value = 0.1949lillie.test(residuos)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuos
## D = 0.16212, p-value = 0.088895. Igualdad de varianzas (Homocedasticidad).
5.1 Prueba de Bartlett.
bartlett.test(resis,peso)##
## Bartlett test of homogeneity of variances
##
## data: resis and peso
## Bartlett's K-squared = 0.93309, df = 4, p-value = 0.91985.2 Prueba de Levene de la mediana.
leveneTest(resis,peso, center=c("median"))## Levene's Test for Homogeneity of Variance (center = c("median"))
## Df F value Pr(>F)
## group 4 0.3179 0.8626
## 205.3. Prueba de Levene de la media.
leveneTest(resis,peso, center=c("mean"))## Levene's Test for Homogeneity of Variance (center = c("mean"))
## Df F value Pr(>F)
## group 4 0.6443 0.6372
## 205.4. Prueba de Hartley.
hartleyTest(resis~peso) ##
## Hartley's maximum F-ratio test of homogeneity of variances
##
## data: resis by peso
## F Max = 2.6047, df = 4, k = 5, p-value = 0.89645.5 Prueba de Cochran.
cochran.test(resis~peso)##
## Cochran test for outlying variance
##
## data: resis ~ peso
## C = 0.27792, df = 5, k = 5, p-value = 1
## alternative hypothesis: Group 15 has outlying variance
## sample estimates:
## 15 20 25 30 35
## 11.2 9.8 4.3 6.8 8.25.6 Prueba de Fligner-Killeen.
fligner.test(resis~peso)##
## Fligner-Killeen test of homogeneity of variances
##
## data: resis by peso
## Fligner-Killeen:med chi-squared = 1.7265, df = 4, p-value = 0.78596. Un modelo de regresión.
x<-c(rep(15,5),rep(20,5),rep(25,5),rep(30,5),rep(35,5))
x1<-x^2
mo1<-lm(resis~x+x1)
summary(mo1)##
## Call:
## lm(formula = resis ~ x + x1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.4686 -1.9714 -0.4686 1.5657 6.9257
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -39.98857 9.78467 -4.087 0.000488 ***
## x 4.59257 0.82771 5.549 1.41e-05 ***
## x1 -0.08857 0.01644 -5.388 2.07e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.439 on 22 degrees of freedom
## Multiple R-squared: 0.5916, Adjusted R-squared: 0.5545
## F-statistic: 15.94 on 2 and 22 DF, p-value: 5.27e-05a1<-mo1$coefficients[1]
a2<-mo1$coefficients[2]
a3<-mo1$coefficients[3]
x2<-x^3
mo2<-lm(resis~x+x1+x2)
summary(mo2)##
## Call:
## lm(formula = resis ~ x + x1 + x2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.4686 -1.4686 -0.4686 2.6457 4.8886
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.611429 39.757436 1.575 0.1302
## x -9.011429 5.196609 -1.734 0.0976 .
## x1 0.481429 0.216046 2.228 0.0369 *
## x2 -0.007600 0.002874 -2.644 0.0152 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.048 on 21 degrees of freedom
## Multiple R-squared: 0.6936, Adjusted R-squared: 0.6499
## F-statistic: 15.85 on 3 and 21 DF, p-value: 1.295e-05b1<-mo2$coefficients[1]
b2<-mo2$coefficients[2]
b3<-mo2$coefficients[3]
b4<-mo2$coefficients[4]
f1<-function(x){
a1+a2*x+a3*x^2}
f2<-function(y){
b1+b2*y+b3*y^2+b4*y^3}
plot(x,resis)
curve(f1,lwd=2,xlim=c(15,40),add=TRUE)
curve(f2,xlim=c(15,40),add=TRUE,lwd=2,col="red")
7. Contrastes.
c1<-c(0,0,0,-1,1)
c2<-c(1,0,1,-1,-1)
c3<-c(1,0,-1,0,0)
c4<-c(-1,4,-1,-1,-1)
matri<-cbind(c1,c2,c3,c4)
contrasts(peso) = matri
mlist<-list("c1"=1,"c2"=2,"c3"=3,"c4"=4)
modelo1<-aov(resis~peso)
summary(modelo1, split=list(peso=mlist))## Df Sum Sq Mean Sq F value Pr(>F)
## peso 4 475.8 118.94 14.757 9.13e-06 ***
## peso: c1 1 291.6 291.60 36.179 7.01e-06 ***
## peso: c2 1 31.3 31.25 3.877 0.062960 .
## peso: c3 1 152.1 152.10 18.871 0.000315 ***
## peso: c4 1 0.8 0.81 0.100 0.754520
## Residuals 20 161.2 8.06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 18. Comparaciones por pares.
8.1 Prueba de Tukey.
tu<-TukeyHSD(anava)
tu## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = resis ~ peso)
##
## $peso
## diff lwr upr p adj
## 20-15 5.6 0.2270417 10.9729583 0.0385024
## 25-15 7.8 2.4270417 13.1729583 0.0025948
## 30-15 11.8 6.4270417 17.1729583 0.0000190
## 35-15 1.0 -4.3729583 6.3729583 0.9797709
## 25-20 2.2 -3.1729583 7.5729583 0.7372438
## 30-20 6.2 0.8270417 11.5729583 0.0188936
## 35-20 -4.6 -9.9729583 0.7729583 0.1162970
## 30-25 4.0 -1.3729583 9.3729583 0.2101089
## 35-25 -6.8 -12.1729583 -1.4270417 0.0090646
## 35-30 -10.8 -16.1729583 -5.4270417 0.0000624plot(tu)
8.2 Prueba de Fisher (Diferencia mínima significativa - LSD).
df<-df.residual(anava) ##agricolae
cme<-deviance(anava)/df
fisher <- LSD.test(resis,peso,df,cme, p.adj="bonferroni", group=TRUE,
main="LSD",alpha=0.05)
fisher## $statistics
## MSerror Df Mean CV t.value MSD
## 8.06 20 15.04 18.87642 3.153401 5.662089
##
## $parameters
## test p.ajusted name.t ntr alpha
## Fisher-LSD bonferroni peso 5 0.05
##
## $means
## resis std r LCL UCL Min Max Q25 Q50 Q75
## 15 9.8 3.346640 5 7.151566 12.44843 7 15 7 9 11
## 20 15.4 3.130495 5 12.751566 18.04843 12 18 12 17 18
## 25 17.6 2.073644 5 14.951566 20.24843 14 19 18 18 19
## 30 21.6 2.607681 5 18.951566 24.24843 19 25 19 22 23
## 35 10.8 2.863564 5 8.151566 13.44843 7 15 10 11 11
##
## $comparison
## NULL
##
## $groups
## resis groups
## 30 21.6 a
## 25 17.6 ab
## 20 15.4 bc
## 35 10.8 c
## 15 9.8 c
##
## attr(,"class")
## [1] "group"plot(fisher)
8.3. Prueba de Duncan
dunca <- duncan.test(resis,peso,df,cme, group=TRUE)
dunca## $statistics
## MSerror Df Mean CV
## 8.06 20 15.04 18.87642
##
## $parameters
## test name.t ntr alpha
## Duncan peso 5 0.05
##
## $duncan
## Table CriticalRange
## 2 2.949998 3.745452
## 3 3.096506 3.931466
## 4 3.189616 4.049682
## 5 3.254648 4.132249
##
## $means
## resis std r Min Max Q25 Q50 Q75
## 15 9.8 3.346640 5 7 15 7 9 11
## 20 15.4 3.130495 5 12 18 12 17 18
## 25 17.6 2.073644 5 14 19 18 18 19
## 30 21.6 2.607681 5 19 25 19 22 23
## 35 10.8 2.863564 5 7 15 10 11 11
##
## $comparison
## NULL
##
## $groups
## resis groups
## 30 21.6 a
## 25 17.6 b
## 20 15.4 b
## 35 10.8 c
## 15 9.8 c
##
## attr(,"class")
## [1] "group"plot(dunca,horiz=TRUE,las=1,main="duncan")
8.4. Prueba de Student-Newman_keuls (SNK).
newman <- SNK.test(resis,peso,df,cme, group=TRUE)
newman## $statistics
## MSerror Df Mean CV
## 8.06 20 15.04 18.87642
##
## $parameters
## test name.t ntr alpha
## SNK peso 5 0.05
##
## $snk
## Table CriticalRange
## 2 2.949998 3.745452
## 3 3.577935 4.542709
## 4 3.958293 5.025630
## 5 4.231857 5.372958
##
## $means
## resis std r Min Max Q25 Q50 Q75
## 15 9.8 3.346640 5 7 15 7 9 11
## 20 15.4 3.130495 5 12 18 12 17 18
## 25 17.6 2.073644 5 14 19 18 18 19
## 30 21.6 2.607681 5 19 25 19 22 23
## 35 10.8 2.863564 5 7 15 10 11 11
##
## $comparison
## NULL
##
## $groups
## resis groups
## 30 21.6 a
## 25 17.6 b
## 20 15.4 b
## 35 10.8 c
## 15 9.8 c
##
## attr(,"class")
## [1] "group"plot(newman)
8.5 Prueba de scheffé.
sche <- scheffe.test(resis,peso,df,cme, group=TRUE,main="dureza")
sche## $statistics
## MSerror Df F Mean CV Scheffe CriticalDifference
## 8.06 20 2.866081 15.04 18.87642 3.385901 6.079555
##
## $parameters
## test name.t ntr alpha
## Scheffe peso 5 0.05
##
## $means
## resis std r Min Max Q25 Q50 Q75
## 15 9.8 3.346640 5 7 15 7 9 11
## 20 15.4 3.130495 5 12 18 12 17 18
## 25 17.6 2.073644 5 14 19 18 18 19
## 30 21.6 2.607681 5 19 25 19 22 23
## 35 10.8 2.863564 5 7 15 10 11 11
##
## $comparison
## NULL
##
## $groups
## resis groups
## 30 21.6 a
## 25 17.6 ab
## 20 15.4 bc
## 35 10.8 c
## 15 9.8 c
##
## attr(,"class")
## [1] "group"plot(sche)
8.6 Ajuste Bonferroni.
pairwise.t.test(resis,peso, p.adj = "bonferroni")##
## Pairwise comparisons using t tests with pooled SD
##
## data: resis and peso
##
## 15 20 25 30
## 20 0.0541 - - -
## 25 0.0031 1.0000 - -
## 30 2.1e-05 0.0251 0.3754 -
## 35 1.0000 0.1859 0.0116 7.0e-05
##
## P value adjustment method: bonferroni8.7. Ajuste de Holm.
pairwise.t.test(resis,peso, p.adj = "holm")##
## Pairwise comparisons using t tests with pooled SD
##
## data: resis and peso
##
## 15 20 25 30
## 20 0.0270 - - -
## 25 0.0025 0.4694 - -
## 30 2.1e-05 0.0151 0.1126 -
## 35 0.5838 0.0744 0.0081 6.3e-05
##
## P value adjustment method: holm9. Tamaño de muestra.
sigma=3;signif=0.01;potencia=0.95
medias<-c(11,12,15,18,19)
muestra<-power.anova.test(group=5,between.var=var(medias),within.var=sigma^2,sig.level=signif,power=potencia)
muestra##
## Balanced one-way analysis of variance power calculation
##
## groups = 5
## n = 5.812129
## between.var = 12.5
## within.var = 9
## sig.level = 0.01
## power = 0.95
##
## NOTE: n is number in each group10. Curva característica.
sd <- seq(2,10,by=0.5)
nn <- seq(3,30,by=3)
beta <- matrix(nr=length(sd),nc=length(nn))
for (i in 1:length(sd))
beta[i,] <- power.anova.test(groups=5,n=nn,between.var=var(medias),
within.var=sd[i]^2,sig.level=0.01)$power
colnames(beta) <- nn
rownames(beta) <- sd
matplot(sd,beta,type="l",xlab=expression(sigma),ylab=expression(1-beta),lty=1)
grid()
text(rep(10,50),beta[length(sd),],as.character(nn),pos=3)
title(" Curva característica de operación\n para 5 medias de tratamientos")
11. Enfoque de regresión.
regre<-lm(resis~peso)
summary(regre)##
## Call:
## lm(formula = resis ~ peso)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.8 -2.6 0.4 1.4 5.2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.0400 0.5678 26.488 < 2e-16 ***
## pesoc1 -5.4000 0.8978 -6.015 7.01e-06 ***
## pesoc2 -1.2500 0.6348 -1.969 0.062960 .
## pesoc3 -3.9000 0.8978 -4.344 0.000315 ***
## pesoc4 0.0900 0.2839 0.317 0.754520
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.839 on 20 degrees of freedom
## Multiple R-squared: 0.7469, Adjusted R-squared: 0.6963
## F-statistic: 14.76 on 4 and 20 DF, p-value: 9.128e-06Se utilzaron los paquetes
library(agricolae)
library(lsmeans)
library(emmeans)
library(nortest)
library(PMCMRplus)
library(car)
library(outliers)
__