Factorial design

Si vuole misurare una certa risposta sperimentale al variare di tre parametri specifici (due di processo e uno qualitativo):

• temperatura
• durata del trattamento
• tipo di legno

Si inizia creando una tabella dove ogni colonna è un fattore e ogni riga è un livello secondo il modello postulato ottenuto con \(2^3\) condizioni sperimentali, ovvero 3 fattori e 2 livelli, prevedendo almeno un paio di repliche.

---

df = data.frame(
  Temperature = c(160,180,160,180,160,180,160,180,160,180,160,180,160,180,160,180),
  Duration = c(20,20,40,40,20,20,40,40,20,20,40,40,20,20,40,40),
  Wood_type = c('A','A','A','A','B','B','B','B','A','A','A','A','B','B','B','B'),
  Y = c(59,74,50,69,50,81,46,79,61,70,58,67,54,85,44,81)
)
#DT::datatable(df)
df %>% kable() %>% kable_material('striped')

Temperature	Duration	Wood_type	Y
160	20	A	59
180	20	A	74
160	40	A	50
180	40	A	69
160	20	B	50
180	20	B	81
160	40	B	46
180	40	B	79
160	20	A	61
180	20	A	70
160	40	A	58
180	40	A	67
160	20	B	54
180	20	B	85
160	40	B	44
180	40	B	81

plot_ly(df, x = ~Temperature, y = ~Duration, z = ~Wood_type) %>%
  add_markers()  

---

Dopodiché creo una tabella con i valori precedenti codificandoli come \(-1\) e \(+1\)

x1 = c()
x2 = c()
x3 = c()

levels = c(-1,1)

for (k in 1:2){
  for (j in 1:2) {
    for (i in 1:2) {
      x1 = c(x1, levels[i])
      x2 = c(x2, levels[j])
      x3 = c(x3, levels[k])
    }
  }
}

x1 = c(x1,x1)
x2 = c(x2,x2)
x3 = c(x3,x3)

df1 = data.frame(x1,x2,x3,Y = c(59,74,50,69,50,81,46,79,61,70,58,67,54,85,44,81))

fig = df1 %>% kable() %>% kable_styling( full_width = T) 

for (m in 1:(ncol(df1)-1)) {
 fig = column_spec(fig, m ,color= 'white', background = spec_color(df1[,m],end = 0.7, option = 'A' ))
}

fig %>%  column_spec( ncol(df1) ,color= 'white', background = 'darkgreen')

x1	x2	x3	Y
-1	-1	-1	59
1	-1	-1	74
-1	1	-1	50
1	1	-1	69
-1	-1	1	50
1	-1	1	81
-1	1	1	46
1	1	1	79
-1	-1	-1	61
1	-1	-1	70
-1	1	-1	58
1	1	-1	67
-1	-1	1	54
1	-1	1	85
-1	1	1	44
1	1	1	81

---

Adesso aggiungo le altre colonne seguendo il modello postulato

\[ Y = b_0+b_1x_1+b_2x_2+b_3x_3+b_{12}x_1x_2+b_{13}x_1x_3+b_{23}x_2x_3+b_{123}x_1x_2x_3 \]

df3 = mutate(df1, x1x2 = x1*x2, x1x3= x1*x3, x2x3=x2*x3, x1x2x3 = x1*x2*x3)
df3 = relocate(df3, Y, .after = last_col())

fig = df3 %>% kable() %>% kable_styling( full_width = T) 

for (m in 1:(ncol(df3)-1)) {
 fig = column_spec(fig, m ,color= 'white', background = spec_color(df3[,m],end = 0.7, option = 'A' ))
}

fig %>%  column_spec( ncol(df3) ,color= 'white', background = 'darkgreen')

x1	x2	x3	x1x2	x1x3	x2x3	x1x2x3	Y
-1	-1	-1	1	1	1	-1	59
1	-1	-1	-1	-1	1	1	74
-1	1	-1	-1	1	-1	1	50
1	1	-1	1	-1	-1	-1	69
-1	-1	1	1	-1	-1	1	50
1	-1	1	-1	1	-1	-1	81
-1	1	1	-1	-1	1	-1	46
1	1	1	1	1	1	1	79
-1	-1	-1	1	1	1	-1	61
1	-1	-1	-1	-1	1	1	70
-1	1	-1	-1	1	-1	1	58
1	1	-1	1	-1	-1	-1	67
-1	-1	1	1	-1	-1	1	54
1	-1	1	-1	1	-1	-1	85
-1	1	1	-1	-1	1	-1	44
1	1	1	1	1	1	1	81

---

Creo un modello di regressione

a = c()
for (i in colnames(df3[-ncol(df3)])){
  a = paste(a , '+', i, sep = '')
}
fmla = as.formula(paste('Y ~ ', a))

model = lm(fmla, df3)
model

## 
## Call:
## lm(formula = fmla, data = df3)
## 
## Coefficients:
## (Intercept)           x1           x2           x3         x1x2         x1x3  
##   6.425e+01    1.150e+01   -2.500e+00    7.500e-01    7.500e-01    5.000e+00  
##        x2x3       x1x2x3  
##   7.216e-16    2.500e-01

## 
## Call:
## lm(formula = fmla, data = df3)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -4.00  -1.25   0.00   1.25   4.00 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.425e+01  7.071e-01  90.863 2.40e-13 ***
## x1           1.150e+01  7.071e-01  16.263 2.06e-07 ***
## x2          -2.500e+00  7.071e-01  -3.536 0.007670 ** 
## x3           7.500e-01  7.071e-01   1.061 0.319813    
## x1x2         7.500e-01  7.071e-01   1.061 0.319813    
## x1x3         5.000e+00  7.071e-01   7.071 0.000105 ***
## x2x3         7.216e-16  7.071e-01   0.000 1.000000    
## x1x2x3       2.500e-01  7.071e-01   0.354 0.732810    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.828 on 8 degrees of freedom
## Multiple R-squared:  0.9763, Adjusted R-squared:  0.9555 
## F-statistic: 47.05 on 7 and 8 DF,  p-value: 7.071e-06

---

Osservo la significatività di P e scarto i valori superiori al mio limite di confidenza (5%). Il P-value può essere visto come la probabilità di sbagliare nel rigettare l’ipotesi nulla.

(Noto che c’è una correlazione tra \(x_1\) e \(x_3\) (temperatura e tipo di legno), quindi tengo come significativo anche il vettore \(x_3\))

---

Ripeto la regressione dopo aver scartatato dalla mia tabella le tre colonne non significative: \(x_1x_2,x_2x_3,x_1x_2x_3\)

model2 = lm(Y ~ x1 + x2 + x3 + x1x3, df3)
summary(model2)

## 
## Call:
## lm(formula = Y ~ x1 + x2 + x3 + x1x3, data = df3)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -4.50  -1.25   0.00   1.50   3.50 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  64.2500     0.6484  99.086  < 2e-16 ***
## x1           11.5000     0.6484  17.735 1.93e-09 ***
## x2           -2.5000     0.6484  -3.855  0.00267 ** 
## x3            0.7500     0.6484   1.157  0.27192    
## x1x3          5.0000     0.6484   7.711 9.25e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.594 on 11 degrees of freedom
## Multiple R-squared:  0.9726, Adjusted R-squared:  0.9626 
## F-statistic: 97.55 on 4 and 11 DF,  p-value: 1.629e-08

---

Il vettore con i valori calcolati sarà

fitted(model2)

##    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
## 59.5 72.5 54.5 67.5 51.0 84.0 46.0 79.0 59.5 72.5 54.5 67.5 51.0 84.0 46.0 79.0

---

Posso a questo punto crare un grafico per vedere i valori osservati sperimentalmente e quelli calcolati dal modello di regressione

---

df_y = data.frame(Y_cal=fitted(model2), Y_sper = df3$Y)

bisettrice = data.frame(x=c(40,90), y=c(40,90))

ggplot(df_y, aes(Y_sper,Y_cal))+ 
  geom_point()+ 
  geom_line(data=bisettrice, aes(x,y, color='red'), show.legend = FALSE)+ 
  ggtitle( 'Y sperimentale vs Y calcolata')+ 
  theme_minimal()+
  theme(plot.title = element_text(hjust=0.5))

modelcoeff = coefficients(model2)
modelcoeff %>% round(digits = 2)

## (Intercept)          x1          x2          x3        x1x3 
##       64.25       11.50       -2.50        0.75        5.00

x = seq(from = -1,  to = 1, by = 0.1)
y = seq(from = -1,  to = 1, by = 0.1)
x3 = c(-1,1)
z = c()

for (k in x3) {
  for (j in y){
    for (i in x){
    z = c(z, modelcoeff[1]+modelcoeff[2]*i+modelcoeff[3]*j+ modelcoeff[4]*k+modelcoeff[5]*i*k)
    }
  }
}
z = array(z, dim = c(length(x), length(y),length(x3)))

Z = melt(z)

Temp = seq(from = 160, to = 180, by = 1) 
Time = seq(from = 20, to = 40, by = 1)
Wood = c(1:2)
Wd = c()
Tm = c()
Tp = c()

for (k in Wood) {
  for (j in Time){
    for (i in Temp){
      Tp = c(Tp, i)
      Tm = c(Tm, j)
      Wd = c(Wd, k)
    
    }
  }
}


plot_ly(x=Tp, y= Tm, z= Wd, color = Z$value) %>% add_markers() %>%
  layout(scene = list(
                 xaxis = list(title = 'Temperature'),
                 yaxis = list(title = 'Time'),
                 zaxis = list(title = 'Wood_type'),
                 aspectmode = 'cube')) %>%
  colorbar(title = "Yield")