Predictive Modelling of Critical Variables for Improving HVOF Coating using Gamma Regression Models

Thermal spray coating is a critical process in many industries, involving the application of coatings to surfaces to enhance their functionality. This paper proposes a framework for modelling and predicting critical target variables in thermal spray coating processes, based on the application of statistical design of experiments (DoE) and the modelling of the data using generalized linear models (GLMs) with a particular emphasis on gamma regression. Experimental data obtained from thermal spray coating trials are used to validate the presented approach, demonstrating that it is able to accurately model and predict critical target variables. As such, the framework has significant potential for the optimization of thermal spray coating processes, and can contribute to the development of more efficient and effective coating technologies in various industries.


Introduction
Thermal spraying is a surface modification process that involves the deposition of a coating material onto a substrate by heating and accelerating a feedstock material through * voestalpine Stahl GmbH, voestalpine-Straße 3, A-4020 Linz, Austria (wolfgang.rannetbauer@voestalpine.at), Corresponding author.
Numerous techniques have been employed to model and forecast the properties of coatings produced via HVOF spraying.These methods include empirical models based on regression analysis, mechanistic models that simulate the physical and chemical processes occurring during spraying [8,25,32], and hybrid models that integrate both approaches [34].Linear and nonlinear regression models have been used to establish a relationship between process variables and coating properties [18,24,30,33], while computational fluid dynamics (CFD) models have been utilized to simulate the gas flow, heat transfer, and particle behavior in the spray gun [12,19].In addition, artificial neural networks (ANNs) [4,20,21,37] and genetic algorithms (GAs) [17] have been implemented to optimize process conditions and predict coating properties.
Despite the notable progress, the prediction of coating properties is still a challenging task, due to the complex interactions among the process variables, material properties, and the microstructure of the coatings.This study proposes a novel approach for modelling HVOF coatings through systematic variation of process variables that have received relatively limited attention in prior research.To achieve this, a Central Composite Design (CCD) of experiments is employed, which enables efficient exploration of a vast parameter space.The subsequent analysis focuses on developing gamma regression models derived from generalized linear models (GLMs), which are particularly well-suited for modeling data with skewed, non-negative distributions.By integrating these adjusted process variables into the models, a promising opportunity arises to identify novel associations between process conditions and coating properties.
Our approach provides insights into the intricate HVOF process, improving predictive models for key coating characteristics.The framework presented in this study supports the development of efficient coating technologies with enhanced attributes like wear resistance, corrosion protection, and oxidation resilience.These advancements have practical applications in industries such as aerospace, automotive, and manufacturing.
The structure of this paper is as follows: Section 2 presents a comprehensive overview of the HVOF process, including a detailed explanation of the main factors influencing coating properties and particle in-flight characteristics.Section 3 introduces the powerful application of generalized linear models (GLMs) and maximum likelihood estimation (MLE) to model and accurately estimate the dependence of coating properties on process conditions.This section also provides an overview of the theoretical foundations related to the asymptotic properties of MLE and hypothesis testing.The assessment of the predictive performance of the proposed GLMs is presented in Section 4. Section 5 examines the statistical design of experiments (DoE) and central composite design (CCD) as a potent approach for efficient data collection on various levels of factors in the thermal spray coating process.The findings of the proposed framework applied to experimental data is presented in Section 6, with a specific focus on the precise prediction of critical target variables.Finally, Section 7 concludes the paper with a discussion on the effectiveness and potential of the proposed framework for predicting critical target variables in thermal spraying processes.It highlights the contribution of the framework to the development of more efficient coating technologies.

Technical Background
Thermal spraying is a versatile and widely used surface engineering technique that involves the deposition of coatings on the surface of a substrate to enhance its functional properties, such as wear resistance, corrosion resistance, and thermal insulation.The thermal spray coating process typically involves the application of thermal and kinetic energy to induce partial liquefaction of the coating material, thereby accelerating its projection towards the substrate surface.The amount of thermal and kinetic energy depends on the thermal spray coating technique.Various techniques, such as flame spraying, plasma spraying, arc spraying, and high-velocity oxygen fuel (HVOF) spraying can be used for coating using different types of coating material such as powder or wire.
In this work we focus on the gas-fuel HVOF technology, which is described in more detail below.
HVOF thermal spraying has gained significant attention in recent years due to its ability to produce high-quality coatings with superior mechanical and chemical properties [15].The gas-fuel HVOF process creates its thermal energy by combustion of a mixture of oxygen and fuel gas, typically propane, methane or hydrogen, in a highpressure chamber [11].The kinetic energy in the thermal spray process is created by its specific geometric convergent-divergent nozzle design which accelerates the gas stream to supersonic velocities.The kinetic energy is transferred to the sprayed material particles, causing them to partially melt and deform upon impact with the substrate.This results in coatings with high density and superior adhesion [7], almost independent of the thermal spray material composition (metallic, ceramic, cermet).Figure 2.1 shows a schematic depiction of the HVOF process, where the distinct symmetric structure of the torch becomes evident [8].The diagram reveals an axial axis of symmetry, i.e., a balanced and mirrored arrangement of components along a horizontal line within the HVOF system.In the combustion chamber, a defined but variable fuel gas reacts with oxygen, leading to combustion.Compressed air is used as a shroud gas and forms a cover inside the nozzle and outside in the spray plume.The powder coating material is inserted axially and is fed by using nitrogen as a carrier gas.Every characteristic of the gases and the powder, as well as the process attributes, influences combustion and consequently impacts the kinetic and thermal energy transferred to the spray particles.Influencing factors are for example: • Characteristics of gases (e.g., pressure, flow rate, temperature), • Additional process variables (e.g., the ratio of the combustion gases, powder feed rate, spraying distance), • Powder characteristics (e.g., size, chemical composition, density).
The influence of these factors on the thermal and kinetic energy of the spray material can be measured by using in-situ camera equipment that is able to simultaneously detect the temperature and velocity of the sprayed particles.Process input conditions not only affect the particle characteristics but also the performance of the process, as well as the quality of the coating.Relevant measures for the performance of coating processes are • Deposition rate, • Deposition efficiency.
These performance indicators are especially important when taking into account economic aspects of the thermal spray coating process.However, the coating characteristics are the most important properties, since they influence the industrial performance.The desirable properties are, e.g., • Specific porosity of coating, • Specific hardness of coating, • Specific phase and chemical composition, • Specific thickness of coating.
In order to obtain best wear resistance, corrosion resistance, or thermal insulation, a specific combination of coating properties is essential.Additionally, proper surface preparation prior to spraying is necessary to ensure maximum adhesion strength and achieve the desired coating performance characteristics.The formulation of a robust mathematical relationship that links the input conditions controlling the spraying process, the dynamics of particles during flight, and the resulting coating characteristics is essential.Such a correlation not only facilitates a deeper understanding of the fundamental physical mechanisms underlying thermal spraying but also enables the optimization of deposition conditions to achieve the desired coating properties.
The accurate prediction of coating properties remains a challenge, primarily due to the complex and non-linear nature of the relationships between process attributes and coating properties.Achieving high accuracy in property prediction is often elusive, considering the multifaceted interactions at play.Hence, the demand for a reliable and precise prediction model becomes apparent, as it can serve as a catalyst for optimizing the process and elevating the overall quality of the resulting coatings.

Predictive Modelling of HVOF Coating Properties
The following section is dedicated to the derivation of mathematical models that enable the prediction of coating properties for the HVOF process.For this, we propose the use of Generalized Linear Models (GLMs) along with Maximum Likelihood Estimation (MLE) as an effective approach for modelling and estimating the expected values of target variables, conditioned on the explanatory variables (= process variables) in the HVOF process.We adopt the theoretical framework and notation of [10], to develop the statistical model used in this study, which can be expressed by the following general equation: where µ = (µ i ) n i=1 is a vector denoting the conditional mean of the response variable y = (y i ) n i=1 , i.e., the coating properties of interest.The coefficient vector β = (β 0 , β 1 , . . ., β k ) encodes the effects of the explanatory variables (x 1 , . . ., x k ), i.e., the potentially influential process conditions with observations The mean vector µ is related to the linear combination x i T β of the process input attributes by a oneto-one mapping g(•), which is often referred to as the link function [9].In regression analysis, the assumption of additive random errors allows for the decomposition of the response y i into a systematic component E(y i |x i ) and a random component ϵ i , yielding the equation: where the measurement error ϵ i is assumed to be independent of the covariates.The primary objective of regression analysis is to use the data (y i , x i ) n i=1 to estimate the systematic component E(y i |x i ).
Based on the broad background presented in Section 2, the use of Bayesian generalized linear models to model coating properties appears to be a logical choice.The Bayesian framework offers distinct advantages over classical inference by accommodating prior knowledge regarding model parameters.This facilitates the incorporation of insights into the effects of process input variables, enhancing parameter estimation accuracy.Even in the absence of specific information, non-informative priors like the uniform distribution or Jeffreys's prior can be specified [16].Nevertheless, a classical frequentist statistical approach is adopted in this study to develop a model that characterizes the relationship between input variables and coating properties because the knowledge of suitable priors for Bayesian modeling is absent in this context, and the use of non-informative priors does not offer substantial advantages over classical approaches.Subsequent investigations may explore the application of Bayesian GLMs in future research.

Generalized Linear Models (GLMs)
GLMs, as defined in [23], include a broad range of useful statistical models and serve as a powerful tool for data analysis in various fields such as engineering, physics, and biology.They extend the concept of linear regression to handle non-normal response variables, such as binary or count data, by introducing a link function that relates the mean of the response variable to the linear predictor.Here, the effectiveness of GLMs in analyzing data from the HVOF process is explored, aiming to establish a comprehensive model equation that effectively captures the conditional dependence of coating characteristics on process input variables.
In the context of GLMs, the response vector y = (y 1 , y 2 , ..., y n ), is modeled as a vector that follows any distribution from the exponential family, where each element y i is distributed with a mean µ i and variance σ 2 i .To model the relationship between the response variable y and the predictor variables x i = (x i 1 , x i 2 , . . ., x i k ), a predictor vector η = (η 1 , η 2 , ..., η n ), with elements η i = x T i β is introduced, which is linked to the mean vector µ = (µ 1 , µ 2 , ..., µ n ) via a link function g, as expressed by This formulation allows for the incorporation of multiple predictors and the estimation of their effects on the response variable.The choice of link function g depends on the distribution of the response variable y and can vary between models.For instance, when the response variable is binary (y i = 0 or y i = 1), the logit link function is frequently employed, connecting the mean of the response µ variable to the logarithm of the odds ratio.This can mathematically be expressed as: Similarly, if the response variable is a non-negative variable, the link function can be logarithmic, relating the mean of the response variable to the linear predictor via After selecting the appropriate link function, the GLM implies that the linear predictor can be represented as a linear combination of the predictor variables where β 0 is called the intercept, and β 1 , β 2 , . . ., β k are the regression coefficients.The estimation of these coefficients will be performed using the method of maximum likelihood estimation, detailed in Section 3.2.2.This method quantifies the regression coefficients that are most probable given the observed data (y i , x i ) n i=1 , considering the assumed conditional distribution of the response variable y and the selected link function g.

Application of Generalized Linear Models to HVOF Data
In the HVOF process, the response variables of interest include the coating properties, such as roughness, porosity, layer thickness, and hardness, as well as the in-flight properties, such as particle temperature and particle velocity (see Table 2.1).Since these response variables are continuous and positive, it is necessary to choose an appropriate probability distribution and link function for modelling them.The gamma distribution with a log link function is a common choice for non-negative continuous data [10] and thus, we will employ it in our analysis.

The Gamma Distribution
A continuous, non-negative random variable Y is said to follow a gamma distribution with shape parameter a > 0 and rate parameter b > 0, denoted as Y ∼ G(a, b), if it has the density function: y a−1 exp(−by) , y > 0 .
The expected value and variance are given by E(Y ) = a b and V(Y ) = a b 2 .An illustrative comparison of gamma distributions with varying shape and rate parameters is provided in Figure 3.1.Occasionally, the gamma distribution is defined via an alternative parameterization.Depending on the expected value µ and the scale parameter ν > 0, the density is then given by: where µ = E(Y ) is the parameter of interest and the variance ν = V(Y ) is considered as a nuisance parameter, meaning that the value of ν is not the main focus of the analysis.
In other words, while ν plays a role in determining the shape of the density function, it is not the parameter that one aims to estimate or draw conclusions about.

Maximum Likelihood Estimation Gamma Regression
The likelihood function is a fundamental concept in statistical inference that quantifies the plausibility of the observed data under a given statistical model.As a consequence of the conditional independence of y i given x i , it is defined as the product of the probability density function of each observation in the sample, conditioned on the parameter values.In other words, the likelihood is the joint probability of the observed data, viewed as a function of the parameters.Therefore, the product of the likelihood contributions, defined in 3.4, yields the likelihood for the observed data, providing a basis for inference on the unknown parameters.However, it is essential to acknowledge that the assumption of conditional independence of the response variable y may not always hold in practical applications, particularly in complex systems or processes where various factors can influence the response variables.
In the context of this study on the coating process, assuming conditional independence implies consistent equipment-and process performance across observations, with response variable variation attributed solely to considered explanatory variables.
Nonetheless, factors like equipment stability, environmental conditions, or procedural variations could introduce dependencies among response variables, challenging this assumption.To address this concern, careful execution of the experiments was performed, encompassing thorough validation of equipment functionality and accurate examination for potential issues, such as instrument cleaning and procedural consistency.The impact of changing environmental conditions can be disregarded due to the operations being conducted within a coating booth equipped with a continuous suction system.These measures were taken to validate the conditional independence assumption and ensure the reliability of the experimental results.
An empirical investigation of the dataset considered in this study reveals a notable right-skewness across all examined response variables.While marginal distributional properties alone do not necessarily dictate the choice of distributional family for modeling conditional means, this observation suggests a departure from symmetric distributional patterns.Furthermore, these variables exhibit non-negativity and continuity.To account for these specific distributional characteristics, the assumption of a gamma regression framework is made (cf. Figure 3.1).In addition to the assumed gamma distribution, other distributions such as the log-normal distribution or the inverse Gaussian distribution can also be considered in this setting.Since similar results are anticipated for these alternative distributions, the emphasis on the specific distribution assumption is relaxed.Therefore, the assumption made in this study is that the response variables are conditionally gamma-distributed, with the expected value depending on the explanatory variables.Future investigations are warranted to explore the implications of alternative distributional assumptions in this domain.
The gamma regression model with a logarithmic link function assumes that the response variable y i for i = 1, . . ., n follows a Gamma distribution with mean µ i and scale parameter ν > 0. The mean µ i is modeled as a function of the covariates x i = (x i 1 , x i 2 , . . ., x i k ) through the logarithmic transformation, which is defined as , where η i is the linear predictor.For simplicity, we consider a gamma regression model with a single covariate x.Nevertheless, it is important to emphasize that extending the model to multiple covariates is possible and adheres to the same theoretical framework outlined here.In particular, the incorporation of additional predictors would require a simple augmentation of the linear predictor η i to account for their effects.Therefore, the model can be readily extended to encompass more intricate predictor configurations, as warranted by the research question at hand.The univariate regression model used in this work, in which a single dependent variable is considered, can be summarized in the form of the following equations: The characterization of HVOF coating quality involves multiple aspects, prompting consideration of a multivariate modeling approach.Such an approach allows for the es-timation of covariate effects and the assessment of a variance-covariance matrix, which quantifies correlations among different quality criteria and potentially enhances predictive capabilities.However, due to practical limitations such as the limited sample size and the need for additional data for variance-covariance matrix estimation, a univariate approach was chosen.Moreover, within the specific applications where the coated material is applied, only a subset of properties listed in Table 2.1 is relevant, thus a univariate approach was considered more appropriate.
Given the univariate model in (3.3) and assuming conditional independence of y i |x i it is now possible to define the likelihood function.
Definition 3.1.The likelihood function L(β|y) of the observed data for the model described in (3.3) is defined as the product of the likelihood contributions L i (β|y i ), i.e.,
In the process of maximizing the log-likelihood function, the score function serves as an essential mathematical tool.It accurately measures the sensitivity of the loglikelihood to changes in the parameters of interest.Furthermore, the maximum likelihood estimator (MLE) β is defined as the solution of s( β|y) = 0 .The score function s(β|y) quantifies the rate of change of the log-likelihood function as the parameter values are varied, and provides a measure of the direction and magnitude of the parameter updates that increase the log-likelihood.It can be computed by either numerically or analytically differentiating ℓ(β|y), which in our case leads to Proposition 3.1.Let s(β|y) represent the score function for a response vector y, as defined in Definition 3.2, and let y consist of observed values y i from a random variable Y i following a gamma distribution, as described in Here, X = (1 x) represents the design matrix, consisting of the explanatory variable x = (x 1 , . . ., x n ) T , where x i ∈ R, y = (y 1 , . . ., y n ) T is the response vector with y i ∈ R + , µ(β) = (µ 1 (β), . . ., µ n (β)) T is the mean vector with µ i (β) ∈ R + , given in (3.3), and 1 = (1, . . ., 1) T is a vector of ones.
Proof.This proof is adapted from [10], with suitable changes accounting for the gamma regression framework considered here.First of all, the first partial derivatives of individual log-likelihoods log L i (β|y i ) are given by Together with the definitions of the vectors x, y, µ(β), and 1, as well as of the design matrix X = (1 x) and the definition of the score function s(β|y) there holds which completes the proof.
Remark 3.1.Note that the specific form of the design matrix X, as described here, applies only to the model under consideration.In general, the design matrix X comprises not only the original explanatory variables but also their higher-order powers and/or products.This expanded form allows for the estimation of higher-order effects or interaction effects between two or more covariates.[10] By setting the score function s(β|y) to zero, a linear system of equations for (β 0 , β 1 ) arises that needs to be solved numerically.The numerical algorithm used in this work (cf.Section 3.2.3)involves the computation of the observed information matrix H(β|y) (= Hessian matrix) or expected information matrix F (β|y) (= Fisher matrix), which is a key component of the algorithm.Note that setting the score function s(β|y) to zero is independent of ν, meaning that the process of finding solutions for (β 0 , β 1 ) is not influenced by the value of ν.While solving the equation s(β|y) = 0 does not inherently guarantee a maximum, the concave nature of the log-likelihood function ℓ(β|y) in this model ensures its maximization [35].This concavity can be confirmed by verifying the positive semi-definiteness of H(β|y), defined in 3.5.For further insights into the existence and uniqueness of the maximum likelihood estimator in generalized linear models, refer to [35].
In general, there is no guarantee that solving the equation function s(β|y) = 0 yields a maximum.However, since the log-likelihood ℓ(β|y) for this model is a concave function, solving the equation will yield a maximum.The concativity for ℓ(β|y), defined in 3.5, can be verified by checking if H(β|y) is positive definite.The expected information matrix F (β|y) is defined as where E[•] denotes the expected value.
The matrices H(β|y) and F (β|y) quantify the amount of information that the observed data provides about the unknown parameters of the model.For our specific gamma regression framework, they can be computed explicitly as described in the following Proposition 3.2.Let y be as in Proposition 3.1, let W = diag(ν y i /µ i (β)) i=1,...,n be a diagonal matrix with elements ν y i /µ i (β) and W = diag(ν) be a diagonal matrix with elements ν.Then the observed information matrix H(β|y) and the expected information matrix F (β|y), defined in (3.5) and (3.6), respectively, can be expressed as H(β|y) = X T WX , and F (β|y) = X T WX .
Proof.The second partial derivatives of individual log-likelihoods log L i (β|y) are given by (3.7) The observed information matrix H(β|y) is obtained through the aggregation of the second partial derivatives of individual log-likelihoods log L i (β|y), i.e.,

H(β|y
. Together with the definition of W we thus obtain H(β|y) = X T WX .
Since E(y i ) = µ i (β), the Fisher matrix F (β|y) is given by which yields the assertion.

Numerical Computation of the Maximum Likelihood Estimator
Numerical algorithms are essential for estimating the maximum likelihood estimator of parameters in a statistical model, particularly when an analytical solution to the likelihood equations is unattainable [13].Frequently, the likelihood function is an intricate, nonlinear function of parameters, lacking a closed-form expression for its maximum, e.g., in gamma regression with a logarithmic link function.In such cases, numerical algorithms such as the Newton-Raphson algorithm are employed to iteratively approximate the solution of the likelihood equations until convergence is reached [10].These methods rely on numerical techniques to estimate the derivatives of the likelihood function, which are used in computing the updates to the parameter estimates.
Newton-Raphson Method [13] is an iterative method used to find a value of β that satisfies the equation s(β|y) = 0, which corresponds to the point where the log-likelihood function is maximized.The Newton-Raphson algorithm achieves this by iteratively approximating the solution of s(β|y) = 0 using Taylor series expansion of s(β|y) around the current estimate of β.Specifically, the expansion can be written as: where β (k) is the estimate of β at the k-th iteration, s(β (k) |y) is the score function evaluated at β (k) , and H(β (k) |y) = −∂ s(β (k) |y)/∂ β T is the observed information matrix evaluated at β (k) .The score function is approximated using a linear tangent line, resulting in an improved approximate solution.This involves finding the root of the tangent line in (3.8).Thus, the algorithm approximates the maximum likelihood estimator of β by solving the equation: for β, which yields: The algorithm iterates until convergence is achieved, which is typically defined as the point at which the change in the estimate of β between two successive iterations falls below a certain threshold.
Fisher Scoring Method [10] is a useful approach for maximum likelihood estimation that involves replacing the observed information matrix H(β (k) |y) by the expected information matrix F (β (k) |y) in the update formula (3.9), i.e., This simplifies the required computations, making it faster and more efficient.

Asymptotic Properties of the Maximum Likelihood Estimator (MLE)
Given the gamma regression model with logarithmic link function, as defined in (3.3), and the MLE procedure presented in the previous section, we now investigate the asymptotic properties of the MLE of the regression coefficients β = (β 0 , . . ., β k ) T .Specifically, under mild regularity conditions introduced below, the MLE can be proven to be a consistent and asymptotically normal estimator, with its asymptotic covariance matrix being equivalent to the inverse of the Fisher information matrix [9].

Linear Hypothesis Testing
By conducting hypothesis tests on the estimated regression coefficients β, one can provide evidence-based justifications for the inclusion or exclusion of specific predictors, ensure the robustness and reliability of a model, and enhance the interpretability and generalizability of the findings.Testing a linear hypothesis on the coefficients of the Gamma GLM can be represented as follows: where C is a known r ×p matrix of rank r, β is the p×1 vector of regression coefficients, and d is the r×1 vector of known constants.This matrix C is used to define the specific hypothesis being tested, and its structure depends on the research question at hand.In the context of our study, C is constructed to examine the significance of certain predictors in relation to the response variable.
Under hypothesis H 0 , the unrestricted maximum likelihood estimator β is not efficient, and therefore we need to consider restricted estimators that take into account the constraints imposed by H 0 [10].For this, we consider the Wald statistic w given in Definition 3.5.The Wald statistic w is defined as: where X is the n × p design matrix, and W is the n × n diagonal matrix with the weights w i on the diagonal.
In the specific case of predictive modelling in HVOF coating, hypothesis testing plays a vital role in determining the relevance of regression coefficients β j , where β j denotes a subvector of β.Specifically, we consider the case where the null hypothesis H 0 : β j = 0 versus the alternative hypothesis H 1 : β j ̸ = 0. Proposition 3.4.Let β j be a subvector of β with dimension r, d = 0, and C be a r × p matrix with 1 at the entries corresponding to the elements of β j and 0 elsewhere.With this choice the Wald statistic w, defined in (3.11), takes the form where A j is the submatrix of the asymptotic covariance matrix A = (X T W X) † corresponding to the elements of β j .
Proof.Assuming that the prerequisites for β j and d, as specified in Proposition 3.4, are met, and with the matrix C taking on the following form: Here, β j represents the initial r regression coefficients, given by (β 1 , . . ., β r ) T .Together with the definition 3.5 there holds which yields (3.12).
In accordance with Proposition 3.4, the assessment of the relevance of a subvector β j is determined by (3.12).If β j is one-dimensional, the Wald statistic w corresponds to the application of a t-test [10].The test statistic, denoted as t j , quantifies the extent to which the estimated coefficient βj deviates from H 0 , taking into account the corresponding standard error, i.e., with a jj the j-th diagonal element of A = (X T W X) † .According to [10], t j is tdistributed with n − p degrees of freedom and H 0 is rejected at significance level α if Alternatively, one can also perform the Likelihood-Ratio test using the likelihood ratio L statistic, defined as where L( β|y) is the likelihood function for the unrestricted estimator, and L( βHo |y) is the likelihood function for the restricted estimator obtained by maximizing the likelihood subject to H 0 .Analogous to the Wald statistic, L follows an asymptotic χ 2distribution with r degrees of freedom under the null hypothesis H 0 .Linear hypothesis testing serves as a tool to assess the significance of estimated regression coefficients within a specified confidence level.This approach enables the determination of whether a particular predictor variable contributes meaningfully to the model's description or if a simpler model could suffice without sacrificing essential information.In contrast, model selection criteria, described in the next subsection, aim to identify the most suitable model for predicting outcomes accurately.

Remark 3.2. Note on Statistical Power and Variable Selection
In statistical inference, it is crucial to consider two types of errors.Type I error occurs when the null hypothesis is incorrectly rejected, mistakenly identifying an effect or relationship that does not exist.This risk is quantified by the significance level α.Conversely, Type II error arises when one fails to detect a genuine effect, incorrectly retaining the null hypothesis.This does not necessarily mean there is no effect; rather, it may reflect the test's limitations.Decisions to exclude terms from a model based solely on statistical significance should be made cautiously.While simple models are preferred for their ease of interpretation, overly strict criteria for variable selection may lead to important predictors being overlooked.Type II error risk is often denoted by β (distinct from regression parameters).The probability that a statistical test will correctly reject a false null hypothesis is known as the power of the test and is represented by 1 − β.High power increases confidence in hypothesis test outcomes, while low power raises doubts about non-significant findings.Statistical power relies on factors such as the significance level α, the sample size n, and the population effect size (ES) [5].
In the context of predictive modeling for HVOF coating, domain expertise is essential in addressing statistical power challenges.Due to the constraints of a small sample size, compounded by the laborious and expensive nature of experiments (cf.Section 5), the statistical power of hypothesis tests is inherently limited.Consequently, a thorough examination of regression coefficients was made in collaboration with thermal coating technicians to assess the relevance of predictors, particularly in cases where the performed test might not achieve statistical significance or the model selection criterion decides to exclude the respective effect.Further techniques for calculating and enhancing statistical power in regression analysis are explored in [6].

Model Selection Criteria
In practice, it is often necessary to compare different models and select the one which provides the best balance between model fit, reflecting the agreement with the observed data, and model complexity.Various criteria can be used for this purpose, including the Akaike Information Criterion (AIC) [2].The AIC is based on the maximized loglikelihood function ℓ(β|y) and is defined by: AIC := −2ℓ( β|y) + 2p ; (3.13) where β is the maximum likelihood estimate of the model parameters, and p is the number of parameters in the model.The AIC penalizes models with many parameters, thus favoring models that fit the data well but are not too complex.Smaller AIC values indicate better models, with a difference of 2 between AIC values suggesting strong evidence in favor of the model with the lower AIC.However, note that the AIC is a relative measure of model fit and should be used for comparing models within the same class.For example, the AIC cannot be used to compare a gamma regression model to a Poisson regression model, since they belong to different classes.The application of model selection criteria such as the AIC is valuable in predicting HVOF coating properties based on process conditions.While it is important to develop accurate prediction models to optimize coating performance and ensure the desired coating properties, it is worth to consider that including too many irrelevant parameters in the model can introduce disturbances and adversely affect its predictive ability.

Assessing Predictive Performance of HVOF Coating Models
To assess the predictive performance of the HVOF regression model, the commonly employed technique of Leave-One-Out-Cross-Validation (LOOCV) is utilized.It allows for a comprehensive evaluation of the model's generalization ability and its accuracy in forecasting coating properties.LOOCV is particularly suitable for evaluating the model's generalization capability when only a limited number of observations is available [36].The LOOCV approach is a computationally intensive procedure, requiring the model to be fit n times, i.e., once for each observation in the dataset.To improve computational efficiency, alternative resampling techniques such as k-fold cross-validation may be used.
The LOOCV procedure involves iteratively fitting the model using all observations except one, and then using the fitted model to predict the response for the left-out observation.This is repeated for each observation in the dataset, resulting in n predicted responses.The predicted response for the i-th observation is denoted as ŷ(−i) , where the superscript (−i) indicates that the i-th observation was left out during the fitting.
The prediction error for the i-th observation is defined as the difference between the predicted response and the observed response, i.e., ϵ i = y i − ŷ(−i) .

Definition 4.1 ( [14]
).The LOOCV estimate of the expected out-of-sample prediction error, i.e., the expected difference between the model's predictions and the true values of new, unseen observations, is defined by: where n is the number of observations in the dataset.
The LOOCV estimate of the expected out-of-sample prediction error is an unbiased estimator of the true out-of-sample prediction error and can be used to compare the predictive performance of different models.The smaller the value of CV (n) , the better the predictive performance of the model.In addition to the LOOCV, we also use the R 2 statistic, which measures the proportion of variance in the observed response that is explained by the model.Definition 4.2.The R 2 statistic is defined as: where n is the number of observations, y i is the observed response for the i-th observation, ŷi is the predicted response for the i-th observation, and ȳ is the mean of the observed responses.
The R 2 statistic can take values between 0 and 1, with higher values indicating a better fit of the model to the data.However, the R 2 statistic can be biased towards models with more predictors, even if the predictors have little or no effect on the response.To address this issue, the adjusted R 2 statistic is used, which adjusts R 2 for the number of predictors in the model.Definition 4.3.The adjusted R 2 statistic is defined as: where p is the number of predictor variables in the model, n is the number of observations, and R 2 is the statistic, defined in (4.1).
The adjusted R 2 takes into account the trade-off between model complexity and model fit, and provides a more reliable measure of the model's predictive performance, compared to the traditional R 2 , since it also accounts for the number of predictors p.

Application to HVOF Coating: Practical Implementation
The HVOF process is influenced by a multitude of variables, making it challenging to identify the most important factors that actually impact coating properties.In this study, a selection of five factors was deliberately chosen, guided by the knowledge of thermal spray experts who identified these variables as significant determinants influencing the HVOF process.Moreover, a well-designed experiment is crucial to efficiently collect data on the effects of various factors on the process outcomes.The selection of an optimal experimental design is essential within the domain of HVOF coating, primarily attributed to the considerable costs and time-intensive nature associated with conducting experiments using coating materials.Furthermore, a carefully planned experimental design enables strategic allocation of available experiments, maximizing information and providing valuable insights within a limited experimental scope.
In industrial processes, statistical design of experiments (DoE) is considered a reliable technique for conducting experiments.A DoE allows for the systematic variation of process variables (= explanatory variables), which enables a more comprehensive understanding of their impact on the outcome.In contrast to the traditional onefactor-at-a-time approach, where interaction effects between two or more explanatory variables cannot be estimated, a DoE approach enables concise mathematical analysis of the resulting data and facilitates the identification of significant factors and their interactions.In addition, DoE allows researchers to investigate complex relationships between explanatory, revealing hidden insights, and supporting the optimization of industrial processes.

Central Composite Design
The central composite design (CCD), a well-established and commonly employed experimental design in the field of industrial process optimization, is utilized in this work to acquire empirical data for the HVOF process.Compared to other designs, the CCD is particularly useful as it can be efficiently used for fitting second-order models [22], i.e., estimating the effects of factors and their interactions in a quadratic form.The design incorporates a systematic variation of pre-defined factors, employing three levels (−1, 0, 1) for each factor.Additional star points are included to enable the inclusion of quadratic terms in the model [22].Explanatory variables, which are given in quantitative form, are transformed into qualitative factors.The center point x 0 , represented by the level 0 for each factor, serves as a reference point and is used to assess the impact of factors on the system.The cube points correspond to the corners of the experimental region, represented by the levels (−1, 1).The star points are additional experimental points that are used to estimate the behavior beyond the linear response and to identify potential quadratic effects of factors.These points are positioned at a value of α, where α is determined for a explanatory variable x as where x 0 is the center point, δ x the difference between x 0 and the quantitative value that corresponds to −1, and k is the number of explanatory variables under consideration.After transformation from quantitative values into qualitative ones, the values of x 0 , x 0 ± δ x , and ±α are replaced by 0, ±1, and ± √ k respectively.These qualitative values are then used in the design matrix X to represent the explanatory variables (= factors).The number of factors k, which represent the process variables under investigation (cf.Section 5.2.1), directly determines the number of cube and star points in the CCD.Specifically, there are 2 k cube points and 2k star points included in the design.Additionally, the CCD consists of n c center points, where n c represents the total number of (potentially repeated) center points.To enhance the efficiency and accuracy of the design, a spherical CCD was employed with the choice of α = √ k concerning the star points.The spherical design allows for the estimation of effects of any factor with equal precision and reduces the risk of overemphasis on any factor.Thus, an optimal balance between precision and stability of the model parameters is obtained, which is important for receiving reliable estimates of the factor effects and their interactions.As recommended in [22], it is essential to randomize the experimental runs to avoid the influence of uncontrolled sources of variation.

Identification of Influencing Factors
Based on a review of the literature [27,31,38], previous one-factor-at-a-time experiments, and expert knowledge by thermal spray coating experts of voestalpine TSM [1], five key factors, which are described in the next paragraph, were identified for systematic variation: powder feed rate (PFR), stand off distance (SOD), lambda (λ), i.e. the stoichiometric ratio of oxygen to fuel, coating velocity (CV), and total gas flow (TGF).The schematic diagram in Figure 5.2 provides a comprehensive visual representation of the considered key process factors.Using these k = 5 factors and conducting n c = 7 replications at the central point, a total of 49 trials were carried out, forming the CCD.The experiments were conducted using a rotational setup that included a turning lathe, allowing for the application of the thermal spray coatings (cf. Figure 5.3).The selected factors play a critical role in the HVOF coating process, exerting significant influence on the quality and performance of the resultant coatings.The PFR governs the amount of coating material supplied, while the SOD regulates the spacing between the spray gun and the substrate.The stoichiometric ratio of oxygen to fuel (λ) ensures specific combustion conditions.Furthermore, the CV, determined by the combined influence of the robot traverse speed and the rotational speed of the turning lathe (cf. Figure 5.2), enables precise control over the deposition process.Finally, the TGF is constituted by the summed gas flow of fuel, oxygen, and air, collectively governing the overall flow rate of the combustion gases.
Each of the five factors is accompanied by a designated set of predefined levels of variation, which are listed in Table 5.1.These levels were determined to cover a range of values that would effectively capture the variability and impact of these factors on the desired coating properties.The chosen levels allow for a systematic and comprehensive exploration of the parameter space.

Factors
Coded values The HVOF coatings were produced using an Oerlikon Metco thermal spraying equipment, namely the DJ 2700 gas-fuel HVOF system with water-cooled gun assembly.The fuel gas used for these tests was propane, its amount and ratio defined by the two key factors TGF and Lambda.For the process preparation, steel plates of type 1.4404 were welded onto an axis mounted on a turning lathe for rotational spraying.All samples were degreased with acetone and sandblasted with alumina before thermal spraying.The powder used for the spraying process was an agglomerated sintered tungsten carbide powder (WC-10Co-4Cr) with a grain size in the range of -45+15 µm, supplied by Oerlikon Metco.The photograph presented in Figure 5.3 showcases the experimental setup employed during the HVOF coating process, wherein the dynamic engagement of the robot, turning lathe, and coating stream can be observed.
Additional thermal spraying attributes like cooling, powder feed gas, pressure and number of passes, i.e., number of times the coating material was applied or sprayed onto the substrate during each experimental run, were kept constant throughout the experiments.In addition to the classic coating properties such as roughness, porosity, layer thickness, and surface hardness, the deposition rate, deposition efficiency, and in-flight particle properties such as particle velocity and particle temperature of the powder particles were measured.
Two different in-situ measurements were performed in the course of these trials.On the one hand, the in-situ particle characterization and on the other hand, the in-situ pyrometric temperature measurement.The particle characteristics were measured using a Spraywatch camera with the software SW4 (supplied by Oseir).The temperature of the sample surface was constantly measured using a Keller pyrometer.
The surface roughness of the sprayed samples was measured using a mobil roughness tester Hommel Etamic Waveline W5.Coating hardness was assessed on the surface using a Cisam-Ernst S.r.l E-Computest mobile hardness tester, using a spheroconical diamond at a load of 5 kg and a testing time of 2 seconds.In addition to the surface characterization, cross-sections of each sample were prepared (according to internal preparation procedure WC) to analyse the coating thickness.The coating thickness and the coating porosity were determined using image analysis software, IMS Client, applied to microscopic images captured with a Zeiss Axio Observer.Z1m. Figure 5.4 provides visual evidence of the observed variations in coating thickness, as captured in the microscopic images acquired from the IMS Client software.

Empirical Results of Experiments in HVOF Thermal Spraying
In this section, the empirical findings derived from the comprehensive analysis of the experimental data are presented, demonstrating the effectiveness and utility of the gamma regression approach in analyzing the relationships between key process variables and coating properties.The analysis and modelling were performed using the statistical software R (version 4.2.2).The gamma regression models were implemented using the glm function [28] with the Fisher scoring algorithm for estimating the regression coefficients β (= ML estimates).
The properties listed in Table 2.1 serve as an overview and exemplify various potentially relevant properties of the HVOF process.This study focuses on analyzing 8 properties, which include in-flight characteristics (particle velocity and particle temperature), performance metrics (deposition rate and deposition efficiency), and coating attributes (thickness, roughness, hardness, and porosity).Analysis of the phase composition is typically not relevant for the coating material used.It should be emphasized, however, that the presented methodology holds equal relevance for the analysis of the other response variables in Table 2.1.Table 6.1 provides a detailed illustrative example of the estimated regression coefficients and their corresponding standard errors of the deposition rate model and the deposition efficiency model.Additional tables containing analogous information for the remaining response variables can be found in the Appendix, specifically Table A.1 and A.2.These tables contain two kinds of models for each property, a full and a reduced version.The full model encompasses all predictor variables that can be estimated by utilizing the CCD methodology, while the reduced model is derived through variable selection based on criteria such as the AIC and hypothesis testing for coefficient relevance, as described in Sections 3.2.5 and 3.2.6.The model selection procedure involved the following steps: Initially, the full model was constructed, and non-significant coefficients were iteratively eliminated in a backward direction using the AIC as the guiding metric.If the removal of a coefficient resulted in a reduction in the AIC, the significance of the respective predictor was reassessed through a hypothesis test in consultation with thermal spray technicians.This consultation aimed to assess the practical significance compared to the full model, the metrics introduced in Section 4 are computed for each model individually.Table 6.2 summarises the outcomes of the gamma regression analysis for in-flight properties (velocity and temperature), performance indicators (deposition rate and deposition efficiency), and coating properties (thickness, roughness, hardness, and porosity).Once again, the outcomes of both the full and reduced models are presented, highlighting their ability to model the studied properties.In addition to the number of coefficients N p and the model selection criterion AIC as in the preceding Table 6.1, this table also incorporates important performance metrics, namely R 2 ,R 2 adj , and CV (n) , to measure the predictive quality of the regression models, as described in Section 4.
Concerning the in-flight properties, both the full and reduced models demonstrate favorable results.The full model for particle velocity exhibits a high R indicating a strong fit to the observed data.However, taking into account the number of predictors N p in the model, it is advisable to consider the adjusted R 2 value of 0.89, which accounts for the model's complexity.Conversely, the reduced model for velocity yields a slightly lower R 2 value of 0.93, yet a higher adjusted R 2 value of 0.92 compared to the full model.These findings, coupled with lower values of the Akaike Information Criterion AIC and reduced out-of-sample prediction error CV (n) , suggest that the reduced model offers superior predictive performance.Similar patterns emerge for the regression models investigating the other target variables in Table 6.2, consistently favoring the reduced models.According to the adjusted R 2 , all reduced models demonstrate a good fit to the data, with values exceeding 0.8.However, the model for coating porosity yields less satisfactory results.This observation may be attributed to the volatile nature of the porosity measurement technique using Image Analysis, as discussed in [3] In addition to the findings in Table 6.2, Figure 6.1 provides a visual representation of the deposition rate predictions obtained from the full and reduced models.Each data point on the scatter plot represents an experimental trial, where the y-axis corresponds to the observed values, and the x-axis represents the LOOCV predictions (refer to Section 4).The color-coded points differentiate between the center points, cube points, and star points obtained from the CCD.
Upon analysis of the scatter plot, it is evident that the reduced model yields pre-dictions that are closer to the diagonal line, indicating a higher degree of agreement between the predicted and observed values.This closer alignment implies a more accurate prediction of the deposition rate by the reduced model compared to the full model.Moreover, as expected, the prediction accuracy varies across different design points.The star points (yellow) demonstrate relatively lower predictive accuracy compared to the center points (blue) and the cube points (green), although this discrepancy is observed only in a subset of star points.This outcome underscores the challenges associated with extrapolating the model's behavior to regions outside the training data, emphasizing the need for caution when interpreting predictions for such points.Overall, Figure 6.1 provides strong evidence supporting the superior predictive performance of the reduced model in estimating the deposition rate.The analysis of these visual results further strengthens the findings presented in Table 6.2, reinforcing the advantages of employing the reduced model in understanding and predicting the deposition rate more accurately.Furthermore, similar findings regarding the superior predictive performance of the reduced model are also observed for the other analyzed target variables.

Conclusion
This study proposed a framework for modelling and predicting critical target variables in HVOF coating processes.By utilizing DoE and GLMs, accurate estimation of model parameters was achieved through maximum likelihood estimation.The framework incorporated a careful selection of predictor variables based on their significance and contribution to the coating properties, enhancing the model's interpretability and predictive performance.The application of this framework to experimental data from thermal spray coating experiments demonstrated its effectiveness in predicting target variables and providing insights into the relationships between factors and coating properties.The systematic variable selection process helps identify the most influential factors and eliminates irrelevant or redundant variables, simplifying the modelling process and improving the accuracy of predictions.The proposed framework has the potential to optimize thermal spray coating processes and contribute to the development of more efficient coating technologies in various industries.By developing a comprehensive understanding of the intricate interplay among process variables, material properties, and coating microstructure, manufacturers can enhance the functionality and performance of coated surfaces.This, in turn, can lead to improved product quality, extended component lifespan, and reduced maintenance costs.
In future investigations, we intend to expand our framework by introducing additional variables and exploring their interactions.This includes varying previously held constant factors, such as process gas pressures, across different levels to assess their impact.Furthermore, we plan to compare gamma regression models with alternative statistical models that require distinct distributional assumptions.While our previous work revealed the effectiveness of machine learning algorithms, particularly support vector machines, in predicting HVOF-related properties [29], we aim to explore a hybrid approach combining these methodologies to further enhance predictive accuracy.Given the dynamic nature of the HVOF process, where process variables often deviate from target values, sensors will be installed within the booth to monitor these variations.Using advanced modeling techniques, this sensor data together with additional quantitative features are used to improve predictive capabilities.Further experiments will be conducted to ensure that a sufficient number of samples is available.Following satisfactory performance in predicting coating quality properties related to WCCoCr, the framework will be extended to other coating materials and their associated characteristics.
Overall, this study provides a systematic and data-driven approach to modeling and predicting coating properties in thermal spray coating.By leveraging this framework, researchers and practitioners can advance the understanding and optimization of thermal spray processes, leading to advancements in surface technology and its applications across industries.The variable selection process improves prediction accuracy and facilitates informed decision-making in the coating optimization process, contributing to the overall improvement of coating methodologies.

A Supplementary Tables of Estimated Regression
Coefficients and Standard Errors

Figure 2 .
Figure 2.1: Schematic depiction of the High Velocity Oxygen Fuel (HVOF) process.

Figure 3 . 1 :
Figure 3.1: Comparison of gamma distributions with varied parameters.The left panel displays the probability density function (PDF) of a gamma-distributed random variable y with a fixed shape parameter a = 5 and varying rate parameters b.The right panel illustrates the PDF of a gamma-distributed random variable y with a fixed rate parameter b = 2 and varying shape parameters a.

Assumption 3 .Definition 3 . 4 .
1 ([9] Regularity Assumptions).Let β ∈ B ⊂ R p denote the ML estimator for the true parameter β, p be the number of predictor variables in the model, and M the image µ(β) of β.Furthermore, the linear combination of the predictor variables η is related to the mean µ(β) of the response y by an injective link function g : M → R p , i.e., η = g(µ(β)) (compare with (3.3)).Additionally, there holds(i) B is open in R p , (ii) The design matrix X has full rank, i.e., rank(X) = p, (iii) g(•) is twice continuously differentiable on M.Note that Assumption 3.1 is valid for our gamma regression model with a logarithmic link function (3.1), i.e., where the response variable follows a gamma distribution(3.3).An estimator β is consistent for the true parameter vector β if, as the sample size n goes to infinity, β converges in probability to β.In other words, for any small positive number ϵ, it holds thatlim n→∞ P (|| β − β|| > ϵ) = 0 .Using the Law of Large Numbers, the sample mean of a sequence of i.i.d.random variables with finite mean converges in probability to the expected value.Since the loglikelihood function ℓ(β|y) in this model(3.3)  is the sum of i.i.d.Gamma distributions, the Law of Large Numbers can be used to establish convergence in probability of the MLE to the true parameter values.In the following proposition, two key properties of the gamma regression model are established without providing formal proof.Proposition 3.3 ( [9]).
(i) In the setting of the gamma regression model (3.3), the MLE β is consistent for β.In particular, under the regularity conditions stated in Assumption 3.1, the ML estimator β converges in probability to the true regression coefficients β for increasing sample size, i.e., β p → β, where p → denotes convergence in probability.(ii) Let the assumptions of Proposition 3.3 hold.Then the gamma regression model defined in (3.3) is asymptotically normal in relation to the maximum likelihood estimator (MLE) β, i.e., √ n( β − β) d → N (0, F † (β|y)), where d → denotes convergence in distribution, and n denotes the sample size.

Figure 5 .
Figure 5.1: Central Composite Design with a) k = 2 factors and b) k = 3 factors.

Figure 5 . 2 :
Figure 5.2: Illustration of the considered key factors in the HVOF coating process.

Figure 5 . 3 :
Figure 5.3: Photograph illustrating the experimental setup during the HVOF coating process, showing the robot, turning lathe, and coating stream in action.

Figure 5 . 4 :
Figure 5.4: Microscopic images obtained from the IMS Client software, showing the observed variations in coating thickness across the cross-sectional profiles of the sprayed samples.

Figure 6 . 1 :
Figure 6.1: Scatter plot showing the comparison between observed values and LOOCV predictions for deposition rate using the full and reduced models.The data points are color-coded based on the corresponding center points (blue), cube points (green), and star points (yellow) from the Central Composite Design.

Table 2
.1 provides a comprehensive overview of the various HVOF process variables and coating characteristics discussed above.

Table 2 .
1: Examples of HVOF process variables and characteristics.

Table 5 .
1: Levels of key factors for HVOF coating experiments depicted in Figure 5.1.

Table 6 .
2value of 0.94, 2: Results of the gamma regression analysis for in-flight properties (velocity and temperature), performance indicators (deposition rate and deposition efficiency), and coating properties (thickness, roughness, hardness, and porosity) using full and reduced models.