Taller Final
Diseño en bloques y diseño factorial
Dayana Narvaez Gomez- Ana Jaramillo Trujillo- Laura Pacheco Arevalo
library(ggplot2)
library(MASS)
library(plotly)##
## Attaching package: 'plotly'## The following object is masked from 'package:MASS':
##
## select## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutlibrary(agricolae) # Libreria para realizar comparaciones
library(lawstat) # Libreria para realizar prueba de corrida
library(nortest) # Pruebas de normalidad
library(gridExtra)Diseño en Bloque Completamente Aleatorizados - DBCA
Punto 1:
Un químico debe probar el efecto de cuatro agentes químicos sobre la resistencia de un tipo particular de tela. Debido a que podría haber variabilidad de un rollo de tela a otro, el químico decide usar un diseño de bloques aleatorizados, con los rollos de tela considerados como bloques. Selecciona cinco rollos y aplica los cuatro agentes químicos de manera aleatoria a cada rollo. A continuación, se presentan las resistencias a la tensión resultantes. Analizar los datos de este experimento (utilizar 𝛼 =0.05).
Ord.Corrida <- c(1:20)
Agente <- c(1,1,1,1,1,2,2,2,2,2,3,3,3,3,3,4,4,4,4,4)
Rollo <- c(73,68,74,71,67,73,67,75,72,70,75,68,78,73,68,73,71,75,75,69)
df1 <- data.frame(Ord.Corrida = Ord.Corrida,
Agente = factor(Agente),
Rollo = Rollo)a) Realice la gráfica de valores individuales.
p <- ggplot(data = df1, mapping = aes(x = Agente , y = Rollo, col = Agente))
p + geom_point() + ggtitle(label = "Grafica de valores individuales")
En promedio en los 4 agentes quimicos los 5 rollos tienen el mismo efecto en terminos de resistencia a la tension.
El modelo estadistico lineal para el DBCA es:
\[ Y_{ij} = \mu + \tau_{i} + \beta_{j}+\varepsilon_{ij} \] ### b) Plantee el modelo y las hipótesis más adecuadas para el problema.
\[ H_{0}: \tau_{1} = \tau_{2} = ... = \tau_{a} = 0 \\ H_{a}: \tau_{i} \neq 0\: \]
\[ H_{0}: \beta_{1} = \beta_{2} = ... = \beta_{b} = 0 \\ H_{a}: \beta_{j} \neq 0\: \] ### c) Realice el ANOVA y saque las conclusiones apropiadas.
lm.DBCA <- lm(formula = Rollo ~ Agente, data = df1) # Primero se construye un modelo lineal
anova(lm.DBCA) # Se genera el ANOVA para el modelo lineal## Analysis of Variance Table
##
## Response: Rollo
## Df Sum Sq Mean Sq F value Pr(>F)
## Agente 3 12.95 4.3167 0.3863 0.7644
## Residuals 16 178.80 11.1750Con un valor p de 0.7644 y un α=0.05 se rechaza la hipótesis nula para los agentes quimicos, por tanto se infiere que el eecto de cada uno de estos agentes no difiere de manera significativa.
Comparar bloques y tratamientos LSD
LSD <- LSD.test(y = aov(lm.DBCA),
trt = "Rollo",
alpha = 0.05, group = TRUE, console = TRUE)##
## Study: aov(lm.DBCA) ~ "Rollo"
##
## LSD t Test for Rollo
##
## Mean Square Error: 11.175
##
## Rollo, means and individual ( 95 %) CI
##
## Rollo std r LCL UCL Min Max
## 67 67 0 2 61.98899 72.01101 67 67
## 68 68 0 3 63.90853 72.09147 68 68
## 69 69 NA 1 61.91336 76.08664 69 69
## 70 70 NA 1 62.91336 77.08664 70 70
## 71 71 0 2 65.98899 76.01101 71 71
## 72 72 NA 1 64.91336 79.08664 72 72
## 73 73 0 4 69.45668 76.54332 73 73
## 74 74 NA 1 66.91336 81.08664 74 74
## 75 75 0 4 71.45668 78.54332 75 75
## 78 78 NA 1 70.91336 85.08664 78 78
##
## Alpha: 0.05 ; DF Error: 16
## Critical Value of t: 2.119905
##
## Groups according to probability of means differences and alpha level( 0.05 )
##
## Treatments with the same letter are not significantly different.
##
## Rollo groups
## 78 78 a
## 75 75 a
## 74 74 ab
## 73 73 ab
## 72 72 ab
## 71 71 ab
## 70 70 ab
## 69 69 ab
## 68 68 b
## 67 67 btratamientos
LSD <- LSD.test(y = aov(lm.DBCA),
trt = "Agente",
alpha = 0.05, group = TRUE, console = TRUE)##
## Study: aov(lm.DBCA) ~ "Agente"
##
## LSD t Test for Rollo
##
## Mean Square Error: 11.175
##
## Agente, means and individual ( 95 %) CI
##
## Rollo std r LCL UCL Min Max
## 1 70.6 3.049590 5 67.43076 73.76924 67 74
## 2 71.4 3.049590 5 68.23076 74.56924 67 75
## 3 72.4 4.393177 5 69.23076 75.56924 68 78
## 4 72.6 2.607681 5 69.43076 75.76924 69 75
##
## Alpha: 0.05 ; DF Error: 16
## Critical Value of t: 2.119905
##
## least Significant Difference: 4.481983
##
## Treatments with the same letter are not significantly different.
##
## Rollo groups
## 4 72.6 a
## 3 72.4 a
## 2 71.4 a
## 1 70.6 aSegundo punto 2k replicado
Un ingeniero esta interesado en los efectos de la velocidad de corte (A), la geometria de la herramienta (B) y el angulo de corte (C) sobre la vida (en horas) de una maquina herramienta. Se eligen dos niveles de cada factor y se corren tres r?plicas de un dise?o factorial 23. Los resultados fueron los siguientes:
n = 2
k = 3
# Se definen los factores y sus niveles
# Se crea un data frame con todas las combinaciones de los niveles de A , B, C
fac2k <- expand.grid(A = c(-1,1),
B = c(-1,1),
C = c(-1,1))
# Se repita el ciclo for n-1 veces
for(i in 1:(n-1)){
fac2k <- rbind(fac2k,fac2k)
}
# Se agrega la variable de respuesta
fac2k$Y <- c(78,104,119,148,127,113, 164,127)
# Se muestra el diseño
fac2k## A B C Y
## 1 -1 -1 -1 78
## 2 1 -1 -1 104
## 3 -1 1 -1 119
## 4 1 1 -1 148
## 5 -1 -1 1 127
## 6 1 -1 1 113
## 7 -1 1 1 164
## 8 1 1 1 127
## 9 -1 -1 -1 78
## 10 1 -1 -1 104
## 11 -1 1 -1 119
## 12 1 1 -1 148
## 13 -1 -1 1 127
## 14 1 -1 1 113
## 15 -1 1 1 164
## 16 1 1 1 127###a) Estimar los efectos de los factores. ¿Que efectos parecen ser grandes? ####Efecto del factor A
ggplot(fac2k, aes(x = A, y = Y)) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Gráfica de efectos del factor A")
Estimacion del efecto de A
EfectoA <- with(data = fac2k, expr = sum(A*Y)/(n*2^(k-1)))
EfectoA## [1] 1####Efecto del factor B
ggplot(fac2k, aes(x = B, y = Y)) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Gr?fica de efectos del factor B")
##### Estimacion del efecto de B
EfectoB <- with(data = fac2k, expr = sum(B*Y)/(n*2^(k-1)))
EfectoB## [1] 34####Efecto del factor C
ggplot(fac2k, aes(x = C, y = Y)) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Gr?fica de efectos del factor C")
##### Estimacion del efecto de C
EfectoC <- with(data = fac2k, expr = sum(C*Y)/(n*2^(k-1)))
EfectoC## [1] 20.5Estimacion de la interaccion AB
ggplot(fac2k, aes(x = A, y = Y, color = factor(B))) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Gráfica de efectos del factor B")
Estimacion de la interaccion AC
ggplot(fac2k, aes(x = A, y = Y, color = factor(C))) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Gráfica de efectos del factor B")
Estimacion de la interaccion BC
ggplot(fac2k, aes(x = B, y = Y, color = factor(C))) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Gráfica de efectos del factor B")
b) Efectuar un análisis de varianza para confirmar los resultados obtenidos en el inciso a
lm2kDep <- lm(Y ~ A + B+ C+A*B+A*C+B*C+A*B*C, data = fac2k)
summary(lm2kDep)## Warning in summary.lm(lm2kDep): essentially perfect fit: summary may be
## unreliable##
## Call:
## lm(formula = Y ~ A + B + C + A * B + A * C + B * C + A * B *
## C, data = fac2k)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.547e-14 -1.205e-15 0.000e+00 1.205e-15 1.547e-14
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.225e+02 2.011e-15 6.092e+16 <2e-16 ***
## A 5.000e-01 2.011e-15 2.486e+14 <2e-16 ***
## B 1.700e+01 2.011e-15 8.454e+15 <2e-16 ***
## C 1.025e+01 2.011e-15 5.097e+15 <2e-16 ***
## A:B -2.500e+00 2.011e-15 -1.243e+15 <2e-16 ***
## A:C -1.325e+01 2.011e-15 -6.589e+15 <2e-16 ***
## B:C -4.250e+00 2.011e-15 -2.113e+15 <2e-16 ***
## A:B:C -3.250e+00 2.011e-15 -1.616e+15 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.044e-15 on 8 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 2.136e+31 on 7 and 8 DF, p-value: < 2.2e-16Con un valor p de 2.2e-16 y un α=0.05 se rechaza la hipótesis nula para los efectos y sus interacciones, interpretando esto como una diferencia significativa en cada uno respecto a la vida util de la maquina sometida al experimento.
c) Escribir y desarrollar todo lo pertinente al modelo de regresión para predecir la vida de la herramienta (en horas) con base en los resultados de este experimento.
lm2k <- lm(Y ~ A + B+A*B+C+A*C+B*C+A*B*C, data = fac2k)
summary(lm2k)## Warning in summary.lm(lm2k): essentially perfect fit: summary may be unreliable##
## Call:
## lm(formula = Y ~ A + B + A * B + C + A * C + B * C + A * B *
## C, data = fac2k)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.547e-14 -1.205e-15 0.000e+00 1.205e-15 1.547e-14
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.225e+02 2.011e-15 6.092e+16 <2e-16 ***
## A 5.000e-01 2.011e-15 2.486e+14 <2e-16 ***
## B 1.700e+01 2.011e-15 8.454e+15 <2e-16 ***
## C 1.025e+01 2.011e-15 5.097e+15 <2e-16 ***
## A:B -2.500e+00 2.011e-15 -1.243e+15 <2e-16 ***
## A:C -1.325e+01 2.011e-15 -6.589e+15 <2e-16 ***
## B:C -4.250e+00 2.011e-15 -2.113e+15 <2e-16 ***
## A:B:C -3.250e+00 2.011e-15 -1.616e+15 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.044e-15 on 8 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 2.136e+31 on 7 and 8 DF, p-value: < 2.2e-16anova(lm2k)## Warning in anova.lm(lm2k): ANOVA F-tests on an essentially perfect fit are
## unreliable## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## A 1 4 4 6.1820e+28 < 2.2e-16 ***
## B 1 4624 4624 7.1464e+31 < 2.2e-16 ***
## C 1 1681 1681 2.5980e+31 < 2.2e-16 ***
## A:B 1 100 100 1.5455e+30 < 2.2e-16 ***
## A:C 1 2809 2809 4.3413e+31 < 2.2e-16 ***
## B:C 1 289 289 4.4665e+30 < 2.2e-16 ***
## A:B:C 1 169 169 2.6119e+30 < 2.2e-16 ***
## Residuals 8 0 0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Análisis de los coeficientes de regresion
X0 <- c(1,1,1,1,1,1,1,1)
X1 <- c(-1,1,-1,1,-1,1,-1,1)
X2 <- c(-1,-1,1,1,-1,-1,1,1)
X3 <-c(-1,-1,-1,-1,1,1,1,1)
X12 <- X1*X2
X13 <- X1*X3
X23 <- X2*X3
X123 <-X1* X2*X3
X <- cbind(X0,X1,X2,X12,X13,X23,X123)
X## X0 X1 X2 X12 X13 X23 X123
## [1,] 1 -1 -1 1 1 1 -1
## [2,] 1 1 -1 -1 -1 1 1
## [3,] 1 -1 1 -1 1 -1 1
## [4,] 1 1 1 1 -1 -1 -1
## [5,] 1 -1 -1 1 -1 -1 1
## [6,] 1 1 -1 -1 1 -1 -1
## [7,] 1 -1 1 -1 -1 1 -1
## [8,] 1 1 1 1 1 1 1Luego se calcula X’X:
t(X)%*%X## X0 X1 X2 X12 X13 X23 X123
## X0 8 0 0 0 0 0 0
## X1 0 8 0 0 0 0 0
## X2 0 0 8 0 0 0 0
## X12 0 0 0 8 0 0 0
## X13 0 0 0 0 8 0 0
## X23 0 0 0 0 0 8 0
## X123 0 0 0 0 0 0 8#####Luego se calcula la inversa de la matriz resultante. Como es una matriz diagonal, esto se puede calcular como el inverso de la matriz diagonal de la matriz:
Cjj <- ginv(t(X)%*%X)
Cjj## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0.125 0.000 0.000 0.000 0.000 0.000 0.000
## [2,] 0.000 0.125 0.000 0.000 0.000 0.000 0.000
## [3,] 0.000 0.000 0.125 0.000 0.000 0.000 0.000
## [4,] 0.000 0.000 0.000 0.125 0.000 0.000 0.000
## [5,] 0.000 0.000 0.000 0.000 0.125 0.000 0.000
## [6,] 0.000 0.000 0.000 0.000 0.000 0.125 0.000
## [7,] 0.000 0.000 0.000 0.000 0.000 0.000 0.125Prediccion con el modelo de regresion
convertir las variables naturales en variables codificadas.
Xcod <- function(X, Bajo, Alto){
Cod <- (X - ((Alto+Bajo)/2))/((Alto-Bajo)/2)
return(Cod)
}A continuación se muestran los valores codificados para diferentes valores de los factores A ,B , C.
XANat <- c(15, 18, 20, 21, 25)
XBNat <- c(1.0, 1.1, 1.5, 1.7, 2.0)
XCNat <- c(1.0, 1.1, 1.5, 1.7, 2.0)
Xpred <- data.frame(
A = Xcod(X = XANat, Bajo = 15, Alto = 25),
B = Xcod(X = XBNat, Bajo = 1, Alto = 2)
)
Xpred## A B
## 1 -1.0 -1.0
## 2 -0.4 -0.8
## 3 0.0 0.0
## 4 0.2 0.4
## 5 1.0 1.0e) Determinar los intervalos de confianza para las betas.
Intervalo de confianza para los B
confint(object = lm2k, level = 0.95)## Warning in summary.lm(object, ...): essentially perfect fit: summary may be
## unreliable## 2.5 % 97.5 %
## (Intercept) 122.50 122.50
## A 0.50 0.50
## B 17.00 17.00
## C 10.25 10.25
## A:B -2.50 -2.50
## A:C -13.25 -13.25
## B:C -4.25 -4.25
## A:B:C -3.25 -3.252k no replicado
Punto 3: Se realizó un experimento en una fábrica de semiconductores en un esfuerzo para incrementar el rendimiento. Se estudiaron cinco factores, cada uno con dos niveles. Los factores (y los niveles) fueron: A=ajuste de apertura (pequeña, grande), B=tiempo de exposición (20% abajo del nominal, 20% arriba del nominal), C=tiempo de desarrollo (30 segundos, 45 segundos), D= tamaño de la máscara (pequeña, grande) y E= tiempo de grabado (14.5 minutos, 15.5 minutos). Se corrió el diseño 25 no replicado que se muestra a continuación:
# Se define el diseño factorial con k = 4 factores y sin replicas
fac2k <- expand.grid(
A = c(-1,1),
B = c(-1,1),
C = c(-1,1),
D = c(-1,1),
E=c(-1,1)
)
Rendimiento <- c(7,9,34,55,16,20,40,60,8,10,32,50,18,21,44,61,8,12,35,52,15,22,45,65,6,10,30,53,15,20,41,63)
# Se agregan los datos de la variable de respuesta
fac2k$Y <- Rendimiento
# Se muestra un resumen del diseño experimental
fac2k## A B C D E Y
## 1 -1 -1 -1 -1 -1 7
## 2 1 -1 -1 -1 -1 9
## 3 -1 1 -1 -1 -1 34
## 4 1 1 -1 -1 -1 55
## 5 -1 -1 1 -1 -1 16
## 6 1 -1 1 -1 -1 20
## 7 -1 1 1 -1 -1 40
## 8 1 1 1 -1 -1 60
## 9 -1 -1 -1 1 -1 8
## 10 1 -1 -1 1 -1 10
## 11 -1 1 -1 1 -1 32
## 12 1 1 -1 1 -1 50
## 13 -1 -1 1 1 -1 18
## 14 1 -1 1 1 -1 21
## 15 -1 1 1 1 -1 44
## 16 1 1 1 1 -1 61
## 17 -1 -1 -1 -1 1 8
## 18 1 -1 -1 -1 1 12
## 19 -1 1 -1 -1 1 35
## 20 1 1 -1 -1 1 52
## 21 -1 -1 1 -1 1 15
## 22 1 -1 1 -1 1 22
## 23 -1 1 1 -1 1 45
## 24 1 1 1 -1 1 65
## 25 -1 -1 -1 1 1 6
## 26 1 -1 -1 1 1 10
## 27 -1 1 -1 1 1 30
## 28 1 1 -1 1 1 53
## 29 -1 -1 1 1 1 15
## 30 1 -1 1 1 1 20
## 31 -1 1 1 1 1 41
## 32 1 1 1 1 1 63a) Realice un análisis de las gráficas de los efectos principales y las interacciones.
lmDoe2k1 <- lm(formula = Y ~ A*B*C*D*E, data = fac2k)
summary(lmDoe2k1)##
## Call:
## lm(formula = Y ~ A * B * C * D * E, data = fac2k)
##
## Residuals:
## ALL 32 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.53125 NA NA NA
## A 5.90625 NA NA NA
## B 16.96875 NA NA NA
## C 4.84375 NA NA NA
## D -0.40625 NA NA NA
## E 0.21875 NA NA NA
## A:B 3.96875 NA NA NA
## A:C 0.21875 NA NA NA
## B:C 0.03125 NA NA NA
## A:D -0.03125 NA NA NA
## B:D -0.34375 NA NA NA
## C:D 0.40625 NA NA NA
## A:E 0.46875 NA NA NA
## B:E 0.28125 NA NA NA
## C:E 0.15625 NA NA NA
## D:E -0.59375 NA NA NA
## A:B:C -0.21875 NA NA NA
## A:B:D 0.15625 NA NA NA
## A:C:D -0.21875 NA NA NA
## B:C:D 0.21875 NA NA NA
## A:B:E -0.09375 NA NA NA
## A:C:E 0.15625 NA NA NA
## B:C:E 0.46875 NA NA NA
## A:D:E 0.40625 NA NA NA
## B:D:E 0.09375 NA NA NA
## C:D:E -0.40625 NA NA NA
## A:B:C:D -0.03125 NA NA NA
## A:B:C:E 0.09375 NA NA NA
## A:B:D:E 0.46875 NA NA NA
## A:C:D:E -0.15625 NA NA NA
## B:C:D:E -0.46875 NA NA NA
## A:B:C:D:E -0.09375 NA NA NA
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 31 and 0 DF, p-value: NALos efectos e interacciones se pueden calcular multiplicando los betas estimados por 2.
# Se extraen los coeficientes del modelo y se multiplican por 2
Terminos <- lmDoe2k1$coefficients[2:32]*2
# Se guarda la información en un data frame
Terminos <- data.frame(
Terminos = factor(names(Terminos)),
Estimado = Terminos,
Signo = factor(ifelse(Terminos >= 0, "Positivo", "Negativo"))
)
# Se reordean los terminos (niveles internos del data frame) segun su valor estimado
Terminos$Terminos <- reorder(Terminos$Terminos, abs(Terminos$Estimado))
Terminos## Terminos Estimado Signo
## A A 11.8125 Positivo
## B B 33.9375 Positivo
## C C 9.6875 Positivo
## D D -0.8125 Negativo
## E E 0.4375 Positivo
## A:B A:B 7.9375 Positivo
## A:C A:C 0.4375 Positivo
## B:C B:C 0.0625 Positivo
## A:D A:D -0.0625 Negativo
## B:D B:D -0.6875 Negativo
## C:D C:D 0.8125 Positivo
## A:E A:E 0.9375 Positivo
## B:E B:E 0.5625 Positivo
## C:E C:E 0.3125 Positivo
## D:E D:E -1.1875 Negativo
## A:B:C A:B:C -0.4375 Negativo
## A:B:D A:B:D 0.3125 Positivo
## A:C:D A:C:D -0.4375 Negativo
## B:C:D B:C:D 0.4375 Positivo
## A:B:E A:B:E -0.1875 Negativo
## A:C:E A:C:E 0.3125 Positivo
## B:C:E B:C:E 0.9375 Positivo
## A:D:E A:D:E 0.8125 Positivo
## B:D:E B:D:E 0.1875 Positivo
## C:D:E C:D:E -0.8125 Negativo
## A:B:C:D A:B:C:D -0.0625 Negativo
## A:B:C:E A:B:C:E 0.1875 Positivo
## A:B:D:E A:B:D:E 0.9375 Positivo
## A:C:D:E A:C:D:E -0.3125 Negativo
## B:C:D:E B:C:D:E -0.9375 Negativo
## A:B:C:D:E A:B:C:D:E -0.1875 NegativoGrafica pareto de los efectos
ggplot(Terminos, aes(x = Terminos, y = abs(Estimado), fill = Signo)) +
geom_bar(stat = "identity") +
coord_flip() +
ggtitle("Gráfica de pareto de los efectos e interacciones")
#### Grafica de normalidad para los efectos
ggplot(Terminos, aes(sample = Estimado, label = Terminos)) +
stat_qq() +
stat_qq_line() +
ggtitle("Gráfica Normal de los Efectos")
Varianza
d) Escribir el modelo de regresión que relaciones el rendimiento con las variables significativas del proceso.
Modelo de regresion
lmDoe2k2 <- lm(formula = Y ~ A + B + C + A*B , data = fac2k)
summary(lmDoe2k2)##
## Call:
## lm(formula = Y ~ A + B + C + A * B, data = fac2k)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7812 -1.2812 -0.3438 1.2500 2.7812
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.5312 0.3021 101.07 < 2e-16 ***
## A 5.9063 0.3021 19.55 < 2e-16 ***
## B 16.9688 0.3021 56.17 < 2e-16 ***
## C 4.8438 0.3021 16.03 2.53e-15 ***
## A:B 3.9688 0.3021 13.14 3.04e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.709 on 27 degrees of freedom
## Multiple R-squared: 0.9932, Adjusted R-squared: 0.9922
## F-statistic: 991.8 on 4 and 27 DF, p-value: < 2.2e-16efectos e interacciones significativas
p1 <- ggplot(fac2k, aes(x = A, y = Y)) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Efecto de A:Ajuste de apertura")
p2 <- ggplot(fac2k, aes(x = B, y = Y)) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Efecto de B: Tiempo de exposicion")
p3 <- ggplot(fac2k, aes(x = C, y = Y)) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Efecto de C: Tiempo de desarrollo")
p4 <- ggplot(fac2k, aes(x = D, y = Y)) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Efecto de D: Tamaño de la mascara")
p5<-ggplot(fac2k, aes(x = E, y = Y)) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Efecto de E: Tiempo de grabado")
grid.arrange(p1, p2, p3, p4,p5, nrow = 2)
p5 <- ggplot(fac2k, aes(x = A, y = Y, color = factor(B))) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Interacción AD")
grid.arrange(p5, nrow = 1)
e) Interpretar cualquier interacción significativa.
Segun los resultados obtenidos, cada efecto e interacicion son significativos para el experimento.
3(k)
Punto 4:
Un investigador médico estudia el efecto de la lidocaína sobre el nivel de enzimas en el músculo cardiaco de perros Beagle. En el experimento se usan tres marcas comerciales de lidocaína (A), tres dosis (B) y tres perros (C). y se corren dos réplicas de un diseño factorial 33. Los niveles de enzimas observados se presentan a continuación.
A <- c(-1,0,1)
B <- c(-1,0,1)
C <- c(-1,0,1)
Base <- expand.grid(A = A,
B = B,
C = C)
for(i in 1:(n-1)){
Doe3k <- rbind(Base,Base)
}
Doe3k$adhesion <- c(96, 94, 101, 85, 95, 108, 84, 95, 105, 84, 99, 106, 84, 98, 114, 83, 97, 100, 85, 98, 98, 86, 97, 109, 81, 93, 106, 84, 95, 105, 80, 93, 110, 83, 92, 102, 85, 97, 104, 82, 99, 102, 80, 96, 111, 86, 90, 103, 84, 95, 100, 79, 93, 108 )
Doe3k## A B C adhesion
## 1 -1 -1 -1 96
## 2 0 -1 -1 94
## 3 1 -1 -1 101
## 4 -1 0 -1 85
## 5 0 0 -1 95
## 6 1 0 -1 108
## 7 -1 1 -1 84
## 8 0 1 -1 95
## 9 1 1 -1 105
## 10 -1 -1 0 84
## 11 0 -1 0 99
## 12 1 -1 0 106
## 13 -1 0 0 84
## 14 0 0 0 98
## 15 1 0 0 114
## 16 -1 1 0 83
## 17 0 1 0 97
## 18 1 1 0 100
## 19 -1 -1 1 85
## 20 0 -1 1 98
## 21 1 -1 1 98
## 22 -1 0 1 86
## 23 0 0 1 97
## 24 1 0 1 109
## 25 -1 1 1 81
## 26 0 1 1 93
## 27 1 1 1 106
## 28 -1 -1 -1 84
## 29 0 -1 -1 95
## 30 1 -1 -1 105
## 31 -1 0 -1 80
## 32 0 0 -1 93
## 33 1 0 -1 110
## 34 -1 1 -1 83
## 35 0 1 -1 92
## 36 1 1 -1 102
## 37 -1 -1 0 85
## 38 0 -1 0 97
## 39 1 -1 0 104
## 40 -1 0 0 82
## 41 0 0 0 99
## 42 1 0 0 102
## 43 -1 1 0 80
## 44 0 1 0 96
## 45 1 1 0 111
## 46 -1 -1 1 86
## 47 0 -1 1 90
## 48 1 -1 1 103
## 49 -1 0 1 84
## 50 0 0 1 95
## 51 1 0 1 100
## 52 -1 1 1 79
## 53 0 1 1 93
## 54 1 1 1 108c) Efectos cuadráticos.
Grafica de efectos
p1 <- ggplot(Doe3k, aes(x = A, y = adhesion)) +
geom_point() +
stat_summary(geom = "line", fun.y = mean, color = "red") +
ggtitle("Efecto de la marca 1")
p2 <- ggplot(Doe3k, aes(x = B, y = adhesion)) +
geom_point() +
stat_summary(geom = "line", fun.y = mean, color = "red") +
ggtitle("Efecto de la marca 2")
p3 <- ggplot(Doe3k, aes(x = C, y = adhesion)) +
geom_point() +
stat_summary(geom = "line", fun.y = mean, color = "red") +
ggtitle("Efecto de la marca 3")
grid.arrange(p1, p2,p3, nrow = 1)
#### Grafica de interacciones}
ggplot(Doe3k, aes(x = A, y = adhesion, color = factor(B))) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Interacción AB")
ggplot(Doe3k, aes(x = A, y = adhesion, color = factor(C))) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Interacción AC")
ggplot(Doe3k, aes(x = B, y = adhesion, color = factor(C))) +
geom_point() +
stat_summary(geom = "line", fun.y = mean) +
ggtitle("Interacción BC")
Analisis de regresión
lmDoe3k1 <- lm(formula = adhesion ~ A + B + C + A*B + A*C + B* C + I(A^2) + I(B^2) + I(C^2), data = Doe3k)
summary(lmDoe3k1)##
## Call:
## lm(formula = adhesion ~ A + B + C + A * B + A * C + B * C + I(A^2) +
## I(B^2) + I(C^2), data = Doe3k)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.7824 -1.8588 0.1968 1.8391 9.1898
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 96.9630 1.2150 79.802 < 2e-16 ***
## A 10.5833 0.5625 18.816 < 2e-16 ***
## B -0.6111 0.5625 -1.087 0.28317
## C -0.4444 0.5625 -0.790 0.43366
## I(A^2) -0.8056 0.9742 -0.827 0.41276
## I(B^2) -1.2222 0.9742 -1.255 0.21625
## I(C^2) -1.2222 0.9742 -1.255 0.21625
## A:B 1.8750 0.6889 2.722 0.00927 **
## A:C 0.1667 0.6889 0.242 0.80995
## B:C 0.5833 0.6889 0.847 0.40169
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.375 on 44 degrees of freedom
## Multiple R-squared: 0.8932, Adjusted R-squared: 0.8713
## F-statistic: 40.88 on 9 and 44 DF, p-value: < 2.2e-16# Modelo Depurado
lmDoe3k2 <- lm(formula = adhesion ~ A + A*B, data = Doe3k)
summary(lmDoe3k2)##
## Call:
## lm(formula = adhesion ~ A + A * B, data = Doe3k)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.644 -1.779 -0.213 1.867 9.301
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 94.7963 0.4557 208.023 < 2e-16 ***
## A 10.5833 0.5581 18.963 < 2e-16 ***
## B -0.6111 0.5581 -1.095 0.27878
## A:B 1.8750 0.6836 2.743 0.00843 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.349 on 50 degrees of freedom
## Multiple R-squared: 0.8805, Adjusted R-squared: 0.8733
## F-statistic: 122.8 on 3 and 50 DF, p-value: < 2.2e-16d) Contornos y superficie de respuesta
Grafica de contornos
Nota: para la gráfica de contorno y superficie de respuesta se debe crear la variable z, que basicamente es la estimación de la variable de respuesta en función de x y y, esto es: z=f(x,y)
XA <- seq(-1,1,0.1)
XB <- seq(-1,1,0.1)
YP <- outer(X = XA,
Y = XB,
FUN = function(X,Y)94.7963 + 10.5833*X + 1.8750 *X*Y)
p4 <- plot_ly(x = XA, y = XB, z = t(YP)) %>%
add_surface(
contours = list(
z = list(
show=TRUE,
usecolormap=TRUE,
highlightcolor="#ff0000",
project=list(z=TRUE)
)
))
p4<div id="htmlwidget-f422f2ca10fe85443c37" style="width:672px;height:480px;" class="plotly html-widget">