Linear regression serves as the fundamental starting point for regression methods and remains a valuable and widely used statistical method. In addition, it acts as a solid foundation for exploring newer approaches. Therefore, the significance of a thorough understanding of linear regression before diving into more complex statistical methods cannot be overstated. Mixed-effects models are employed to describe the relationship between a response variable and one or more covariates in grouped data, structured according to factors such as longitudinal observations, repeated measures, hierarchical organization, or block designs [1]. These models extend linear models by incorporating random effects, which introduce an additional error term to account for the correlation between observations within the same group. Mixed models demonstrate broader applicability and greater generality than fixed or random models, making them particularly suitable for analyzing complex data structures with multiple sources of variability. Mixed models can be orthogonal (e.g. [2]) or non-orthogonal (e.g. [3]). Orthogonal models occur when the fixed and random effects are independent of each other, which simplifies the estimation of parameters. Non-orthogonal mixed models arise when there is a correlation between the fixed and random effects.

Moreover, biadditive regression models, which extend the analysis of variance (ANOVA) to models with quadratic terms, have frequent applications in ecology, where the experimental units within a block or stratum are considered a random sample from a population of units and the blocks or strata themselves are viewed as a random sample drawn from a population of blocks or strata [4].

Biadditive regression models are a flexible statistical framework designed to account for both fixed and random effects. The model is given by the expression

where Y is a vector of N random variables Y₁, …, Y_N, Xβ expresses the fixed effects, and ∑i=1wXiZi represents the sum of w independent random terms, each associated with a specific source of variability. The covariance matrix of Y is structured additively as

where Mi=XiCiXit, and σi2 and C_i, i = 1, …, w, unknown and known and invertible, respectively. Random terms Z₁, …, Z_w were assumed to have independent and identically distributed (i.i.d.) components with null mean values and higher-order cumulants (_cr)1, …, (_cr)w, r = 1, 2, 3, 4. However, these conditions have been refined. Currently, Z_i1, …, Z_iw are treated as independent, with null mean vectors and covariance matrices σ2C1,…,σ2Cw. This revised assumption allows for a more flexible and realistic modeling of variability in the data, accommodating non-identical covariance structures across random effects.

Alexandre et al. [5] conducted a detailed analysis to estimate the covariance components σ12,…,σw2 and the coefficient vectors in the biadditive regression models. This estimation process involved modeling the variability and dependence structure of the data through the covariance terms, which represent the scale of variation for each random effect Z_i. Additionally, the coefficient vectors β were estimated to quantify the contributions of fixed effects, capturing the systematic relationships between covariates and the response variable.

In this paper, we present biadditive models with two extended frameworks: the orthogonal block structure and the commutative orthogonal block structure. These models introduce a novel perspective, enabling a more detailed exploration of underlying covariance structures. We establish commutativity conditions and derive optimal estimators, focusing specifically on best linear unbiased estimators (BLUE) and best quadratic unbiased estimators (BQUE). Furthermore, we extend the concept of biadditive models to encompass families of biadditive models, offering additional possibilities and broadening the applicability of this approach.

This paper is structured as follows. In Section 2, we present biadditive regression models. In Section 3, estimation and inference procedures for two extensions of the biadditive model; OBS and COBS, with controlled heteroscedasticity, including parameter estimation, variance estimation, and unbiased estimation for the model's parameters, are presented. Section 4 presents families of biadditive regression models with OBS and with COBS and results for estimable functions and derives chi-square tests. In Section 5, we developed a simulation study to compare the estimates of variance components, obtained through the proposed approach with those derived from analysis of variance and Markov Chain Monte Carlo methods. Finally, in Section 6, we conclude the paper with some final remarks.

Let us consider a linear mixed model given by

where Z₁, …, Z_w, i = 1, …, w are independent random vectors with covariance matrices σ2C1,…,σ2Cw. If the vectors Z₁, …, Z_w have mean vectors μ₁, …, μ_w, we can introduce the centered vectors Z.i=Zi-μi for i = 1, …, w and the extended coefficients vector β.=[βt,μ1t,…,μwt]t. This approach simplifies the treatment by assuming that the centered vectors Z.i1,…,Z.iw have null mean vectors. Currently, going into inference for Y, given the independence of Z₁, …, Z_w we consider its covariance matrix

where the matrices Mi=XiCiXit for i = 1, …, w provide profound insights into the relationships between variables within our dataset. This analysis reveals crucial information about how each component contributes to the overall variance and covariance observed in the data.

Moreover, we consider the orthogonal basis {α₁, …, αn.} for Ω^⊥ = ℝ(X)^⊥, the orthogonal complement of the range space, ℝ(X) of matrix X. Currently, the vectors with components

have null mean vectors and variance given by

where hℓi=αitMiαℓ, ℓ = 1, …, n., i = 1, …, w, thus expressing the transformation's impact on the covariance structure.

Putting Y.[2]=(Y.12,…,Y.n.2), we obtain

where H = [h_ℓi], ℓ = 1, …, n., i = 1, …, w, and σ2=(σ12,…,σw2)t. The least-square estimator (LSE) for the variance components is vector σ~2, given by

where the symbol (·)⁺ indicates the Moore–Penrose inverse matrix [6].

Similarly, as the expected value of Y is Xβ., the LSE for β. will be

We also obtain the estimators for Σ(Y) given by

and the generalized least square estimator (GLSE) for β.~,

as shown in Kariya and Kurata [7].

This method facilitates the estimation of variance components in mixed models where the random-effects factors may follow various distributions, including non-normal ones. This flexibility is achieved by focusing on the structure of the covariance matrix, rather than imposing a specific distributional assumption for random effects.

Let us now consider two classes of models with orthogonal properties, the OBS and COBS models. Our approach is general case for these models.

We consider families of models sharing the matrices X, X₁, …, X_w, so

where the vectors Z_i(h), i = 1, …, w, h = 1, …, d, are independent and have null mean vectors. We also have,

and the random vectors are

The components Z_i1(h), …, Z_ici(h) are i.i.d with cumulants c_ri(h), i = 1, …, w, h = 1, …, d, r = 2, 3, 4, the same for all models. In addition, the matrices M = XX^t, and

are the same for all models. In the homogeneous case, in which the matrices Mi=XiXit, i=1,…,w, are null, we have the GLSE given by

where

and

A simulation study was conducted to assess the performance of the proposed estimation method. The R programming language was used to generate the simulation data, following the procedure outlined below. The process was repeated a total of N = 1, 000 times to ensure robust and statistically reliable results. In each iteration, random values for the model parameters were generated according to the specified distributions for the random effects and fixed effects. The corresponding observation vectors were then calculated using the model equation. For each simulated dataset, the variance components were estimated and performance metrics such as bias, standard deviation (SD), and efficiency of the estimators were calculated. This repetition allowed for a comprehensive evaluation of the accuracy and precision of the method across a variety of random configurations. Simulate the observation vectors

where X₂ and β_2j represent an additional design matrix and random effects term, respectively.

Random effects were generated according to the following distributions:

where σ1j2 is the variance component for the first random effect,

where a = j and b = 10−j.

The true variance components (σj2) for j = 1, …, 10 were estimated using LSE (Equation 2), where

We estimate variance components and evaluate bias, standard deviation (SD), and δ. to evaluate the performance of methods, analysis of variance (ANOVA), and Markov chain Monte Carlo (MCMC) (Table 1).

For the cases j = 1, …, 10, our estimator consistently exhibited the smallest δ. The probability of this occurring by chance, assuming no superior precision among the methods, would be q=(13)10. This provides strong evidence that our method demonstrates significantly greater precision.

In this paper, we consider biadditive models, often used in studies of manuring and other agronomic applications. We incorporated two extensions of these models: orthogonal block structures and commutative orthogonal block structures, which allow for a more detailed analysis of the additive structure of covariance matrices. In addition to individual models, we considered families of biadditive regression models. By incorporating commutative conditions, we derived optimal estimators, including best linear unbiased estimators and best quadratic unbiased estimators. The proposed methodology extends classical results, such as the Hsu theorem, while providing a robust framework for hypothesis testing within families of models. This framework also emphasizes the estimation of covariance components and coefficients, offering researchers valuable tools for investigating variability across models. The chi-square tests presented here establish a solid statistical foundation for evaluating variability and ensuring precision in agronomic studies. Furthermore, we highlighted the relevance of commutative orthogonal structures for factorial models, particularly those based on prime basis factorials. Such models, often used in studies of manuring and other agronomic applications, showcase the versatility of our approach. To evaluate the performance of the proposed methods, we conducted a simulation study. In this study, we simulated observation vectors based on biadditive regression models with predefined covariance structures and random effects. Using simulations conducted for our model, as well as for the ANOVA and MCMC methods, we estimated variance components and computed performance metrics, such as bias, SD, and δ. Simulation results demonstrated the enhanced precision of our estimation approach.