## disponibilidad de ggplot2
if(!require(ggplot2)) install.packages('ggplot2')
require(ggplot2)
## disponibilidad de ggridges
if(!require(ggridges)) install.packages('ggridges')
require(ggridges)
## disponibilidad de knitr
if(!require(knitr)) install.packages('knitr')
require(knitr)
## disponibilidad de lmtest
if(!require(lmtest)) install.packages('lmtest')
require(lmtest)
## disponibilidad de lme4 para componentes de varianza
if(!require(lme4)) install.packages('lme4')
require(lme4)
## availability of pwr
if(!require(pwr)) install.packages('pwr')
require(pwr)
## check availability of the 'beanplot' package
if(!require(beanplot)) install.packages('beanplot')
require(beanplot)
## check availability of 'car' package
if(!require(car)) install.packages('car')
require(car)
if(!require(effectsize)) install.packages('effectsize')
require(effectsize)
if(!require(broom)) install.packages('broom')
require(broom)Diseño y Análisis de Experimentos
Resultados Trabajo #1
##Paquetes
3.1
An experimenter has conducted a single-factor experiment with four levels of the factor, and each factor level has been replicated six times. Find the \(P-value\) if the \(F-statistic\) is \({F_0=3.26}\)`
Rta/. Datos:
\({a=4}\)
\({n=6}\)
\(F_{calc}=3.26\)
\({gl_1= a - 1 = 3}\)
\({gl_2= N - a = 24-4=20}\)
El valor de \({P_{value}~(F=3.26,3,20)}\) es:
F <- 3.26
gl1 <- 3
gl2 <- 20
Pvalue <- pf(F,gl1,gl2,lower.tail = FALSE )
Pvalue[1] 0.04300116Conclusión: el factor controlado explica los efectos sobre la respuesta mejor que los efectos de los errores sobre la misma. Suponiendo un nivel de significancia de \(\alpha=0.05\) Rechazamos \(H_0\), por lo que conlcuimos que existe almenos un nivel del factor controlado para el que la respuesta son diferentes a las demás. En otras palabras, esto indica que la variable estudiada afecta la respuesta.
3.8
The ANOVA of an experiment is show below. Fill in the blanks. Is there any way to approach the observed power? Explain.
MSSA <- 246.93
SST <- 1174.24
SSE <- 186.53
gl2 <- 25
glt <- 29
Z <- 186.53/25
X <- glt-gl2
N <- glt+1
a <- X+1
n <- N/a
Y <- SST-SSE
SSA <- SST - SSE
W <- MSSA/Z
V <- pf(W,X,gl2, lower.tail = FALSE)\(X=a-1= (N-1)-(N-a) = 29 - 25 =\) 4
\(a=X+1=\) 5
\(N-1 = 29\)
\(N=\) 30
\(n=N/a=30/5=\) 6
\(SSA =SST - SSE = 1174.24 -186.53\)
\(SSA =\) 987.71
\(Y=SSE\) 987.71
\(Z=MSS_E=SSE/gl2\)
\(Z=\) 7.4612
\(W=F_{calc.}=MSS_A/MSS_E=\) 33.0952126
El valor de \(P_{value}\) para \(F=\) 33.0952126, con 4 y 25 grados de libertad es: 1.1846546^{-9}
Para obtener una apróximación de la potencia observada debemos utilizar la función pwr.anova.test(K=,n=,f=,Sig.level=), con \(K=\) 5,\(n=\) 6, \(f^2=R^2/(1-R^2)\) y un nivel de significancia asumido utilizado en el diseño de la prueba de \(\alpha=0.05\)
R2 <- SSA/SST
f <- ((R2)/(1-(R2)))^(0.5)
pwr.anova.test(a,n,f,sig.level = 0.05)
Balanced one-way analysis of variance power calculation
k = 5
n = 6
f = 2.301126
sig.level = 0.05
power = 1
NOTE: n is number in each groupConclusión: el factor controlado explica los efectos sobre la respuesta mejor que los efectos de los errores sobre la misma. Para casi cualquier nivel de significancia que asumamos superior a \(\alpha=\) 1.1846546^{-9} Rechazamos \(H_0\), por lo que concluimos que existe almenos un nivel del factor controlado para el que la respuesta son diferentes a las demás. En otras palabras, esto indica que la variable estudiada afecta la respuesta. Adicionalmente, podemos encontrar el efecto que buscamos de los factores sobre la respuesta en este diseño de experimento con una probabilidad de 1 contemplando una relación de \(f~(R) ->\) \(f^2=R^2/(1-R^2)\) y un nivel de significancia asumido en el diseño de la prueba de \(\alpha=0.05\).
3.28
Four different designs for a digital computer circuit are being studied to compare the amount of noise present. The following data have been obtained
Para el análisis del experimento se seguirán los siguientes pasos:
- Lectura de datos
- Análisis gráfico
- Construcción de la tabla ANOVA
- Construcción del modelo de Regresión Lineal Simple/Múltiple
- Pruebas simultáneas
- Validación de supuestos
- Conclusiones
A continuación se desarolla cada paso.
Lectura de datos
Los datos con los que trabajaremos son los siguientes:
## read the data set
datos <- read.table('HomeWork2_files/3.28 Ejercicio.txt',
header = TRUE)
## convert 'Circuit_Design' to factor "Convertir la columna Circuit_Design rn un factor"
datos$Circuit_Design <- as.factor(datos$Circuit_Design)
## portion of the data "mostrar una parte de los datos"
head(datos, 10) Circuit_Design Noise_observed
1 1 19
2 1 20
3 1 19
4 1 30
5 1 8
6 2 80
7 2 61
8 2 73
9 2 56
10 2 80Análisis gráfico
Una vez leemos los datos en R, el primer paso es realizar un análisis gráfico de los resultados del experimento.
## boxplots for main effects
boxplot(Noise_observed ~ Circuit_Design, data = datos, las = 1, col = 2:5, xlab = 'Circuit_Design', ylab = 'Noise_observed')
#"Noise_observed ~ Circuit_Design" y en funcion de x= y(x)
# insertar una línea horizontal en la gráfica
abline(h= mean(datos$Noise_observed),col="orange",lwd =3)
En nuestro caso, el gráfico puede obtenerse haciendo:
## nice colors
cols <- c("#0080ff", "#ff00ff", "darkgreen", "#ff0000",
"orange", "#00ff00", "brown")
with(datos,
beanplot(Noise_observed ~ Circuit_Design, las = 1, horizontal = FALSE,
col = c('black', cols[1], 3, cols[2]),
log = "", what = c(1, 1, 1, 1), cutmin = 100, cutmax = 200,
cex.axis = 1,
lwd = 2, main = "",
xlab = "Circuit_Design", ylab = 'Noise_observed'))
mean(datos$Noise_observed)[1] 51.4En el caso de visualizar los datos utilizando un violin plot podemos proceder de la siguiente manera:
## violin plot
ggplot(datos, aes(x = Circuit_Design, y = Noise_observed, fill = Circuit_Design)) +
geom_violin(alpha = .5) +
geom_point() +
geom_jitter() +
theme(legend.position = "none") +
ylab("Noise_observed") +
xlab("Circuit_Design") +
theme_minimal()
Ahora, procedamos a calcular algunas medidas de tendencia central, posición y dispersión:
## measures using tapply
with(datos, tapply(Noise_observed, Circuit_Design, mean)) 1 2 3 4
19.2 70.0 36.6 79.8 with(datos, tapply(Noise_observed, Circuit_Design, sd)) 1 2 3 4
7.79102 11.02270 11.58879 20.51097 with(datos, tapply(Noise_observed, Circuit_Design, summary))$`1`
Min. 1st Qu. Median Mean 3rd Qu. Max.
8.0 19.0 19.0 19.2 20.0 30.0
$`2`
Min. 1st Qu. Median Mean 3rd Qu. Max.
56 61 73 70 80 80
$`3`
Min. 1st Qu. Median Mean 3rd Qu. Max.
25.0 26.0 35.0 36.6 47.0 50.0
$`4`
Min. 1st Qu. Median Mean 3rd Qu. Max.
46.0 78.0 83.0 79.8 95.0 97.0 Construcción de la tabla ANOVA
En este paso hacemos uso de las funciones lm() y anova():
fit <- lm(Noise_observed ~ Circuit_Design, data = datos)
anova(fit)Analysis of Variance Table
Response: Noise_observed
Df Sum Sq Mean Sq F value Pr(>F)
Circuit_Design 3 12042.0 4014.0 21.78 6.797e-06 ***
Residuals 16 2948.8 184.3
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Rta. a.
a Is the same amount of noise present for all four designs? Use \({\alpha}=\) 0.05.
A partir del análisis gráfico y la tabla ANOVA se puede evidenciar que no hay evidencia suficiente, con un nivel de significancia de \({\alpha}=\) 0.05, que nos lleve a concluir que las medias son iguales. \(Pr(>F)=\) 6.797e-06, por lo cual se rechaza la hipótesis núla, y se puede afirmar que existe al menos un Diseño de Circuito que en el que los efectos sobre el Ruido observado es diferente a las demás.
Validación de Supuestos
En el 1FEF, al igual que en otros DOE, se tienen los siguientes supuestos sobre el error aleatorio, \(\epsilon\):
- \(\mu_{\epsilon} = 0\), donde \(\mu_{\epsilon}\) es la media de \(\mathbf{\epsilon}\),
- \(\epsilon_1, \epsilon_2, \ldots, \epsilon_n\) tiene varianza constante \(\sigma^2\),
- \(\mathbf{\epsilon} \sim N(0, \sigma^2)\),
- \(\epsilon_i\) y \(\epsilon_j\) son independientes para \(i\neq j\).
En R dicha validación se realiza sobre los residuales del modelo ajustado. En nuestro caso, este modelo está contenido en el objeto fitdidisplays.
Inicialmente calculamos los residuales del modelo ajustado:
## residuals
r <- rstudent(fit)El supuesto de media cero podemos validarlo usando una prueba \(t\) de Student, aunque no es necesario:
## media cero
t.test(r, mu = 0)
One Sample t-test
data: r
t = -0.1731, df = 19, p-value = 0.8644
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
-0.5978119 0.5064808
sample estimates:
mean of x
-0.04566556 El supuesto de Normalidad puede validarse de varias maneras. En este curso usaremos una prueba formal como la Shapiro-Wilk y un método gráfico basado en el QQ-plot.
## Shapiro-Wilk normality test
shapiro.test(r)
Shapiro-Wilk normality test
data: r
W = 0.8682, p-value = 0.01093Puesto que el valor \(p\) de la prueba es \(<0.05\), rechazamos \(H_0\) y concluimos que los errores del modelo ajustado no siguen una distribución Normal.
El QQ-plot puede construirse de la siguiente forma:
## QQ-plot
qqnorm(r, las = 1)
qqline(r, col = 1, lty = 2)
Cuando se cumple el supuesto de normalidad de los residuales, se espera que los puntos en el gráfico se encuentren alrededor de la línea punteada. En este caso en particular, el supuesto no parece ser razonable, alejados significativamente de la línea punteada y con mayor densidad hacia un costado de la línea que del otro.
Rta.b
Analyze the residuals from this experiment. Are the analysis of variance assumptions satisfied?.
El supuesto de varianza constante puede validarse, esencialmente, utilizando dos aproximaciones: una prueba formal y métodos gráficos. En en este caso usaré la prueba de Breusch-Pagan implementada a través de la función ncvTest() del paquete car.
## constant variance using Breusch-Pagan
car:::ncvTest(fit)Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 3.220545, Df = 1, p = 0.072719Dado que \(P_{value}>\) 0.05 no se rechaza \(H_0\) por lo cual se puede concluir que no hay evidencia significativa que nos sugiera la existencia de varianzas diferentes.
Gráficamente podemos proceder de la siguiente manera:
## constant variance (graphical method)
# r
par(mfrow = c(1, 2))
plot(r, las = 1, ylab = 'Residual', xlab = "Observation #",
main = 'Residuals')
abline(h = 0, col = 2, lty = 2)
# r vs. display
boxplot(r ~ Circuit_Design, col = cols,
las = 1,
ylab = 'Residual',
xlab = "Circuit_Design (W)",
main = "Circuit_Design entry",
data = datos)
abline(h = 0, col = 2, lty = 2)
La conclusión general es que los residuales del modelo ajustado tienen varianza constante.
El supuesto de independencia puede validarse, al igual que los supuestos de varianza constante y de normalidad, puede validarse de manera gráfica o haciendo una prueba formal. Gráficamente se utiliza la función de autocorrelación o ACF en inglés. En R procedemos de la siguiente manera para construir la ACF:
## independence using ACF
acf(r, las = 1, main = "")
mtext("ACF", side = 3, line = .2)
La idea general es que, a excepción de la primera barra vertical, ninguna de las barras verticales debe estar por fuera de las líneas horizontales de color azul.
Para probar independencia de manera formal, usamos la prueba de Durbin-Watson. En R esta prueba se encuentra implementada en la función durbinWatsonTest del paquete car:
## independence test using Durbin-Watson
set.seed(123)
durbinWatsonTest(fit) lag Autocorrelation D-W Statistic p-value
1 -0.2969208 2.493502 0.652
Alternative hypothesis: rho != 0Puesto que el valor \(p\) de la prueba es \(>0.05\), concluimos que los residuales del modelo ajustado son independientes.
Rta c.
Which circuit design would you select for use? Low noise is best.
Estos resultados indican que los errores del modelo ajustado son independientes, siguen una distribución normal y tienen varianza constante. Por lo tanto, el modelo y las conclusiones que se deriven de él son válidas.
Por lo tanto, podemos sugerir con un nivel de confiabilidad del \({(1-\alpha)=}\) ´95%´el Diseño de circuito que mejor condiciones de ruido presenta (Niveles más bajo), y con el cual se sugiere trabajar es el Circuit_Design1
Construcción del modelo de Regresión
Aunque el resultado anterior es importante, por lo general interesa determinar la magnitud del efecto de cada nivel de Circuit_Design sobre Noise_Observed. Esto es posible a partir del objeto fit. Para ello utilizamos los estimadores de máxima verosimilitud obtenidos al emplear el método de mínimos cuadrados ordinarios en el modelo de Regresión Lineal Simple:
## SLR with categorical predictor
summary(fit)
Call:
lm(formula = Noise_observed ~ Circuit_Design, data = datos)
Residuals:
Min 1Q Median 3Q Max
-33.8 -9.4 0.3 10.1 17.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.200 6.071 3.162 0.00604 **
Circuit_Design2 50.800 8.586 5.917 2.17e-05 ***
Circuit_Design3 17.400 8.586 2.027 0.05971 .
Circuit_Design4 60.600 8.586 7.058 2.71e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.58 on 16 degrees of freedom
Multiple R-squared: 0.8033, Adjusted R-squared: 0.7664
F-statistic: 21.78 on 3 and 16 DF, p-value: 6.797e-06El modelo que nos indica explica los efectos de los Circuit_designsobre el Noise_Observedestá dado por \(Y=19.2+50.8*Circuit_{Design2}+17.4*Circuit_{Design3}+60.6*Circuit_{Design4}+13.58\)
Sin embargo es posible que el nivel Circuit_Design3 pueda ser excluido del modelo, debido a que para un valor de significancia del 5% no resulta significativo, para lo cual se vuelve a calcular los estimadores de los efectos excluyendo el nivel 3:
# Excluye los datos con Circuit_Design == "Circuit_Design3"
datos_filtrados <- datos[datos$Circuit_Design != "3", ]
# Vuelve a estimar el modelo con los datos filtrados
fit <- lm(Noise_observed ~ Circuit_Design, data = datos_filtrados)
summary(fit)
Call:
lm(formula = Noise_observed ~ Circuit_Design, data = datos_filtrados)
Residuals:
Min 1Q Median 3Q Max
-33.8 -5.4 0.8 10.0 17.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19.200 6.340 3.028 0.010496 *
Circuit_Design2 50.800 8.966 5.666 0.000105 ***
Circuit_Design4 60.600 8.966 6.759 2.02e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.18 on 12 degrees of freedom
Multiple R-squared: 0.8144, Adjusted R-squared: 0.7835
F-statistic: 26.33 on 2 and 12 DF, p-value: 4.088e-05Se evidencia un ligero incremento en el error residual estandar. Por lo cual el modelo ajustado sin el nivel “Circuit_Design3” es:
\(Y=19.2+50.8*Circuit_{Design2}+60.6*Circuit_{Design4}+14.18\)
Pruebas simultáneas
A partir del modelo de Regresión, podemos concluir que utilizar cualquier diseñor diferente al Circuit_Design 1 incrementará considerablemente el Noise_Observed. En particular, utilizar Circuit_Design4 en lugar del tipo Circuit_Design 1 incrementa en 60.6 el Noise_Observed
Qué podemos decir acerca de Circuit_Design 1 vs. Circuit_Design 3? Para responder a esto utilizamos una comparación ad hoc. Esta comparación puede hacerse utilizando la prueba de mínima diferencia significativa de Tukey, también conocida como Tukey’s Honestly Significant Difference (HSD):
## Honestly Significant Difference (HSD) test
anova <- aov(fit)
(hsd <- TukeyHSD(anova, "Circuit_Design", ordered = TRUE)) Tukey multiple comparisons of means
95% family-wise confidence level
factor levels have been ordered
Fit: aov(formula = fit)
$Circuit_Design
diff lwr upr p adj
2-1 50.8 26.8803 74.7197 0.0002850
4-1 60.6 36.6803 84.5197 0.0000556
4-2 9.8 -14.1197 33.7197 0.5361030Los resultados indican que las diferencias Circuit_Design 1 vs. Circuit_Design 3 no es estadísticamente significativa después de aplicar la correción por múltiples comparaciones, es decir, que la columna p adj es \(>0.05\).
Para visualizar dichas diferencias usamos la función plot sobre el objeto hsd:
## HSD plot
plot(hsd, las = 1)
La conclusión es la misma que obtuvimos al analizar los resultados numéricos.
Calculo de la Potencia Observada
Para obtener una apróximación de la potencia observada debemos utilizar la función pwr.anova.test(K=,n=,f=,Sig.level=), con \(K=\) 5,\(n=\) 6, \(f^2=R^2/(1-R^2)\) y un nivel de significancia asumido utilizado en el diseño de la prueba de \(\alpha=0.05\)
fit <- lm(Noise_observed ~ Circuit_Design, data = datos)
anova(fit)Analysis of Variance Table
Response: Noise_observed
Df Sum Sq Mean Sq F value Pr(>F)
Circuit_Design 3 12042.0 4014.0 21.78 6.797e-06 ***
Residuals 16 2948.8 184.3
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1SSA <- 12042
SST <- 12042+2948.8
R2 <- SSA/SST
f <- ((R2)/(1-(R2)))^(0.5)
pwr.anova.test(k = 4,n = 5,f = f,sig.level = 0.05)
Balanced one-way analysis of variance power calculation
k = 4
n = 5
f = 2.020815
sig.level = 0.05
power = 0.9999999
NOTE: n is number in each groupRta. d
Estimate the effect of the experiment and determine the observed power of the design. What is your conclusion?.
Conclusión: Los efectos de los circuit_designsobre la variable de respuestos Noise_observed son significantes para los niveles 1, 2 y 4 con una con \(\alpha=0.05\). Adicionalmente la potencia observada en la prueba es de 0.999 -> asumiendo \(f^2=R^2/(1-R^2)\) y un nivel de significancia asumido utilizado en el diseño de la prueba de \(\alpha=0.05\). Por lo tanto podemos decir que el diseño de experimento tiene una posibilidad muy alta de encotrar el efecto que se está buscando. Esto está dado principalmente a que el ajuste del modelo \(R^2=0.8033\) indica un buen ajuste.
4.8
A chemist wishes to test the effect of four chemical agents on the strength of a particular type of cloth. Because there might be variability from one bolt to another, the chemist decides to use a RCBD, with the bolts of cloth considered as blocks. She selects five bolts and applies all four chemicals in random order to each bolt. The resulting tensile strengths follow. Analyze the data from this experiment and draw appropriate conclusions. Use \(\alpha=0.05\).
Para el desarrollo de este ejercicio utilizaremos el modelo de DBCA, teniendo como caracteristica que no hay réplicas lo que significa que cada rollo va a sufrir todos los tratamientos de cada uno de los niveles de los químicos. No existe interacción entre los rollos.
Lectura de datos
Para leerlos en R procedemos de la siguiente manera:
datos <- read.table('HomeWork2_files/4.8 Ejercicio.txt',
header = TRUE)
datos$A <- factor(datos$A)
datos$block <- factor(datos$block)
## portion of the data "mostrar una parte de los datos"
head(datos, 10) A block y
1 1 1 73
2 1 2 68
3 1 3 74
4 1 4 71
5 1 5 67
6 2 1 73
7 2 2 67
8 2 3 75
9 2 4 72
10 2 5 70Análisis del Experimento
Visualmente, la relación entre Bolt y el factor A: Chemical puede establecerse haciendo
ggplot(datos, aes(x = A, y = y, fill = A)) +
geom_boxplot() + geom_hline(yintercept=mean(datos$y), linetype="dashed", color ="orange") + xlab('Chemical') + ylab('tensile strengths') + theme_minimal()
Graficamente no se evidencia que el tipo del químico aplicado tenga un efecto significativo sobre la resistencia a la tracción.
En el caso del block se tiene
## boxplot for block using ggplot2
ggplot(datos, aes(x = block, y = y, fill = block)) +
geom_boxplot() +
xlab('bolt') + ylab('tensile strengths') + theme_minimal()
Visualmente se puede sugerir que si existe en al menos uno de los rollos diferencia significativa, en la resistencia a la tracción, si el químico se aplica sobre un rollo u otro.
Construcción de la tabla ANOVA
Ahora construimos la tabla ANOVA incluyendo el término de bloqueo:
## ANOVA including the blocking factor "block"
fit <- aov(y ~ A + block, data = datos)
summary(fit) Df Sum Sq Mean Sq F value Pr(>F)
A 3 12.95 4.32 2.376 0.121
block 4 157.00 39.25 21.606 2.06e-05 ***
Residuals 12 21.80 1.82
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La conclusión de este modelo es que el factor A: Chemical, con un nivel de significancia \(\alpha=0.05\), no influencia la respuesta tensile strenghts, sin embargo el rollo sobre el cual se aplique el químico si tiene efecto sobre la variable de tensile strenghts.
Dado que los rollos sobre los cuales se aplique el químico explican la mayoría de los efectos sobre la respuesta, al no considerarlos como un bloque, se puede concluir que el factor A: Chemical,para cualquier \(\alpha <0.764\) no es dintinguible su efecto sobre la resistencia a la tensión. Lo cual se evidencia a continuación:
## ANOVA with no blocking
fitnoblock <- aov(y ~ A, data = datos)
summary(fitnoblock) Df Sum Sq Mean Sq F value Pr(>F)
A 3 12.95 4.317 0.386 0.764
Residuals 16 178.80 11.175 Cuantificación del efecto
Para cuantificar la magnitud del efecto de cada nivel de A-Chemical sobre Tensile strengths utilizamos un modelo de Regresión Lineal Múltiple:
## MLR for the tubes experiment
fit <- lm(y ~ A + block, data = datos)
summary(fit)
Call:
lm(formula = y ~ A + block, data = datos)
Residuals:
Min 1Q Median 3Q Max
-1.3500 -0.7375 -0.3500 0.7000 1.8500
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 72.3500 0.8524 84.873 < 2e-16 ***
A2 0.8000 0.8524 0.938 0.366507
A3 1.8000 0.8524 2.112 0.056374 .
A4 2.0000 0.8524 2.346 0.036968 *
block2 -5.0000 0.9531 -5.246 0.000206 ***
block3 2.0000 0.9531 2.098 0.057699 .
block4 -0.7500 0.9531 -0.787 0.446584
block5 -5.0000 0.9531 -5.246 0.000206 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.348 on 12 degrees of freedom
Multiple R-squared: 0.8863, Adjusted R-squared: 0.82
F-statistic: 13.36 on 7 and 12 DF, p-value: 8.335e-05A partir de este resultado concluimos que utilizar el químito tipo 4 aumentará en ~2 unidades la variable y: tensile strengths comparado con usar el químico tipo 1
Para calcular el valor \(P\) de esta diferencia, después de ajustar por el factor block, podemos utilizar la HSD de Tukey:
## Tukey HSD
(thsd <- TukeyHSD(aov(fit), "A", ordered = TRUE)) Tukey multiple comparisons of means
95% family-wise confidence level
factor levels have been ordered
Fit: aov(formula = fit)
$A
diff lwr upr p adj
2-1 0.8 -1.7308322 3.330832 0.7852734
3-1 1.8 -0.7308322 4.330832 0.2042593
4-1 2.0 -0.5308322 4.530832 0.1417326
3-2 1.0 -1.5308322 3.530832 0.6540138
4-2 1.2 -1.3308322 3.730832 0.5182726
4-3 0.2 -2.3308322 2.730832 0.9952030Sin embargo, al realizar el comparativo por pares, nos percatamos que no es significativo realizar el cambio debido a que \(P_{adj}>\alpha\), para un \(\alpha=0.05\).
## Tukey's HSD
par(mfrow = c(1, 1), mar = c(5, 6, 3, 1))
plot(thsd, las = 1)
Validación de supuestos
Realizamos la validación de los supuestos sobre de los residuos del modelo:
## cálculo de residuales
r <- residuals(fit)A partir del objeto r, podemos calcular el valor \(P\) de las pruebas de Normalidad, independencia y varianza constante:
## valor p de la prueba de Normalidad de Shapiro-Wilks
shapiro.test(r)$p.value[1] 0.04053571## prueba de varianza constante
ncvTest(fit)$p[1] 0.9464394## valor p de la prueba de independencia
lmtest::dwtest(fit)$p.value[1] 0.7492284Conclusión: Como los valores \(P\) de las pruebas realizadas no son todos \(>0.05\), los errores del modelo ajustado no cumplen con todos los supuestos. Por lo tanto, el modelo y las conclusiones que se deriven de él no son estadisticamente significativas.
Rta
Si se requiriera de una inversión para realizar un cambio de químicos, se recomienda no realizar dicha inversión y utilizar el químico que se habitúa aplicar. Si se está Tener presente que las propiedades mecánicas de un material dependen principalmente de su composición del material, su microestructura y temperatura a la que es expuesto, por lo cual la aplicación de un producto que normalmente solo tiene ingenerencia sobre una capa superficial, no cambiaría significativamente las propiedades.
5.12
Johnson and Leone (Statistics and Experimental Design in Engineering and the Physical Sciences, Wiley, 1977) describe an experiment to investigate warping of cop- per plates. The two factors studied were the temperature and the copper content of the plates. The response variable was a measure of the amount of warping. The data are as follows:
Lectura de datos
## read the data set
datos <- read.table('HomeWork2_files/5.12 Ejercicio.txt',
header = TRUE)
## convert 'Temperature' and 'CopperContent' to factors
datos$Temperature <- as.factor(datos$Temperature)
datos$CopperContent <- as.factor(datos$CopperContent)
## see first 3 rows of the data
head(datos, 16) Temperature CopperContent y
1 50 40 17
2 50 40 20
3 50 60 16
4 50 60 21
5 50 80 24
6 50 80 22
7 50 100 28
8 50 100 27
9 75 40 12
10 75 40 9
11 75 60 18
12 75 60 13
13 75 80 17
14 75 80 12
15 75 100 27
16 75 100 31El número de réplicas puede obtenerse haciendo:
## replicates per combination
with(datos, tapply(y, list(Temperature,CopperContent), length)) 40 60 80 100
50 2 2 2 2
75 2 2 2 2
100 2 2 2 2
125 2 2 2 2El número total de unidades experimentales es
## how many runs?
NROW(datos)[1] 32## boxplots for main effects
par(mfrow = c(1, 2))
boxplot(y ~ Temperature, data = datos, las = 1, col = 2:4,
xlab = 'Temperature', ylab = 'warping')
boxplot(y ~ CopperContent, data = datos, las = 1, col = 5:7,
xlab = 'CopperContent', ylab = 'warping')
Construcción de la tabla ANOVA
Ahora construimos la tabla ANOVA incluyendo el término de interacción:
# Ajustar el modelo ANOVA de dos vías
modelo <- aov(y ~ CopperContent * Temperature, data = datos)
# Obtener el resumen del modelo ANOVA
anova(modelo)Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
CopperContent 3 698.34 232.781 34.3272 3.35e-07 ***
Temperature 3 156.09 52.031 7.6728 0.002127 **
CopperContent:Temperature 9 113.78 12.642 1.8643 0.132748
Residuals 16 108.50 6.781
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Rta. a.
Is there any indication that either factor affects the amount of warping? Is there any interaction between the factors? Use ().
Dado que los valores de \(Pr(>F)\) tanto para CopperContent yTemperatureson \(<0.05\) lo cual indica que tanto los factores de la temperatura como el porcentaje de cobre afectan la deformación
La conclusión de este primer modelo es que la interacción entre CopperContent y Temperature, definida como Temperature:CopperContent en el resultado de anova(fit), es no significativa al 5%.
Puesto que la interacción no es significativa al 5%, debemos remover este término del modelo y recalcularlo:
## ANOVA with no interaction
fit_no_interaction <- aov(y ~ Temperature + CopperContent, data = datos)
summary(fit_no_interaction) Df Sum Sq Mean Sq F value Pr(>F)
Temperature 3 156.1 52.03 5.852 0.00359 **
CopperContent 3 698.3 232.78 26.181 6.98e-08 ***
Residuals 25 222.3 8.89
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La conclusión de la tabla ANOVA es que existe al menos un tipo de Temperature para el que la deformación es diferente de la media global. Además, existe también un nivel de CopperContent para el cual la deformación difiere sustancialmente de la media global. En otras palabras, que ambos factores seleccionados influencian considerablemente la cantidad de deformación
Rta b.
Can we draw conclusions for the population based on the residuals?. Explain.
Para sacar conclusiones de la población primero validaremos los supuestos a partir de los residuales. Para luego proceder a calcular los intervalos de confianza para cada factor y concluir sobre la respuesta.
Cuantificación del efecto
Ahora procedemos a determinar la magnitud del efecto de cada nivel de Temperature y CopperContent sobre la deformación. Para ello utilizamos los estimadores obtenidos al emplear el método de mínimos cuadrados ordinarios en un modelo de Regresión Lineal Múltiple:
##
fit_lm <- lm(y ~ Temperature + CopperContent, data = datos)
summary(fit_lm)
Call:
lm(formula = y ~ Temperature + CopperContent, data = datos)
Residuals:
Min 1Q Median 3Q Max
-5.4688 -1.5313 0.0313 1.9375 6.2812
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 16.469 1.395 11.809 1.01e-11 ***
Temperature75 -4.500 1.491 -3.018 0.00578 **
Temperature100 -0.875 1.491 -0.587 0.56254
Temperature125 1.500 1.491 1.006 0.32401
CopperContent60 3.375 1.491 2.264 0.03252 *
CopperContent80 5.500 1.491 3.689 0.00110 **
CopperContent100 12.750 1.491 8.552 6.82e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.982 on 25 degrees of freedom
Multiple R-squared: 0.7936, Adjusted R-squared: 0.744
F-statistic: 16.02 on 6 and 25 DF, p-value: 1.744e-07Validación de supuestos
## load function
source('https://www.dropbox.com/scl/fi/1z0gsv6g4a9eowcoqs4xf/validar_supuestos.R?rlkey=ioqrt3r3sl3vh7wa7sc82gjqh&dl=1')Esta función fue creada por el profesor Jorge Vélez, PhD
<https://jorgeivanvelez.netlify.app/> para la validación
de los supuestos básicos en RLM.
En caso de presentar algún inconveniente, por favor
repórtelo a <jvelezv@uninorte.edu.co>
Última modificación: Agosto 13, 2024
Ahora validamos el modelo:
## validación de supuestos
validar_supuestos(fit_lm) normalidad homocedasticidad independencia
"cumple" "cumple" "cumple" Intervalos de confianza
Confint(fit_lm) Estimate 2.5 % 97.5 %
(Intercept) 16.46875 13.5964793 19.341021
Temperature75 -4.50000 -7.5705865 -1.429413
Temperature100 -0.87500 -3.9455865 2.195587
Temperature125 1.50000 -1.5705865 4.570587
CopperContent60 3.37500 0.3044135 6.445587
CopperContent80 5.50000 2.4294135 8.570587
CopperContent100 12.75000 9.6794135 15.820587Conclusiones: Dado que el modelo cumple con todos los supuestos validados con una significancia estadistica del 5%. Por lo tanto, el modelo y las conclusiones que se deriven de él son válidas.
Podemos inferir de los intervalos de confianza para los estimadores de los efectos de cada nivel de temperatura, existe una temperatura cuyo efecto sobre la respuesta no contiene el \(0\) entre su intervalo de confianza, por lo cual siempre que se trabaje con una temperatura de 40 se obtendrá un incremento en la deformación. También existe una temperatura cuyo efecto sobre la deformación, está por debajo de \0 siempre en su intervalo de confianza, por lo cual siempre que se trabaje con una temperatura de 125 obtendremos un decrimiento en la cantidad de deformación.
Rta c.
If low warping is desirable, what level of copper content would you specify?
Dado que para cualquier valor de CopperContentpara cualquier valor por encima de 40% el efecto es un incremento en la cantidad de deformación, entonces se sugiere utilizar una placa con un porcentaje de cobre del 40%
Rta d.
Suppose that temperature cannot be easily controlled in the environment in which the copper plates are to be used. Does this change your answer for part c?
A través de un gráfico en el cual fijamos la temperatura y vemos el efecto que tendría unicamente el porcentaje de cobre sobre la cantidad de deformación.
## fixing temperature
with(datos,
interaction.plot(CopperContent,Temperature,y,
las = 1,
lty = 1,
lwd = 2,
col = 2:5,
main = "Warping, fixed Temperature"))
Conclusión: Para cualquier valor de temperatura, un porcentaje de 40% de cobre garantiza tener una menor cantidad de deformación.
Por lo cual si no se pudiera controlar la temperatura facilmente, se sugiere usar un porcentaje de cobre del 40%
Rta e.
Based on the results of the experiment, approximate the observed power. Use the pwr.2way of the pwr2 package in R.
resume <- summary(fit_lm)
R2 <- resume$r.squared
R2[1] 0.7935568f <- eta2_to_f(R2)
f[1] 1.960599pwr.anova.test(k=4,n=2,f=f,sig.level=0.05)
Balanced one-way analysis of variance power calculation
k = 4
n = 2
f = 1.960599
sig.level = 0.05
power = 0.8056964
NOTE: n is number in each groupLa potencia observada es de \(power =0.805\) lo cual nos indica que tenemos un 80.5% de probabilidad de encontrar el efecto que estamos buscando sobre la variable de respuesta.
Rta f.
Suppose that the experiment was conducted assuming that both factors were random. Compute the variance components.
Análisis gráfico
Ahora, procedemos a realizar un análisis gráfico de los resultados del experimento.
En el caso del factor A se tiene que:
## plot for A
p1 <- ggplot(datos, aes(x = y, y = Temperature, fill = Temperature,
colour = Temperature)) +
geom_density_ridges() +
theme_ridges() +
theme(legend.position = "none") +
ylab("Temperature") +
xlab("Warping") +
ggtitle("Distribución del Warping por Temperature")
p1
Para el factor B: CopperContent se tiene los siguiente:
## plot for B
p2 <- ggplot(datos, aes(x = y, y = CopperContent, fill = CopperContent,
colour = CopperContent)) +
geom_density_ridges() +
theme_ridges() +
theme(legend.position = "none") +
ylab("CopperContent") +
xlab("Warping") +
ggtitle("Distribución del Warping por CopperContent")
p2
Tabla ANOVA con interacción
## ANOVA including Temperature and CopperContent as random
fit <- aov(y ~ Error(Temperature*CopperContent), data = datos)
summary(fit)
Error: Temperature
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 3 156.1 52.03
Error: CopperContent
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 3 698.3 232.8
Error: Temperature:CopperContent
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 9 113.8 12.64
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 16 108.5 6.781 MSS1 <- broom::tidy((aov(y ~ Error(Temperature*CopperContent), data = datos)))[[5]][1]
MSS2 <- broom::tidy((aov(y ~ Error(Temperature*CopperContent), data = datos)))[[5]][2]
MSS3 <- broom::tidy((aov(y ~ Error(Temperature*CopperContent), data = datos)))[[5]][3]
MSS4 <- broom::tidy((aov(y ~ Error(Temperature*CopperContent), data = datos)))[[5]][4]A partir de estos resultados, los estadísticos de prueba para cada factor pueden obtenerse como:
\(F_{Temperature x CopperContent}=\) \({MS_{Temperature x CopperContent}/MSE}=\) 12.64/6.78\(=\) 1.86
\(F_{Temperature}=\) \({MS_{Temperature}/MSE}=\) 52.03/6.78\(=\) 7.67
\(F_{CopperContent}=\) \({MS_{TCopperContent}/MSE}=\) 232.78/6.78\(=\) 34.33
\(F_{Error}=\) \({MSE}=\) 6.78
Estimación de las Componentes de Varianza
A continuación se muestra la sintaxis para estimar los componentes de varianza. Para ello utilizaremos la función lmer() del paquete lme4. De los resultados, los componentes de varianza se encuentran en la columna Variance de Random Effects:
## full variance components models
vcm <- lmer(y ~ 1 + (1 | Temperature) + (1|CopperContent) + (1|Temperature:CopperContent),
data = datos)
summary(vcm)Linear mixed model fit by REML ['lmerMod']
Formula:
y ~ 1 + (1 | Temperature) + (1 | CopperContent) + (1 | Temperature:CopperContent)
Data: datos
REML criterion at convergence: 173.1
Scaled residuals:
Min 1Q Median 3Q Max
-1.74722 -0.43270 0.02776 0.65133 1.55530
Random effects:
Groups Name Variance Std.Dev.
Temperature:CopperContent (Intercept) 2.931 1.712
CopperContent (Intercept) 27.517 5.246
Temperature (Intercept) 4.924 2.219
Residual 6.781 2.604
Number of obs: 32, groups:
Temperature:CopperContent, 16; CopperContent, 4; Temperature, 4
Fixed effects:
Estimate Std. Error t value
(Intercept) 20.906 2.916 7.169Por lo tanto,
\[ \begin{eqnarray} \hat{\sigma}^2_{\text{Temperature}} &=& 27.517 \\\nonumber \hat{\sigma}^2_{\text{CopperContent}} &=& 4.924 \\\nonumber \hat{\sigma}^2_{\text{Temperature x CopperContent}} &=& 2.931\\\nonumber \hat{\sigma}^2_{\text{Error}} &=& 6.781 \end{eqnarray} \]
Finalmente,
\(\hat{\sigma}^2_y =\) 27.517 + 4.924 + 2.931 + 6.781 = 42.153
De esta variabilidad total, el factor Temperature contribuye con un 65.28%, CopperContent con un 11.68%, la interacción Temperature:CopperContent con un 6.95% y finalmente el Error un un 16.09%.
Pruebas de hipótesis
Una forma probar si las componentes de varianza asociadas a CopperContent, Temperature y CopperContent:Temperature son estadísticamente significativas para una probabilidad de Error Tipo I del 5% es a través del cálculo de intervalos de confianza. Para ello utilizaremos la función confint():
## intervalos de confianza del 95%
confint(vcm, oldNames = FALSE) 2.5 % 97.5 %
sd_(Intercept)|Temperature:CopperContent 0.000000 3.991723
sd_(Intercept)|CopperContent 2.341513 11.657819
sd_(Intercept)|Temperature 0.000000 6.333092
sigma 1.908673 3.833401
(Intercept) 14.583244 27.229261Así las cosas,
\({\sigma^2}_{Temperature x CopperContent}\in(\) 0\(^2,\) 3.99 \(^2) \rightarrow {\sigma^2}_{Temperature x CopperContent} \in (\) 0 \(,\) 15.93\()\)
\({\sigma^2}_{CopperContent}\in(\) 2.34\(^2,\) 11.66 \(^2) \rightarrow {\sigma^2}_{CopperContent} \in (\) 5.48 \(,\) 135.9 \()\)
\({\sigma^2}_{Temperature}\in(\) 0\(^2,\) 6.33 \(^2) \rightarrow {\sigma^2}_{Temperature} \in (\) 0 \(,\) 40.11 \()\)
\({\sigma^2}_{Error}\in(\) 1.91\(^2,\) 3.83 \(^2) \rightarrow {\sigma^2}_{Error} \in (\) 3.64 \(,\) 14.69 \()\)
A partir de este resultado podemos concluir que ni CopperContent ni la interacción Temperature:CopperContent son fuentes significativas de la variabilidad de y.
Validación de Supuestos
Con el fin de determinar si las conclusiones obtenidas hasta ahorap pueden o no extenderse para toda la población, es necesario [validar los supuestos.
Los residuales del modelo puede estimarse haciendo:
## cálculo residuales
r <- residuals(vcm)A partir del objeto r podemos aplicar diferentes funciones para validar los supuestos de Normalidad, Independencia y varianza constante:
## valores p para la validación de supuestos
validacion <- c('Normalidad' = shapiro.test(r)$p.value,
'Independencia (CopperContent)' = dwtest(r ~ CopperContent, data = datos)$p.value,
'Independencia (Temperature)' = dwtest(r ~ Temperature, data = datos)$p.value,
'Varianza Constante (CopperContent)' = bptest(r ~ CopperContent, data = datos)$p.value,
'Varianza Constante (Temperature)' = bptest(r ~ Temperature, data = datos)$p.value)
validacion <- round(validacion, 3)
as.list(validacion)$Normalidad
[1] 0.419
$`Independencia (CopperContent)`
[1] 0.953
$`Independencia (Temperature)`
[1] 0.941
$`Varianza Constante (CopperContent).BP`
[1] 0.927
$`Varianza Constante (Temperature).BP`
[1] 0.154Conclusión: Como los valores \(p\) de todas las pruebas son \(>0.05\), podemos afirmar que los resultados obtenidos son generalizables para la población de la que proviene la muestra.
Rta. g
Use linear regression to estimate the warping at 80° and 90% copper content.
##Volver a leer datos
datos <- read.table('HomeWork2_files/5.12 Ejercicio.txt',
header = TRUE)
fit_lm <- lm(y ~ Temperature + CopperContent, data = datos)
## convert 'Temperature' and 'CopperContent' to factors
broom::tidy(summary(fit_lm))# A tibble: 3 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 3.93 2.97 1.32 0.196
2 Temperature 0.0325 0.0234 1.39 0.176
3 CopperContent 0.202 0.0293 6.89 0.000000143b1 <-broom::tidy(summary(fit_lm))[[2]][1]
b2 <-broom::tidy(summary(fit_lm))[[2]][2]
b3 <-broom::tidy(summary(fit_lm))[[2]][3]
head(datos,10) Temperature CopperContent y
1 50 40 17
2 50 40 20
3 50 60 16
4 50 60 21
5 50 80 24
6 50 80 22
7 50 100 28
8 50 100 27
9 75 40 12
10 75 40 9summary(fit_lm)
Call:
lm(formula = y ~ Temperature + CopperContent, data = datos)
Residuals:
Min 1Q Median 3Q Max
-10.5188 -1.7875 0.7813 2.3781 6.3687
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.93125 2.97332 1.322 0.196
Temperature 0.03250 0.02344 1.387 0.176
CopperContent 0.20188 0.02930 6.891 1.43e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.706 on 29 degrees of freedom
Multiple R-squared: 0.6301, Adjusted R-squared: 0.6046
F-statistic: 24.7 on 2 and 29 DF, p-value: 5.456e-07coefficients <- coef(fit_lm)
print(coefficients) (Intercept) Temperature CopperContent
3.931250 0.032500 0.201875 interce <- b1
beta_temp <- b2
beta_copper <- b3
temperature <- 80
copper_content_90 <- 90
prediction_90 <- interce + beta_temp * temperature + beta_copper * copper_content_90 + 3.706
print(paste("Predicción para 80° y 90% Copper Content:", prediction_90))[1] "Predicción para 80° y 90% Copper Content: 28.406"Superficie de Respuesta
A partir del modelo de RLM ajustado, contenido en el objeto fitlm, es posible utilizar grid-search para construir la superficie de respuesta asociada y determinar las condiciones de operación óptimas que garantizan un valor minimo de Deformación.
Inicialmente creamos la rejilla de evaluación:
# grid-based optimization
ngrid <- 30
yLm <- lm(y ~ Temperature + CopperContent, data = datos)
Agrid <- Bgrid <- seq(from = 60, to = 100, length = ngrid)
yhat <- predict(yLm, expand.grid(Temperature = Agrid, CopperContent = Bgrid))
yhat <- matrix(yhat, length(Agrid), length(Bgrid))y posteriormente procedemos a graficar la superficie y los contornos:
## disponibilidad del paquete plot3D
if (!require(plot3D)) install.packages("plot3D")
require(plot3D)
## rearranging the window
par(mfrow = c(1, 2))
## 3D plot (same as Fig 6.3a)
persp3D( Agrid, Bgrid, yhat,
ticktype = "detailed",
colkey = FALSE, phi = 30, theta = -40, bty = "b2",
space = 0.1, d = 1.5,
xlab = "\n Temperature",
ylab = "\n CopperContent",
zlab = "\n\nPredicted Deformación", contour = TRUE,
zlim = c(18, 30))
##Agregar planos en Gráfica 3D
x_plane <- c(80, 80, 80, 80)
y_plane <- c(min(Bgrid), max(Bgrid), max(Bgrid), min(Bgrid))
z_plane <- c(min(yhat), min(yhat), max(yhat), max(yhat))
polygon3D(x = x_plane, y = y_plane, z = z_plane,
col = "blue", add = TRUE)
x_plane_2 <- c(min(Agrid), max(Agrid), max(Agrid), min(Agrid))
y_plane_2 <- c(90, 90, 90, 90)
z_plane_2 <- c(min(yhat), min(yhat), max(yhat), max(yhat))
polygon3D(x = x_plane_2, y = y_plane_2, z = z_plane_2,
col = "red", add = TRUE)
## 2D plot (same as Fig 6.3b)
contour2D(yhat, Agrid, Bgrid,
xlab = "Temperature (°)",
ylab = "CopperContent (%)", colkey = FALSE, las = 1)
# Agregar líneas a la gráfica
# Por ejemplo, agregar una línea horizontal en y = 80
abline(h = 90, col = "red", lwd = 2)
# Agregar una línea vertical en x = 75
abline(v = 80, col = "blue", lwd = 2)
# Definir los valores para los cuales se desea hacer la predicción
new_data <- data.frame(Temperature = 80, CopperContent = 90)
# Predecir el valor de 'y' para los valores especificados
predicted_value <- predict(yLm, newdata = new_data)
# Mostrar el valor predicho
predicted_value 1
24.7 El valor esperado de Warping cuando la temperatura es y El porcentaje de cobre 80, 90 es 24.7