Journal of Dairy Science and Biotechnology
Korean Society of Dairy Science and Biotechnology
REVIEW

우유 균질 조건 예측을 위한 반응표면방법론의 활용

임성수1https://orcid.org/0000-0001-7009-3343, 오세종2,*https://orcid.org/0000-0002-5870-3038
Sungsue Rheem1https://orcid.org/0000-0001-7009-3343, Sejong Oh2,*https://orcid.org/0000-0002-5870-3038
1고려대학교 빅데이터사이언스학부
2전남대학교 동물자원학부
1Division of Big Data Science, Korea University, Sejong, Korea
2Division of Animal Science, Chonnam National University, Gwangju, Korea
*Corresponding author : Sejong Oh, Division of Animal Science, Chonnam, National University, Gwangju, Korea, Tel : +82-62-530-2116, Fax : +82-62-530-2129, E-mail : soh@jnu.ac.kr

© Copyright 2023, Korean Society of Dairy Science and Biotechnology. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Received: Feb 16, 2023; Revised: Feb 22, 2023; Accepted: Feb 23, 2023

Published Online: Mar 31, 2023

Abstract

Response surface methodology (RSM) is a statistical approach widely used in food processing to optimize the formulation, processing conditions, and quality of food products. The homogenization process is achieved by subjecting milk to high pressure, which breaks down fat globules and disperses fat more evenly throughout milk. This study focuses on an application of RSM including the logit transformation to predict the efficiency of milk homogenization, which can be maximized by minimizing the relative difference in fat percentage between the top part and the remainder of milk. To avoid a negative predicted value of the minimum of this proportion, the logit transformation is used to turn the proportion into the logit, whose possible values are real numbers. Then, the logit values are modeled and optimized. Subsequently, the logistic transformation is used to turn the predicted logit into the predicted proportion. From our model, the optimum condition for the maximized efficiency of milk homogenization was predicted as the combination of a homogenizer pressure of 30 MPa, a storage temperature of 10°C, and a storage period of 10 days. Additionally, with a combination of a homogenizer pressure of 30 MPa, a storage temperature of 10°C, and a storage period of 50 days, the level of milk homogenization was predicted to be acceptable, even with the problem of extrapolation taken into account.

Keywords: response surface methodology; optimization; logit transformation; logistic transformation; homogenization

서 론

균질공정은 식품 산업에서 많이 사용되는 공정으로 혼합된 첨가물을 균일하게 분산시키는 역할을 수행한다. 전세계적으로 시판되는 우유 대부분은 균질화되고 열처리된 이후에 포장해서 판매된다. 균질공정은 원유내 함유된 지방구의 크기를 물리적으로 작게 하여 상층부로 지방구가 부상하여 크림층을 형성하는 것을 방지하고자 실시한다. 균질 공정 후 우유의 지방구는 직경 3–5 µm에서 직경 1 µm 이하로 미세화되어서 우유 상층부로 부상하는 현상이 일어나지 않는다. 우유 균질기는 1900년 파리세계박람회에서 처음 소개되었는데 우유와 크림의 흐름을 작은 관을 통해 밀어낼 수 있도록 설계된 3개의 피스톤 펌프로 구성되어 있었다.

Gaulin 균질기는 현재까지 시간에 따라 분리되는 성질의 액체를 균일하게 유화시키는 가장 간단한 기계장치로 사용되고 있다. 우유 균질 조건은 약 60°C 정도에서 실시하지만 압력과 온도 조건은 장치와 밸브 유형에 따라 다르다. 균질화 효율은 온도가 42°C–72°C일 때 증가하고 72°C–77°C 부근에서 안정화된다. 일반적으로 우유 균질 온도는 55°C–80°C의 범위이며, 균질 압력은 10–25 MPa의 범위이다[1,2].

반응표면방법론(response surface methodology, RSM)은 실험 설계와 최적화를 위한 통계적 방법으로, 실험 횟수를 최소화하면서 최적조건을 찾을 수 있다. 실제로 여러 요인들이 우유 품질에 미치는 영향을 조사하고 최적 조건을 찾는 데에 활용될 수 있다. 반응표면방법론을 사용하여 우유 품질에 미치는 영향을 조사한 연구는 해외에서 많이 보고되고 있으나, 국내에서는 Kang et al.[3]의 연구가 유일한 것이다. 이 연구에서는 균질 압력, 저장온도 및 저장기간에 따른 우유의 지방구 크기와 크림화 정도 등을 반응변수들로 해서 균질 효율에 대한 예측 모형을 구축하였는데, 그 결과 우유를 10°C에서 15일의 저장 기간 동안 적절한 균질 효율을 유지하기 위해서는 16.8 MPa의 균질 압력이 필요한 것으로 나타났다[3].

2023년도부터 시행하고 있는 식품 소비기한 표시제도는 소비자에 식품을 판매해도 되는 최종기간을 표기하는 제도이다. 식품폐기물을 줄이기 위한 방편으로 식품에 표기되던 유통기한을 소비기한으로 변경하는 것으로 2021년 관련법이 개정되고 준비기간을 거쳐 새롭게 시행되는 제도이다. 냉장식품의 경우에는 품목에 따라 유예기간을 두었는데 우유류는 2031년부터 시행된다. 식품의약품안전처에서는 식품유형별 설정 실험 결과를 예시하였는데, 우유류의 경우 현행 15일의 유통기한이 평균 24일로 권장하고 있으나 제품 유형별 유통기한은 생산업체가 자율로 결정한다. 소비기한 설정에는 미생물 및 이화학적 검사를 통한 우유의 안전성이 보장되어야 하고 관능검사에 따른 품질의 변화가 없는 기한을 고려해야 한다. 우유의 균질은 품질과 관련이 있는데, 균질이 잘 된 우유는 유성분들이 안정적으로 분산되어 있어 상미기간동안 우유의 맛의 변화가 거의 일어나지 않는다. 그러나 유통기한이 길 경우에 현행의 균질조건에서 시간이 지남에 따라 균질된 지방구가 부상할 수 있기 때문에 균질화 정도를 파악하는 것은 중요하다.

따라서, 본 논문은 Kang et al.[3]이 실험한 연구자료를 바탕으로 보다 통계학적으로 발전된 모형을 구축하였으며, 이 모형을 사용하여 우유를 저장온도 10°C에서 30 MPa의 균질압력으로 50일까지 저장하였을 때의 균질효율 측정변수의 값을 예측하고자 하였다.

본 론

1. 반응표면방법론 활용 조건

지금까지 사용되고 있는 최적화 방법으로는 탐색적 방법, EVOP(evolutionary operation), 수학적 기법 등을 들 수 있으며, SAS, Minitab 등의 컴퓨터 프로그램이 최적화를 수행할 수 있다[4,5].

유가공학을 포함하는 식품공학 및 생물공학 분야에서 폭넓게 사용되는 반응표면방법론(RSM)은 실험의 설계와 분석을 통한 반응의 모형화와 반응 및 요인수준 조합의 최적화를 위한 통계적 방법들의 집합이다. 반응표면방법론에 사용되는 실험설계로는 주로 중심 합성 설계(central composite design)[6]가, 그리고 최적화를 위한 분석모형으로는 주로 2차 다항 모형(second-order polynomial model)이 사용되어 왔다[7].

반응을 함수적으로 이해하고 최적화를 수행하기 위하여, RSM에서는 적절한 모형의 사용이 매우 중요한데, 여기서 적절한 모형은 다음의 네 가지 기준들을 모두 충족시키는 모형을 의미한다:

  • (1) 모형이 5% 수준에서 유의함(모형의 p-value≤0.05).

  • (2) 적합 결여(lack of fit)가 5% 수준에서 유의하지 않거나(적합 결여 p-value>0.05) 적합 결여가 유의하더라도, 적합 결여 통계량의 분모인 순수오차 평균제곱(pure error mean square)이 매우 작은 경우에는, 적합 결여 유의성의 원인이 동일 실험 조건에서의 반응치들의 값들이 아주 비슷한 현상이어서, 적합 결여의 유의성이 문제가 되지 않음[8].

  • (3) 수정 결정계수가 0.8 이상임(adjusted R-square≥0.8)[7].

  • (4) 반응의 최적화에서 예측된 반응 최적치가 발생 가능한 값이어야 함.

그러나, 반응값이 비율로 표시되어 그 가능한 값의 범위가 0 이상 1 이하인 경우, 최적화의 목적이 최대화일 때 예측된 반응 최적치로 1보다 큰 값이 나오거나 최적화의 목적이 최소화일 때 예측된 반응 최적치로 0보다 작은 값이 나오는 상황이 발생될 수 있다. 이러한 상황은 반응값인 비율을 로짓 변환(logit transformation)하여 모형화를 실행함으로써 방지될 수 있다. 로짓 변환에 대해서는 3절에서 설명한다.

2. 균질효율과 분석대상 자료

균질 공정은 우유를 일반적으로 약 2,500–3,000 psi의 고압에 노출시켜 우유의 지방구를 더 작은 입자로 미세화 시키는 것으로 지용성 비타민 A와 같은 특정 영양소의 생물학적 가용성을 증가시켜 우유의 영양적 품질을 향상시키고 바람직한 조직감(texture)과 관능특성을 유지하는 데 도움을 준다.

우유 균질 효율의 평가를 위하여 본 논문에서는 미국 공중보건국(United States Public Health Service, USPHS)에서 제시하는 평가방법을 사용하는데, 이 방법은 우유 상층부와 나머지 부분 간의 지방 비율에 있어서의 상대적 차이(relative difference in fat percentage, RDIF)로 표시한다. RDIF의 계산식은 1L 시료에서 상층부 100 mL의 지방 비율을 Ft라고 하고 나머지 부분의 지방 비율을 Fb라고 할 때, RDIF=[(Ft–Fb)/Ft]×100%이다[3]. RDIF는 백분율로 표기되며, 그 값이 낮을수록 우유 균질화가 잘된 것으로 통상적으로 RDIF 값이 10% 이하면 균질 수준이 적절한 것으로 인정된다[3].

Kang et al.[3]의 연구에서는 균질 압력, 저장 온도 및 저장 기간을 설명요인들, 즉 독립변수들로 하고 균질화 효율에 대한 여러 측정값들을 반응변수, 즉 종속변수로 하여 예측 모형을 구축하였다.

Kang et al.[3]이 사용한 반응변수들 중에 가장 의미있는 변수가 퍼센트 값으로 표시되는 RDIF이므로 본 논문에서는 RDIF를 100으로 나누어 0과 1 사이의 값(P)으로 재표현한 다음에, 소수(decimal)로 표기되는 비율 P의 함수로서 “logit(P)=log(P/[1–P])”로 정의되는 logit(로짓) 함수를 사용하여 P를 변환시켜 모형의 반응변수로 사용하였다.

Kang et al.[3]이 실시한 실험 시행(run)에서 퍼센트 값으로 표기된 RDIF, 소수로 표기되는 P, 그리고 P의 logit 함수 값을 Table 1에 제시하였으며, Table 1에서 X1=–1, 0, 1은 균질 압력=10, 20, 30 MPa, X2=–1, 0, 1은 저장 온도=10°C, 15°C, 20°C, X3=–1, 0, 1은 저장 기간=10, 15, 20 days를 각각 뜻한다.

Table 1. Experimental design in coded factors and responses
Case Design point Run order X1 X2 X3 RDIF (%, from USPHS code)* P=RDIF/100 L=logit(P)=log(P/[1–P])
1 1 2 –1 –1 –1 30.12 0.3012 –0.84159
2 2 5 –1 –1 1 36.36 0.3636 –0.55977
3 3 20 –1 1 –1 48.69 0.4869 –0.05241
4 4 10 –1 1 1 59.62 0.5962 0.38966
5 5 11 1 –1 –1 0.73 0.0073 –4.91255
6 6 15 1 –1 1 1.96 0.0196 –3.91243
7 7 13 1 1 –1 4.09 0.0409 –3.15487
8 8 17 1 1 1 8.55 0.0855 –2.36986
9 9 1 –1 0 0 38.91 0.3891 –0.45110
10 10 7 1 0 0 2.63 0.0263 –3.61153
11 11 4 0 –1 0 3.55 0.0355 –3.30208
12 12 16 0 1 0 15.73 0.1573 –1.67846
13 13 18 0 0 –1 5.65 0.0565 –2.81536
14 14 6 0 0 1 7.84 0.0784 –2.46429
15 15 3 0 0 0 8.55 0.0855 –2.36986
16 15 8 0 0 0 6.99 0.0699 –2.58823
17 15 9 0 0 0 6.56 0.0656 –2.65633
18 15 12 0 0 0 6.78 0.0678 –2.62099
19 15 14 0 0 0 5.88 0.0588 –2.77301
20 15 19 0 0 0 6.98 0.0698 –2.58977

X1=–1, 0, 1 stands for homogenizer pressure=10, 20, 30 MPa; X2=–1, 0, 1, storage temperature=10°C, 15°C, 20°C; and X3=–1, 0, 1, storage period=10, 15, 20 days.

The lower RDIF is, the better homogenization efficiency is.

USPHS, United States Public Health Service.

Download Excel Table
3. Logit 변환을 사용한 2차 모형에 대한 평가 및 모형 구축

0과 1 사이의 값을 갖는 비율 P에 대한 변환으로는 logit 변환이 자주 사용되는데[9], logit 함수는 다음과 같이 표기된다.

log it ( P ) = log ( P / [ 1 P ] )

여기서 log는 ln으로도 표기되는 자연 로그로 그 밑은 e=2.71828…이다. 비율(P)과 비율의 logit 함수 (logit[P]) 간의 관계를 나타내는 그래프를 SAS와 Excel을 이용하여 Fig. 1A로 제시하였다. Fig. 1A에서 가로축의 P는 0.001, 0.002, …, 0.998, 0.999의 값들을 갖고, 세로축의 logit은 대략 –6.9와 6.9 사이의 값들을 갖는다. 비율 P의 범위가 0과 1 사이인 것에 비해, logit(P)의 범위는 이론상으로는 실수 구간 전체이어서, logit(P)를 2차 다항 모형으로 표현할 경우, 그 모형이 산출해 내는 값에 제한이 없는 것이 허용된다. 따라서, logit 함수의 값을 반응치로 사용하여, 다음과 같은 2차 모형을 추정하였다. 아래 모형식에서, c0는 절편(intercept)을 뜻하고, c1, c2, …, c23는 각 모형 항(model term)의 계수(coefficient)를 의미한다.

jmsb-41-1-1-g1
Fig. 1. Logit and logistic functions. (A) Logit function of the proportion P: Logit=log(P/[1–P]); (B), Logistic function of the logit L: P=exp(L)/(1+exp[L]).
Download Original Figure
Predicted logit = c 0 + c 1 X 1 +c 2 X 2 +c 3 X 3 +c 11 X 1 2 +c 22 X 2 2 +c 33 X 3 2 +c 12 X 1 X 2 +c 13 X 1 X 3 +c 23 X 2 X 3

그리고, 이 2차 모형으로부터의, logit으로 표현된 최적 반응치를 비율 값으로 돌려놓는 것이다[9]. 이를 위해서는, P를 비율이라고 할 때, 다음 공식이 사용된다.

P=e logit / [ 1 + e logit ] = exp ( logit ) / [ 1 + exp ( logit ) ]

위 공식이 나타내는 함수를 logistic(로지스틱) 함수라고 하며, logit과 비율 간의 관계를 나타내는 logistic 함수의 그래프를 SAS와 Excel을 이용하여 Fig. 1B로 제시하였다. Fig. 1B에서 가로축의 logit은 –7.00, –6.99, …, 6.99, 7.00의 값들을 갖고, 세로축의 P는 0과 1 사이의 값들을 갖는다.

Logit 값들을 얻기 위하여 본 논문에서는 먼저 퍼센트 값인 USPHS code를 100으로 나누어 소수로 표현된 비율 P를 얻고 이 P로부터 logit=log(P/[1-P])의 값을 구하였다. 이렇게 얻어진 logit을 모형의 반응변수로 하여 2차 다항 모형을 적합시켰다.

Logit 변환 후 얻어진 2차 모형은 분산분석 결과, 적합한 모형의 기준들을 모두 충족한 것으로 나타났다.

  • (1) 모형이 5% 수준에서 유의(모형의 p-value=0.0001 미만<0.05): 기준 충족.

  • (2) 적합 결여가 5% 수준에서 유의하지 않음(적합 결여 p-value=0.3930>0.05): 기준 충족.

  • (3) 수정 결정계수가 0.8 이상임(adjusted R-square=0.9892>0.8): 기준 충족.

Logit에 대한 이 2차 모형은, 기준 (1), (2), (3)을 충족시키고 R-square도 0.9943으로 매우 높아서, 적합한 모형인 것으로 보인다. 이 모형을 식으로 표기하면, 다음과 같다.

Predicted logit = 2.605 1.645 X 1 + 0.666 X 2 + 0.286 X 3 + 0.581 X 1 2 + 0.122 X 2 2 0.027 X 3 2 + 0.195 X 1 X 2 + 0.133 X 1 X 3 0.007 X 2 X 3

Table 2는 이 logit에 대한 2차 모형에서의 coefficient estimate, SE, t-value, p-value를 담고 있다. Table 2에 의하면, 5% 유의수준에서, 1차항들은 모두 유의하고, 2차항들 중에서는 X12 항만 유의하며, 교차항들 중에서는 X1 X2 항과 X1 X3 항만 유의하다. Table 2t-value, p-value들을 관찰하면, 요인 영향력의 크기 순서는 X1, X2, X3임을 알 수 있다.

Table 2. Coefficient estimates and related statistics in the logit model
Model term Coefficient estimate SE t-value p-value
Intercept –2.6048 0.0485 –53.75 0.000
X1 –1.6446 0.0446 –36.89 0.000
X2 0.6662 0.0446 14.94 0.000
X3 0.2860 0.0446 6.42 0.000
X12 0.5812 0.0850 6.84 0.000
X22 0.1223 0.0850 1.44 0.181
X33 –0.0273 0.0850 –0.32 0.755
X1 X2 0.1952 0.0498 3.92 0.003
X1 X3 0.1327 0.0498 2.66 0.024
X2 X3 –0.0069 0.0498 –0.14 0.893
Download Excel Table
4. 측정치에 대한 logit(로짓) 변환을 사용한 2차 모형에 의한 최적화

Minitab의 반응 최적화 도구를 이용하여 logit의 최적화 목적을 최소화로 지정해서 앞에서 추정된 2차 다항 모형에 근거한 최적화를 실시하면, logit을 최소화하는 요인수준 조합은 (X1, X2, X3)=(1, –1, –1)로 추정되는데, 이 요인수준 조합에서 예측된 logit 값이 Predicted logit=–4.8602로 계산된다. Logit 함수가 취할 수 있는 값에는 제한이 없으므로, 이 값 –4.8602는 받아들일 수 있는 값이고, 이제 위에 주어진 logistic 함수를 사용하여, 아래에서와 같이, 이 logit 값 –4.8602를 비율 값으로 변환시킨다.

P = e 4.8602 / [ 1 + e 4.8602 ] = exp ( 4.8602 ) / [ 1 + exp ( 4.8602 ) ] = 0.0077

이 비율 값 0.0077은 비율의 최소치로 충분히 받아들일 수 값이다. 그리고, 코드화된 최적조건 추정점 (X1, X2, X3)=(1, –1, –1)에 대응되는 실제 최적조건 추정점은 균질 압력 30 MPa, 저장 온도 10°C, 저장 기간 10 days이다. 이 조건은 상층부 지방함량과 나머지 부분 지방함량의 상대 차이가 0에 가까워지는, 가장 최대의 균질 효율을 보이는 조건으로 추정되는데, 이 조건에서의 실제 확인실험을 실시해서 이를 검증할 필요가 있다.

그런데, 우리의 경우, 코드화된 최적조건 추정점 (X1, X2, X3)=(1, –1, –1)은 이미 실험설계에 들어 있는 설계점(design point)이므로, 확인실험을 추가적으로 시행할 필요 없이, Table 1에서 (X1, X2, X3)=(1, –1, –1)인 설계점에서의 Y 값을 확인함으로써, 위에 제시된 실제 최적조건 추정점의 타당성을 평가할 수 있을 것이다.

Table 1에서 (X1, X2, X3)=(1, –1, –1)인 설계점은 design point 5이고, 여기서의 Y 값은 0.73%=0.0073으로 매우 작은 값으로 모형으로부터 예측된 비율 최소치 0.0077에 무척 가까워서, 본 연구에서 도출된 코드화된 최적조건 추정점 (X1, X2, X3)=(1, –1, –1)은 Y 값의 최소화를 성취해내는 최적 조건으로의 타당성이 인정된다고 판단된다.

코드화된 요인 2개씩 짝을 지어 각 가로축에 배치하고 logit 예측치를 세로축에 배치한 3차원 표면도(3-dimensional surface plot)들을 Fig. 2에 나타내었다. 여기서 각 그림에 표시되지 않은 X 요인의 값은 0으로, Fig. 2에서 보는 바와 같이, logit 값이 최소화되는 점은 (X1, X2, X3)=(1, –1, –1)임을 알 수 있다.

jmsb-41-1-1-g2
Fig. 2. Three-dimensional surface plots of the logit model for each pair of coded factors, where X1=–1, 0, 1 stands for homogenizer pressure=10, 20, 30 MPa; X2=–1, 0, 1, storage temperature= 10°C, 15°C, 20°C; and X3=–1, 0, 1, storage period=10, 15, 20 days.
Download Original Figure

그러나, 이 조건은 상층부의 지방함량과 나머지 부분의 지방함량의 차이가 거의 0이 되는 최적조건으로 우유가 거의 100% 균질하게 되는 조건을 제시한 것으로 실제 공정에서 우유가 100% 균질성을 갖는 것보다 저장기간 10일 이후에서의 균질효율의 변화가 중요하다고 판단된다.

여기서 얻어진 모형에서 우유의 저장기간을 50일까지 연장시켰을 때의 RDIF 예측값들의 변화 추이를 Fig. 3에 나타내었다. Kang et al.[3]은 실제 실험에서 저장 20일까지 조사하였으나, 저장 20일을 넘어 50일까지로 기간을 확장시켜 예측한 결과, 우유를 30 MPa로 균질한 경우 저장온도 10°C에서의 RDIF의 예측값은 45일 저장 시에는 5.54%로, 50일 저장 시에는 5.92%로 계산되었다. 실험에서의 저장 기간의 최대치인 20일을 넘어선 기간에서의 예측에 따른, 외삽(extrapolation)에 의한 예측의 불안정성을 감안하더라도, 저장기간 45일과 50일에서의 RDIF의 예측값 5.54%와 5.92%는 적절한 RDIF 값의 한계치인 10%보다 훨씬 아래에 있으므로, 저장 50일까지 우유내 지방은 적절한 균질 상태를 유지할 것으로 보인다.

jmsb-41-1-1-g3
Fig. 3. Predicted values of RDIF against milk storage periods. Homogenization efficiency is acceptable if RDIF is no higher than 10%.
Download Original Figure

결 론

본 논문은 우선 비율 반응변수의 logit 변환을 통한 반응표면모형 구축 방법론, 즉 비율로 표시되는 반응값을 logit 함수를 사용하여 logit 값들로 변환시키고, 이 logit 값들에 대한 모형화를 실시하여 그 모형으로 최적화를 실시한 후에 그 최적화에서의 logit 값을 logistic 함수를 이용하여 비율 값으로 변환시켜서 그 비율 값을 보고 최적화의 성취 정도를 평가하는 방법론을 제시하였고, 이 방법론을 활용하여 우유 균질 조건의 예측을 가능하게 하는 반응표면모형을 구축하였으며, 이 모형을 통하여 최적 균질압력, 최적 저장온도 및 최적 저장기간을 추정하였고, 최적 균질압력과 최적 저장온도에서의 장기간 저장 시의 우유 균질화 정도의 적절성을 평가하였다.

Kang et al.[3]의 연구에서의 실험자료로부터의 logit 모형을 활용한 결과, 30 MPa의 균질압력으로 우유를 균질하여 10°C에서 50일까지 저장해도, 상층부와 하층부 간의 지방함량의 상대 차이가 5.92% 정도로 크지 않게 예측되어, 우유의 균질 정도는 적절한 수준을 유지할 것으로 보인다.

Conflict of Interest

The authors declare no potential conflict of interest.

감사의 글

본 논문은 2022년도 고려대학교 공공정책대학 교내지원연구비에 의해 연구되었음.

References

1.

Trout GM. The nutritive value of homogenized milk: a review. J Dairy Sci. 1948;31: 627-655.

2.

Michalski MC, Januel C. Does homogenization affect the human health properties of cow's milk? Trends Food Sci Technol. 2006;17:423-437.

3.

Kang HJ, Kang SH, Shin YK. Prediction of homogenization efficiency using response surface methodology. J Dairy Sci Biotechnol. 2017;35:202-207.

4.

Minitab. Minitab 19. State College, PA: Minitab; 2019.

5.

SAS. SAS 9.4. Cary, NC: SAS Institute; 2016.

6.

Box GEP, Wilson KB. On the experimental attainment of optimum conditions. J R Stat Soc Series B Methodol. 1951;13:1-38.

7.

Myers RH, Montgomery DC, Anderson-Cook CM. Response surface methodology: process and product optimization using designed experiments. 3rd ed. Hoboken, NJ: John Wiley & Sons; 2009.

8.

Rheem S, Rheem I. The model-fit proportion and the coefficient of pure error variation, the measures to supplement the lack-of-fit test in response surface analysis. J Korean Data Anal Soc. 2012;14:2989-2994.

9.

Lee YJ, Rheem S. Using a logit model to estimate and compare proportions in the presence of Simpson’s paradox. J Korean Data Anal Soc. 2014;16:3037-3047.