Thickeners for gelato mixtures

DOE and linear regression model

Setting up the experimental design

We wanted to evaluate the effects of three thickeners of a generic ice cream mix on the viscosity of the mix and on the overrun of the finished ice cream, two parameters that influences in an important way the characteristics of the finished product.

We examined three thickeners commonly used as food additives:

• E410 - Locust bean gum or carob seed flour (FSC)
• E417 - Tara gum (TARA)
• E412 - Guar gum (GUAR)

We set up a mixture design of the Simplex Lattice type with k=3 and m=3, whose possible combinations led to carrying out 10 experimental conditions. We performed 3 replicates in 3 different weeks, although for this analysis we chose to consider only the last two sets of replicates..

---

df = read.csv("AddensantiDOE.csv", sep = ";")
week = rep(c(1,2,3), each=10 ) %>% factor()
df = add_column(df, Week = week, .after = "ID")


kable(df, row.names = T) %>% kable_styling('striped', fixed_thead = T, full_width = FALSE) %>%
  scroll_box( height = "500px")

	ID	Week	FSC	TARA	GUAR	Viscosity_cP	Overrun
1	4_0_8	1	0.3333	0.0000	0.6667	4760	0.4845
2	4_4_4	1	0.3333	0.3333	0.3333	8750	0.3459
3	0_12_0	1	0.0000	1.0000	0.0000	9900	0.3039
4	0_8_4	1	0.0000	0.6667	0.3333	9800	0.2744
5	12_0_0	1	1.0000	0.0000	0.0000	6500	0.3690
6	8_0_4	1	0.6667	0.0000	0.3333	7400	0.3447
7	4_8_0	1	0.3333	0.6667	0.0000	9910	0.3163
8	0_0_12	1	0.0000	0.0000	1.0000	5480	0.5092
9	0_4_8	1	0.0000	0.3333	0.6667	5750	0.4532
10	8_4_0	1	0.6667	0.3333	0.0000	8990	0.3372
11	8_0_4	2	0.6667	0.0000	0.3333	7050	0.3572
12	0_8_4	2	0.0000	0.6667	0.3333	7200	0.4339
13	0_4_8	2	0.0000	0.3333	0.6667	6810	0.5000
14	4_8_0	2	0.3333	0.6667	0.0000	8370	0.3800
15	4_4_4	2	0.3333	0.3333	0.3333	7080	0.3968
16	0_0_12	2	0.0000	0.0000	1.0000	5790	0.5429
17	4_0_8	2	0.3333	0.0000	0.6667	6330	0.4619
18	0_12_0	2	0.0000	1.0000	0.0000	9330	0.2734
19	8_4_0	2	0.6667	0.3333	0.0000	7950	0.3647
20	12_0_0	2	1.0000	0.0000	0.0000	6330	0.3572
21	0_0_12	3	0.0000	0.0000	1.0000	4272	0.5653
22	8_4_0	3	0.6667	0.3333	0.0000	6780	0.3800
23	0_8_4	3	0.0000	0.6667	0.3333	5610	0.4935
24	4_8_0	3	0.3333	0.6667	0.0000	7500	0.3911
25	4_4_4	3	0.3333	0.3333	0.3333	6510	0.3968
26	12_0_0	3	1.0000	0.0000	0.0000	5730	0.3403
27	0_12_0	3	0.0000	1.0000	0.0000	8130	0.3498
28	4_0_8	3	0.3333	0.0000	0.6667	6450	0.5185
29	0_4_8	3	0.0000	0.3333	0.6667	5610	0.4496
30	8_0_4	3	0.6667	0.0000	0.3333	5610	0.4082

Postuleted model

Once the experiments have been carried out the data can be processed according to a simplified Scheffè model (special cubic), which allows to take into account the effects given by the possible interactions between the 3 variables FSC, TARA and GUAR.

The model is as follows:

\[ 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 \]

Additional columns related to 2-components interaction and 3-component interaction are then added to the starting dataset.

	ID	Week	FSC	TARA	GUAR	F_T	T_G	F_G	F_T_G	Viscosity_cP	Overrun
1	4_0_8	1	0.3333	0.0000	0.6667	0.0000	0.0000	0.2222	0.000	4760	0.4845
2	4_4_4	1	0.3333	0.3333	0.3333	0.1111	0.1111	0.1111	0.037	8750	0.3459
3	0_12_0	1	0.0000	1.0000	0.0000	0.0000	0.0000	0.0000	0.000	9900	0.3039
4	0_8_4	1	0.0000	0.6667	0.3333	0.0000	0.2222	0.0000	0.000	9800	0.2744
5	12_0_0	1	1.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.000	6500	0.3690
6	8_0_4	1	0.6667	0.0000	0.3333	0.0000	0.0000	0.2222	0.000	7400	0.3447
7	4_8_0	1	0.3333	0.6667	0.0000	0.2222	0.0000	0.0000	0.000	9910	0.3163
8	0_0_12	1	0.0000	0.0000	1.0000	0.0000	0.0000	0.0000	0.000	5480	0.5092
9	0_4_8	1	0.0000	0.3333	0.6667	0.0000	0.2222	0.0000	0.000	5750	0.4532
10	8_4_0	1	0.6667	0.3333	0.0000	0.2222	0.0000	0.0000	0.000	8990	0.3372
11	8_0_4	2	0.6667	0.0000	0.3333	0.0000	0.0000	0.2222	0.000	7050	0.3572
12	0_8_4	2	0.0000	0.6667	0.3333	0.0000	0.2222	0.0000	0.000	7200	0.4339
13	0_4_8	2	0.0000	0.3333	0.6667	0.0000	0.2222	0.0000	0.000	6810	0.5000
14	4_8_0	2	0.3333	0.6667	0.0000	0.2222	0.0000	0.0000	0.000	8370	0.3800
15	4_4_4	2	0.3333	0.3333	0.3333	0.1111	0.1111	0.1111	0.037	7080	0.3968
16	0_0_12	2	0.0000	0.0000	1.0000	0.0000	0.0000	0.0000	0.000	5790	0.5429
17	4_0_8	2	0.3333	0.0000	0.6667	0.0000	0.0000	0.2222	0.000	6330	0.4619
18	0_12_0	2	0.0000	1.0000	0.0000	0.0000	0.0000	0.0000	0.000	9330	0.2734
19	8_4_0	2	0.6667	0.3333	0.0000	0.2222	0.0000	0.0000	0.000	7950	0.3647
20	12_0_0	2	1.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.000	6330	0.3572
21	0_0_12	3	0.0000	0.0000	1.0000	0.0000	0.0000	0.0000	0.000	4272	0.5653
22	8_4_0	3	0.6667	0.3333	0.0000	0.2222	0.0000	0.0000	0.000	6780	0.3800
23	0_8_4	3	0.0000	0.6667	0.3333	0.0000	0.2222	0.0000	0.000	5610	0.4935
24	4_8_0	3	0.3333	0.6667	0.0000	0.2222	0.0000	0.0000	0.000	7500	0.3911
25	4_4_4	3	0.3333	0.3333	0.3333	0.1111	0.1111	0.1111	0.037	6510	0.3968
26	12_0_0	3	1.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.000	5730	0.3403
27	0_12_0	3	0.0000	1.0000	0.0000	0.0000	0.0000	0.0000	0.000	8130	0.3498
28	4_0_8	3	0.3333	0.0000	0.6667	0.0000	0.0000	0.2222	0.000	6450	0.5185
29	0_4_8	3	0.0000	0.3333	0.6667	0.0000	0.2222	0.0000	0.000	5610	0.4496
30	8_0_4	3	0.6667	0.0000	0.3333	0.0000	0.0000	0.2222	0.000	5610	0.4082

---

Viscosity

Regression

A regression is made to find the optimal coefficients according to the model, evaluating their significance (95% confidence limit).

Being a model relative to a ternary mixture I don’t have the value relative to the intercept as for a traditional linear regression.

---

model = lm(Viscosity_cP~ 0+FSC+TARA+GUAR+F_T+F_G+T_G+F_T_G, df_complete)
summary(model)

## 
## Call:
## lm(formula = Viscosity_cP ~ 0 + FSC + TARA + GUAR + F_T + F_G + 
##     T_G + F_T_G, data = df_complete)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1856.37  -752.17    98.34   657.45  2333.63 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## FSC     6300.4      575.3  10.952 1.34e-10 ***
## TARA    9102.2      575.3  15.822 7.44e-14 ***
## GUAR    5084.8      575.3   8.839 7.44e-09 ***
## F_T     2469.3     2689.2   0.918    0.368    
## F_G     2583.5     2689.2   0.961    0.347    
## T_G    -1335.7     2689.2  -0.497    0.624    
## F_T_G   5544.8    19725.4   0.281    0.781    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1048 on 23 degrees of freedom
## Multiple R-squared:  0.9838, Adjusted R-squared:  0.9789 
## F-statistic: 200.1 on 7 and 23 DF,  p-value: < 2.2e-16

It is clear that all the coefficients related to the 4 interactions are statistically not significant therefore I remove the variable with the greater P-value and I repeat the procedure (this because the three starting vectors are linearly dependent in this kind of model).

---

model = lm(Viscosity_cP~ 0+FSC+TARA+GUAR+F_T+F_G+T_G, df_complete)
summary(model)

## 
## Call:
## lm(formula = Viscosity_cP ~ 0 + FSC + TARA + GUAR + F_T + F_G + 
##     T_G, data = df_complete)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1891.56  -730.02    92.47   622.26  2298.44 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## FSC    6276.9      558.2  11.245 4.74e-11 ***
## TARA   9078.7      558.2  16.265 1.84e-14 ***
## GUAR   5061.3      558.2   9.068 3.21e-09 ***
## F_T    2733.3     2471.1   1.106    0.280    
## F_G    2847.5     2471.1   1.152    0.261    
## T_G   -1071.8     2471.1  -0.434    0.668    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1027 on 24 degrees of freedom
## Multiple R-squared:  0.9838, Adjusted R-squared:  0.9797 
## F-statistic: 242.7 on 6 and 24 DF,  p-value: < 2.2e-16

Also F_T is removed

model = lm(Viscosity_cP~ 0+FSC+TARA+GUAR+F_G+T_G, df_complete)
summary(model)

## 
## Call:
## lm(formula = Viscosity_cP ~ 0 + FSC + TARA + GUAR + F_G + T_G, 
##     data = df_complete)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1992.78  -699.19   -65.49   499.26  2197.22 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## FSC    6580.6      488.1  13.481 5.70e-13 ***
## TARA   9382.4      488.1  19.220  < 2e-16 ***
## GUAR   5061.3      560.7   9.027 2.42e-09 ***
## F_G    2391.8     2447.4   0.977    0.338    
## T_G   -1527.4     2447.4  -0.624    0.538    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1032 on 25 degrees of freedom
## Multiple R-squared:  0.983,  Adjusted R-squared:  0.9796 
## F-statistic: 288.4 on 5 and 25 DF,  p-value: < 2.2e-16

Then T_G

## 
## Call:
## lm(formula = Viscosity_cP ~ 0 + FSC + TARA + GUAR + F_G, data = df_complete)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2177.04  -701.62   -50.49   563.24  2012.96 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## FSC    6609.7      480.2   13.77 1.89e-13 ***
## TARA   9237.0      423.9   21.79  < 2e-16 ***
## GUAR   4886.8      480.2   10.18 1.47e-10 ***
## F_G    2610.1     2393.7    1.09    0.286    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1020 on 26 degrees of freedom
## Multiple R-squared:  0.9827, Adjusted R-squared:   0.98 
## F-statistic: 369.1 on 4 and 26 DF,  p-value: < 2.2e-16

And eventually F_G .

## 
## Call:
## lm(formula = Viscosity_cP ~ 0 + FSC + TARA + GUAR, data = df_complete)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2206.03  -820.70     2.15   706.55  1983.97 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## FSC    6870.7      417.7   16.45 1.36e-15 ***
## TARA   9150.0      417.7   21.90  < 2e-16 ***
## GUAR   5147.8      417.7   12.32 1.35e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1023 on 27 degrees of freedom
## Multiple R-squared:  0.9819, Adjusted R-squared:  0.9799 
## F-statistic: 488.4 on 3 and 27 DF,  p-value: < 2.2e-16

Coefficients and residuals

The coefficients of my model will therefore be

coefficients(model)

##      FSC     TARA     GUAR 
## 6870.677 9149.960 5147.763

And a value of r² of

summary(model)$r.squared

## [1] 0.9819043

I can also draw a scatter plot that relates the response obtained experimentally with the one calculated by the regression model

residui = residuals(model)

df_visc = data.frame('ID' = as.character(df$ID), 'Week' = as.character(df$Week), 'Visc' = df$Viscosity_cP, 'Visc_calc' = df$Viscosity_cP-residui )
bisettrice = data.frame(x= c(4000,10500),y=c(4000,10500))

fig = ggplot(df_visc, aes(x=Visc, y=Visc_calc, color = Week))+ geom_point(aes(label = ID))+
  #geom_text(label=df_visc$id, size=3, hjust=0, vjust=2)+
  geom_line(data=bisettrice, aes(x,y), colour = 'grey')+
  ggtitle('Sperimentali vs calcolati')+
  theme_minimal()
  #theme(plot.title = element_text(hjust = 0.5))

  ggplotly(fig)

---

Response surface

Starting from the coefficients obtained from the model I can visualize the response surface.

fsc = df$FSC
tara = df$TARA
guar = df$GUAR

par(mar = rep(0.2, 4))

TernaryPlot(alab='tara',blab =  'guar',clab =  'fsc')


Fu = function(tara, guar, fsc) {
        model$coefficients['FSC']*fsc+
        model$coefficients['GUAR']*guar+
        model$coefficients['TARA']*tara
        
}

values = TernaryPointValues(Fu)
ColourTernary(values)

TernaryContour(Fu, resolution = 36L)

It is clear that the tara gum is the thickener that gives greater viscosity to the mixture at the same weight while the guar gum is the one that has a lower thickening power.

---

Overrun

Regression

As made for the viscosity, a regression is made to find the optimal coefficients according to the model, evaluating their significance (95% confidence limit).

Being a model relative to a ternary mixture I don’t have the value relative to the intercept as for a traditional linear regression.

---

model = lm(Overrun~ 0+FSC+TARA+GUAR+F_T+F_G+T_G+F_T_G, df_complete)
summary(model)

## 
## Call:
## lm(formula = Overrun ~ 0 + FSC + TARA + GUAR + F_T + F_G + T_G + 
##     F_T_G, data = df_complete)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.120859 -0.028096  0.005774  0.023291  0.098241 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## FSC    0.34487    0.02480  13.908 1.10e-12 ***
## TARA   0.31290    0.02480  12.618 8.05e-12 ***
## GUAR   0.54590    0.02480  22.015  < 2e-16 ***
## F_T    0.14700    0.11591   1.268    0.217    
## F_G   -0.07299    0.11591  -0.630    0.535    
## T_G    0.02117    0.11591   0.183    0.857    
## F_T_G -0.86216    0.85023  -1.014    0.321    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04515 on 23 degrees of freedom
## Multiple R-squared:  0.9907, Adjusted R-squared:  0.9879 
## F-statistic: 351.3 on 7 and 23 DF,  p-value: < 2.2e-16

It is clear that all the coefficients related to the 4 interactions are statistically not significant therefore I remove the variable with the greater P-value and I repeat the procedure (this because the three starting vectors are linearly dependent in this kind of model).

---

model = lm(Overrun~ 0+FSC+TARA+GUAR+F_T+F_G+T_G, df_complete)
summary(model)

## 
## Call:
## lm(formula = Overrun ~ 0 + FSC + TARA + GUAR + F_T + F_G + T_G, 
##     data = df_complete)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.115387 -0.022624 -0.000817  0.025591  0.103713 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## FSC   0.34852    0.02455  14.197 3.57e-13 ***
## TARA  0.31654    0.02455  12.894 2.78e-12 ***
## GUAR  0.54954    0.02455  22.386  < 2e-16 ***
## F_T   0.10596    0.10868   0.975    0.339    
## F_G  -0.11403    0.10868  -1.049    0.305    
## T_G  -0.01987    0.10868  -0.183    0.856    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04518 on 24 degrees of freedom
## Multiple R-squared:  0.9903, Adjusted R-squared:  0.9879 
## F-statistic: 409.2 on 6 and 24 DF,  p-value: < 2.2e-16

Also F_T is removed

model = lm(Overrun~ 0+FSC+TARA+GUAR+F_G+T_G, df_complete)
summary(model)

## 
## Call:
## lm(formula = Overrun ~ 0 + FSC + TARA + GUAR + F_G + T_G, data = df_complete)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.119311 -0.022005  0.002922  0.025843  0.099789 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## FSC   0.36029    0.02135  16.874 3.57e-15 ***
## TARA  0.32832    0.02135  15.376 2.99e-14 ***
## GUAR  0.54954    0.02452  22.408  < 2e-16 ***
## F_G  -0.13169    0.10706  -1.230    0.230    
## T_G  -0.03754    0.10706  -0.351    0.729    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04514 on 25 degrees of freedom
## Multiple R-squared:  0.9899, Adjusted R-squared:  0.9879 
## F-statistic: 491.8 on 5 and 25 DF,  p-value: < 2.2e-16

Then F_G

## 
## Call:
## lm(formula = Overrun ~ 0 + FSC + TARA + GUAR + T_G, data = df_complete)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.120148 -0.018338 -0.003353  0.028751  0.098952 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## FSC   0.34775    0.01895  18.355  < 2e-16 ***
## TARA  0.33082    0.02146  15.413 1.36e-14 ***
## GUAR  0.53449    0.02146  24.902  < 2e-16 ***
## T_G  -0.01872    0.10700  -0.175    0.862    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04558 on 26 degrees of freedom
## Multiple R-squared:  0.9893, Adjusted R-squared:  0.9877 
## F-statistic: 602.5 on 4 and 26 DF,  p-value: < 2.2e-16

And eventually T_G .

## 
## Call:
## lm(formula = Overrun ~ 0 + FSC + TARA + GUAR, data = df_complete)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.122436 -0.018130 -0.003145  0.030207  0.096664 
## 
## Coefficients:
##      Estimate Std. Error t value Pr(>|t|)    
## FSC   0.34838    0.01827   19.07   <2e-16 ***
## TARA  0.32895    0.01827   18.00   <2e-16 ***
## GUAR  0.53262    0.01827   29.15   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04475 on 27 degrees of freedom
## Multiple R-squared:  0.9893, Adjusted R-squared:  0.9881 
## F-statistic: 833.2 on 3 and 27 DF,  p-value: < 2.2e-16

Coefficients and residuals

The coefficients of my model will therefore be

coefficients(model)

##       FSC      TARA      GUAR 
## 0.3483774 0.3289526 0.5326228

And a value of r² of

summary(model)$r.squared

## [1] 0.9893134

I can also draw a scatter plot that relates the response obtained experimentally with the one calculated by the regression model

residui = residuals(model)

df_overrun = data.frame('ID' = as.character(df$ID), 'Week' = as.character(df$Week), 'Overrun' = df$Overrun, 'Overrun_calc' = df$Overrun-residui )
bisettrice = data.frame(x= c(0.25,0.6),y=c(0.25,0.6))

fig = ggplot(df_overrun, aes(x=Overrun, y=Overrun_calc, color = Week))+ geom_point(aes(label = ID))+
  #geom_text(label=df_visc$id, size=3, hjust=0, vjust=2)+
  geom_line(data=bisettrice, aes(x,y), colour = 'grey')+
  ggtitle('Sperimentali vs calcolati')+
  theme_minimal()
  #theme(plot.title = element_text(hjust = 0.5))

  ggplotly(fig)

---

Response surface

Starting from the coefficients obtained from the model I can visualize the response surface.

fsc = df$FSC
tara = df$TARA
guar = df$GUAR

par(mar = rep(0.2, 4))

TernaryPlot(alab='tara',blab =  'guar',clab =  'fsc')


Fu = function(tara, guar, fsc) {
        model$coefficients['FSC']*fsc+
        model$coefficients['GUAR']*guar+
        model$coefficients['TARA']*tara
        
}

values = TernaryPointValues(Fu)
ColourTernary(values)

TernaryContour(Fu, resolution = 36L)

As for the overrun, there seems to be an inverse correlation with the viscosity. // This assumption can be verified by correlating the two response variables.

df2 = data.frame(df[c("ID", "Week", "Viscosity_cP", "Overrun")])

fig = ggplot(df2, aes(x=Viscosity_cP, y=Overrun, color = Week))+ geom_point(aes(label = ID))+
  xlim(2500, 10500) + 
  ggtitle('Viscosity vs Overrun')+
  theme_minimal()  
  #theme(plot.title = element_text(hjust = 0.5))

  ggplotly(fig)