library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.5.2
## āā Attaching core tidyverse packages āāāāāāāāāāāāāāāāāāāāāāāā tidyverse 2.0.0 āā
## ā dplyr 1.1.4 ā readr 2.1.5
## ā forcats 1.0.0 ā stringr 1.5.1
## ā ggplot2 4.0.1 ā tibble 3.3.0
## ā lubridate 1.9.4 ā tidyr 1.3.1
## ā purrr 1.1.0
## āā Conflicts āāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāāā tidyverse_conflicts() āā
## ā dplyr::filter() masks stats::filter()
## ā dplyr::lag() masks stats::lag()
## ā¹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(car)
## Cargando paquete requerido: carData
##
## Adjuntando el paquete: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
library(lmtest)
## Cargando paquete requerido: zoo
##
## Adjuntando el paquete: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(rsm)
## Warning: package 'rsm' was built under R version 4.5.2
library(ggpubr)
library(plotly)
##
## Adjuntando el paquete: 'plotly'
##
## 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':
##
## layout
library(effectsize)
## Warning: package 'effectsize' was built under R version 4.5.2
library(patchwork)
library(MASS)
##
## Adjuntando el paquete: 'MASS'
##
## The following object is masked from 'package:patchwork':
##
## area
##
## The following object is masked from 'package:plotly':
##
## select
##
## The following object is masked from 'package:dplyr':
##
## select
dat <- read.csv("galletas_BB_sim.csv")
glimpse(dat)
## Rows: 15
## Columns: 9
## $ run_order <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15
## $ Temp <int> 180, 165, 180, 180, 165, 165, 180, 195, 165, 180, 195, 180, ā¦
## $ Time <int> 10, 12, 10, 8, 8, 10, 8, 10, 10, 12, 10, 10, 8, 12, 12
## $ Butter <int> 15, 15, 15, 20, 15, 10, 10, 20, 20, 10, 10, 15, 15, 20, 15
## $ x1 <int> 0, -1, 0, 0, -1, -1, 0, 1, -1, 0, 1, 0, 1, 0, 1
## $ x2 <int> 0, 1, 0, -1, -1, 0, -1, 0, 0, 1, 0, 0, -1, 1, 1
## $ x3 <int> 0, 0, 0, 1, 0, -1, -1, 1, 1, -1, -1, 0, 0, 1, 0
## $ Firmness <dbl> 34.27917, 23.36742, 31.00729, 46.14353, 43.40646, 43.59924, ā¦
## $ Lstar <dbl> 58.22189, 68.39864, 61.09426, 56.53890, 59.93838, 67.86715, ā¦
dat <- dat %>%
mutate(
x1 = (Temp - 180)/15,
x2 = (Time - 10)/2,
x3 = (Butter - 15)/5
)
table(dat$x1, dat$x2, dat$x3)
## , , = -1
##
##
## -1 0 1
## -1 0 1 0
## 0 1 0 1
## 1 0 1 0
##
## , , = 0
##
##
## -1 0 1
## -1 1 0 1
## 0 0 3 0
## 1 1 0 1
##
## , , = 1
##
##
## -1 0 1
## -1 0 1 0
## 0 1 0 1
## 1 0 1 0
count_centers <- dat %>%
filter(x1==0 & x2==0 & x3==0) %>%
nrow()
count_centers
## [1] 3
mod_firm_q <- lm(Firmness ~ x1 + x2 + x3 +
I(x1^2) + I(x2^2) + I(x3^2) +
x1:x2 + x1:x3 + x2:x3,
data = dat)
summary(mod_firm_q)
##
## Call:
## lm(formula = Firmness ~ x1 + x2 + x3 + I(x1^2) + I(x2^2) + I(x3^2) +
## x1:x2 + x1:x3 + x2:x3, data = dat)
##
## Residuals:
## 1 2 3 4 5 6 7 8 9 10
## 2.5749 -0.4778 -0.6969 1.1182 -0.6479 1.5960 -0.9481 -1.5960 -0.4703 -1.1182
## 11 12 13 14 15
## 0.4703 -1.8780 0.4778 0.9481 0.6479
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.7042 1.2048 26.314 1.48e-06 ***
## x1 5.3727 0.7378 7.282 0.000764 ***
## x2 -6.1501 0.7378 -8.336 0.000406 ***
## x3 -3.2631 0.7378 -4.423 0.006874 **
## I(x1^2) 4.7402 1.0860 4.365 0.007258 **
## I(x2^2) 2.8780 1.0860 2.650 0.045422 *
## I(x3^2) 5.0198 1.0860 4.622 0.005724 **
## x1:x2 3.9545 1.0434 3.790 0.012759 *
## x1:x3 2.6486 1.0434 2.538 0.051992 .
## x2:x3 -2.5363 1.0434 -2.431 0.059326 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.087 on 5 degrees of freedom
## Multiple R-squared: 0.9768, Adjusted R-squared: 0.935
## F-statistic: 23.38 on 9 and 5 DF, p-value: 0.001442
anova(mod_firm_q)
mod_L_q <- lm(Lstar ~ x1 + x2 + x3 +
I(x1^2) + I(x2^2) + I(x3^2) +
x1:x2 + x1:x3 + x2:x3,
data = dat)
summary(mod_L_q)
##
## Call:
## lm(formula = Lstar ~ x1 + x2 + x3 + I(x1^2) + I(x2^2) + I(x3^2) +
## x1:x2 + x1:x3 + x2:x3, data = dat)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -1.93875 -0.02217 0.93362 0.90817 -0.28897 0.93035 -0.64138 -0.93035
## 9 10 11 12 13 14 15
## -0.61920 -0.90817 0.61920 1.00513 0.02217 0.64138 0.28897
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 60.16063 0.84770 70.969 1.05e-08 ***
## x1 -1.58865 0.51911 -3.060 0.02809 *
## x2 2.85241 0.51911 5.495 0.00273 **
## x3 -2.17683 0.51911 -4.193 0.00854 **
## I(x1^2) 1.72847 0.76410 2.262 0.07316 .
## I(x2^2) 0.84632 0.76410 1.108 0.31845
## I(x3^2) 1.18479 0.76410 1.551 0.18170
## x1:x2 -1.24432 0.73413 -1.695 0.15085
## x1:x3 0.09743 0.73413 0.133 0.89960
## x2:x3 1.53179 0.73413 2.087 0.09130 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.468 on 5 degrees of freedom
## Multiple R-squared: 0.9352, Adjusted R-squared: 0.8185
## F-statistic: 8.013 on 9 and 5 DF, p-value: 0.0169
anova(mod_L_q)
par(mfrow=c(2,2))
plot(mod_firm_q)
par(mfrow=c(1,1))
shapiro.test(residuals(mod_firm_q))
##
## Shapiro-Wilk normality test
##
## data: residuals(mod_firm_q)
## W = 0.96551, p-value = 0.7869
bptest(mod_firm_q)
##
## studentized Breusch-Pagan test
##
## data: mod_firm_q
## BP = 8.4602, df = 9, p-value = 0.4885
infl <- influence.measures(mod_firm_q)
summary(infl)
## Potentially influential observations of
## lm(formula = Firmness ~ x1 + x2 + x3 + I(x1^2) + I(x2^2) + I(x3^2) + x1:x2 + x1:x3 + x2:x3, data = dat) :
##
## dfb.1_ dfb.x1 dfb.x2 dfb.x3 dfb.I(1^ dfb.I(2^ dfb.I(3^ dfb.x1:2 dfb.x1:3
## 1 1.30_* 0.00 0.00 0.00 -0.72 -0.72 -0.72 0.00 0.00
## 2 0.00 0.30 -0.30 0.00 -0.20 -0.20 0.20 0.42 0.00
## 3 -0.26 0.00 0.00 0.00 0.15 0.15 0.15 0.00 0.00
## 4 0.00 0.00 -0.77 0.77 -0.52 0.52 0.52 0.00 0.00
## 5 0.00 0.41 0.41 0.00 -0.28 -0.28 0.28 -0.58 0.00
## 6 0.00 -1.33_* 0.00 -1.33_* 0.90 -0.90 0.90 0.00 1.88_*
## 8 0.00 -1.33_* 0.00 -1.33_* -0.90 0.90 -0.90 0.00 -1.88_*
## 9 0.00 0.29 0.00 -0.29 -0.20 0.20 -0.20 0.00 0.41
## 10 0.00 0.00 -0.77 0.77 0.52 -0.52 -0.52 0.00 0.00
## 11 0.00 0.29 0.00 -0.29 0.20 -0.20 0.20 0.00 -0.41
## 13 0.00 0.30 -0.30 0.00 0.20 0.20 -0.20 -0.42 0.00
## 15 0.00 0.41 0.41 0.00 0.28 0.28 -0.28 0.58 0.00
## dfb.x2:3 dffit cov.r cook.d hat
## 1 0.00 1.30 0.03 0.11 0.33
## 2 0.00 -0.72 24.27_* 0.06 0.75
## 3 0.00 -0.26 9.94_* 0.01 0.33
## 4 -1.09_* 1.89 2.74 0.34 0.75
## 5 0.00 -1.00 16.70_* 0.12 0.75
## 6 0.00 3.25 0.07 0.70 0.75
## 8 0.00 -3.25 0.07 0.70 0.75
## 9 0.00 -0.71 24.60_* 0.06 0.75
## 10 1.09_* -1.89 2.74 0.34 0.75
## 11 0.00 0.71 24.60_* 0.06 0.75
## 13 0.00 0.72 24.27_* 0.06 0.75
## 15 0.00 1.00 16.70_* 0.12 0.75
cook <- cooks.distance(mod_firm_q)
which(cook > 4/nrow(dat))
## 4 6 8 10
## 4 6 8 10
bc_firm <- boxcox(mod_firm_q, plotit=TRUE)
lambda_opt <- bc_firm$x[which.max(bc_firm$y)]
lambda_opt
## [1] 1.515152
grid <- expand.grid(
x1 = seq(-1,1,length=60),
x2 = seq(-1,1,length=60),
x3 = 0
)
grid$pred_firm <- predict(mod_firm_q, newdata=grid)
grid <- grid %>%
mutate(
temp = 180 + 15*x1,
tiempo = 10 + 2*x2
)
ggplot(grid, aes(x=temp, y=tiempo, z=pred_firm)) +
geom_raster(aes(fill=pred_firm)) +
geom_contour(col="white", bins=12) +
labs(
title="Superficie predicha (Firmness) ā mantequilla = 15%",
x="Temp (Ā°C)", y="Tiempo (min)"
) +
theme_minimal()
plot_ly(
x = grid$temp,
y = grid$tiempo,
z = grid$pred_firm,
type = "surface"
) %>% layout(
scene = list(
xaxis=list(title="Temp (Ā°C)"),
yaxis=list(title="Tiempo (min)"),
zaxis=list(title="Firmness (N)")
)
)
cat("\n===== RESUMEN INTERPRETATIVO (Firmness) =====\n")
##
## ===== RESUMEN INTERPRETATIVO (Firmness) =====
summary(mod_firm_q)
##
## Call:
## lm(formula = Firmness ~ x1 + x2 + x3 + I(x1^2) + I(x2^2) + I(x3^2) +
## x1:x2 + x1:x3 + x2:x3, data = dat)
##
## Residuals:
## 1 2 3 4 5 6 7 8 9 10
## 2.5749 -0.4778 -0.6969 1.1182 -0.6479 1.5960 -0.9481 -1.5960 -0.4703 -1.1182
## 11 12 13 14 15
## 0.4703 -1.8780 0.4778 0.9481 0.6479
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.7042 1.2048 26.314 1.48e-06 ***
## x1 5.3727 0.7378 7.282 0.000764 ***
## x2 -6.1501 0.7378 -8.336 0.000406 ***
## x3 -3.2631 0.7378 -4.423 0.006874 **
## I(x1^2) 4.7402 1.0860 4.365 0.007258 **
## I(x2^2) 2.8780 1.0860 2.650 0.045422 *
## I(x3^2) 5.0198 1.0860 4.622 0.005724 **
## x1:x2 3.9545 1.0434 3.790 0.012759 *
## x1:x3 2.6486 1.0434 2.538 0.051992 .
## x2:x3 -2.5363 1.0434 -2.431 0.059326 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.087 on 5 degrees of freedom
## Multiple R-squared: 0.9768, Adjusted R-squared: 0.935
## F-statistic: 23.38 on 9 and 5 DF, p-value: 0.001442
cat("\nANOVA:\n")
##
## ANOVA:
anova(mod_firm_q)
pred_both <- function(xx){
newd <- data.frame(x1=xx[1], x2=xx[2], x3=xx[3])
f <- predict(mod_firm_q, newdata=newd)
L <- predict(mod_L_q, newdata=newd)
return(c(f=f, L=L))
}