## 1 Introduction  
The classical theory of statistical analysis is predicated on the assumption that model errors have a normal distribution. When the error has a normal distribution, the Ordinary Least Squares (OLS) method can be used to estimate a linear regression model. However, minor differences in the distribution of model errors could significantly affect the resulting inferences [1]. The multivariate t-distribution approach is more suitable than the multivariate Pareto distribution approach when the data has heavy tails and is not bounded below. The multivariate t-distribution can handle such data, while the multivariate Pareto distribution is only suitable for non-negative variables [2]. Furthermore, the multivariate t-distribution can be utilized for a robust estimate of means, regression coefficients, and variance-covariance matrices in multivariate linear models, even in missing data [3]. The multivariate t-distribution approach matches various data conditions examined in several studies [4, 5, 6]. Although these studies use the t-distribution approach, they do not use data obtained from the experimental design or for the second-order model. The Maximum Likelihood Estimation (MLE) method with the Expectation Maximization (EM) algorithm can be used to estimate the linear regression model parameters with the t-distribution approach [7]. It is recommended to treat the t-distribution's degrees of freedom as a known parameter if the sample size is small. If the sample size is large, estimating the degree of freedom is recommended. It is demonstrated that the EM algorithm has a relatively high computing efficiency and that, from an analytical standpoint, it is pretty straightforward and converges monotonically to an approximated local maximum [3]. The multivariate t-distribution approach matches various data conditions, and utilizing the EM algorithm makes this approach a simple and decent alternative in practice to get robust results [5]. These previous studies show that the estimator obtained using the t-distribution approach is more efficient than that obtained under the normality assumption when the data is in the form of a heavy tail. After obtaining the estimator parameter MLE, the next step is using the Likelihood Ratio Test (LRT) method to test the parameters simultaneously. These methods have also been used in several other studies related to modeling using a non-normal distribution approach [8; 9]. The Response Surface Methodology (RSM) was first introduced by [PERSON] and [PERSON] [10]. RSM includes experimental strategies to explore experimental regions, empirical statistical modeling to obtain an approach to the relationship between responses and input variables (factor levels), then optimization methods to find the level of factors that optimize response values [11]. In this paper, data collection used a Central Composite Design (CCD) for the second-order model strategy. Then, the main objective is to develop a Second-Order Multiresponse Surface (SOMRS) model with a multivariate t-distribution approach and the tools for statistical inference, including parameter estimation and hypothesis testing. Since the RSM is intended to determine the optimal conditions of the system, the optimization technique used in this paper is the Dual Response System (DRS). In the DRS technique, one of the responses (considered the most important) is determined as a primary response, while other responses are secondary. Furthermore, this primary response is optimized, subject to conditions that have been determined from the secondary response [12; 13] and using the mean squared error (MSE) criterion [14]. RSM has been used in several studies in the pharmaceutical industry, especially in optimizing the product's quality characteristics.