Small area estimation using reduced rank regression models

Abstract

Small area estimation techniques have got a lot of attention during the last decades due to their important applications in survey studies. Mixed linear models and reduced rank regression analysis are jointly used when considering small area estimation. Estimates of parameters are presented as well as prediction of random effects and unobserved area measurements.

Keywords:

• Growth curve model
• mixed linear model
• reduced rank regression

1. Introduction

Small area estimation is an active research area in statistics with many applications. Today there exist several approaches to deal with small area problems. For recent contributions to the topic when considering discrete data (e.g. see Chandra et al. 2018), and for continuous data (see Baldermann 2018). In both articles many works of other authors are presented and discussed. For those who want to get a solid introduction to the subject the book (Rao and Molina 2015) can be recommended. Usually a common basis for small area estimation problems is that a survey study (finite population case) has been conducted and based on the survey one intends to extract information concerning small domains by adding specific local information which was not accounted for in the comprehensive survey.

In this article a case is considered where units of a survey sample are investigated several times, leading to the existence of dependent observations. It is assumed that observed background information exist which then is implemented via a linear model, giving a model which includes variance components. Moreover, it is assumed that latent processes exist which often show up when different types of systems are studied, for example when weather impact is studied. In this case often a huge number of variables are studied but we do not know the effects of the individual variables, for example we cannot directly relate temperature to plant growth. This leads us to exploit reduced rank regression models. Thus, the model which will be discussed in this article will include time dependent fix effects, i.e. structured mean responses, random parameters and rank restrictions on parameters, which to our knowledge is a new type of model. In Section 2 more motivating details are given, including basic references, and in Section 3 the assumptions are summarized together with all technical details. Briefly, in Section 4 the estimation of all unknown parameters is considered and finally two prediction results are presented in Section 5.

2. Background

In this section, via an example, we will introduce our way of thinking about extracting information from a survey sample. In particular, a number of different sources of uncertainty are presented which then will be included in a statistical model. It turns out that the model will comprise three different types of uncertainty. Section 3 includes all details of the model and in the subsequent section likelihood-based estimators are derived. In the last section these estimators are utelized when finding predictors.

Example 2.1.

Let there be a national survey which collects data about some specific production from a certain type of companies. The survey is followed-up once per year and it is repeated, say, five years.

Altogether there exist 20 regions and each region consists of 10 subregions. The survey sample (units) consists of four companies from each subregion. Thus, in the survey, in total 800 companies are followed five years. With the models presented in this article the aim is to obtain information about the companies in the subregions. Note that in reality the companies should be classified into strata and then at least one company should be sampled from each stratum and each subregion. To simplify the presentation we will only have one stratum. However, in each subregion we assume that there is a lot of background information which we would like to take into account when predicting production in subregions or predicting production of unobserved units.

Small area models, also called small domain models, can be classified into two categories. One category is the area-level models category (Fay and Herriot 1979) and the other category is the unit-level models category (Battese et al. 1988). Both model categories are discussed in detail in a well-written book (Rao and Molina 2015). In this article we are thinking of area levels but if unit-level information would be available the proposed model can also be used.

Example 2.1 indicates the following use of notation. Let $${\widehat{\theta }}_{\mathit{ij}}, i=1,,\dots ,m, j=1,2,\dots ,n_i,\hspace{1em}n=\sum _{i=1}^m{n}_i$$ be the direct estimates from the survey for the jth subarea within the ith area. The estimates are based on the survey design where units have been sampled on the subarea level. Note that we will not discuss the underlying survey in this article, only use the outcome of the survey, i.e. the direct estimators of the subareas.

Statistics is partly about quantifying uncertainty. Often this is carried out by assuming that an estimate is a representation of an estimator which follows some distribution. However, we also have cases where uncertainty is not described via distributions. Examples are rounding errors, specific types of truncations, identification of influential observations, outlier detection, model comparisons, the number of observations in a sample, etc. One choice of distribution for the estimator ${\widehat{\theta }}_{\mathit{ij}}$ would be the one generated by the survey design, i.e. how many and under what model the units have been sampled would determine the distribution. However, for our purposes we have difficulties to work with the probability space generated by the survey design. Instead, we now think of the direct estimates ${\widehat{\theta }}_{\mathit{ij}}$ being misspecified quantities since the survey sample distribution does not take into account those covariables we are going to use in the adjustment of ${\widehat{\theta }}_{\mathit{ij}}.$ The covariables can be used to correct for bias as well as diminish the uncertainty in the estimate. For example, returning to our example, if the survey design does not take into account the size of companies small companies can be overrepresented in the sample and then θ_ij is wrongly estimated.

In order to derive an applicable model lets start with (2.1) $$y_{\mathit{ij}}={\theta }_{\mathit{ij}}+{ϵ}_{\mathit{ij}}$$(2.1) where for notational convenience ${\widehat{\theta }}_{\mathit{ij}}$ has been replaced by y_ij. The joint distribution for $\left\{{ϵ}_{\mathit{ij}}\right\}$ is specified via (2.2) $$ϵ\sim N_n(\mathbf{0},\mathit{V})$$(2.2) where $ϵ=({ϵ}_{11},\dots ,{ϵ}_{1,n_1},\dots ,{ϵ}_{m,n_m})$ is a row vector, V is supposed to be known and determined by the survey design. Often V is solely a function of the number of sampled units. Moreover, $ϵ$ is defined on another probability space than the one generated by the survey design. However, the uncertainty in the direct estimates, due to the survey, is incorporated in the dispersion of $ϵ.$

Furthermore, one basic assumption in this article is that some of the θ_ij will be the same, i.e. the survey estimators are the same for clusters of subregions. To make such an assumption is often realistic but has to be checked in practise. Under this assumption, instead of (2.1), the model can be written in matrix notation (2.3) $$\mathit{y}=\mathit{\beta }\prime \mathit{C}+ϵ$$(2.3) where $\mathit{y}=(y_{11},\dots ,y_{1,n_1},\dots ,y_{m,n_m}),({\theta }_{11},\dots ,{\theta }_{1,n_1},\dots ,{\theta }_{m,n_m})=\mathit{\beta }\prime \mathit{C},$ and $\mathit{\beta }:$ $k\times 1$ consists of unknown parameters and C: k × n is a usual design matrix describing the relations among the elements $\left\{{\theta }_{\mathit{ij}}\right\}.$ Usually C consists of blocks of one-vectors but in principle it can be arbitrary as long as it fits together with the covariate structures which will be presented later. For example, with three areas and four subareas in each area C equals $$\mathit{C}=\left(\begin{array}{cccccccccccc}1& 1& 1& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 1& 1& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 1& 1& 1\end{array}\right)$$

The error in Equation (2.3)(2.3) $$\mathit{y}=\mathit{\beta }\prime \mathit{C}+ϵ$$(2.3) is defined by Equation (2.2)(2.2) $$ϵ\sim N_n(\mathbf{0},\mathit{V})$$(2.2) . The model means that a probability space is only connected to the modeling of $\left\{{\widehat{\theta }}_{\mathit{ij}}\right\},$ i.e. the set of estimated θ_ij, but we include in the model extra variation due to the survey design.

In the introductory example it was stated that the survey was repeated five times. When survey studies are repeated it will for this article be assumed that when there exists a linear trend it is of the same form for all subareas. In this case there will be a matrix Y where each row in the matrix follows Equation (2.3)(2.3) $$\mathit{y}=\mathit{\beta }\prime \mathit{C}+ϵ$$(2.3) , besides an unknown scalar multiplier. Thus, Equation (2.3)(2.3) $$\mathit{y}=\mathit{\beta }\prime \mathit{C}+ϵ$$(2.3) implies (2.4) $$\mathit{Y}=\mathbf{\text{ABC}}+\mathit{E}$$(2.4) where A is a design matrix which models p repeated surveys (within subarea design matrix), B is a matrix of unknown parameters and $\mathit{E}\sim N_{p,n}(\mathbf{0},\mathbf{\Sigma },\mathit{V}),$ i.e. E is matrix normally distributed with a separable dispersion structure $D[\mathit{E}]=D[\text{vec}\mathit{E}]=\mathit{V}\otimes \mathbf{\Sigma },$ where ⊗ stands for the Kronecker product and $\mathbf{\Sigma }$ is an unknown positive definite matrix. The parameter $\mathbf{\Sigma }$ is used to model the dependency due to the repeated measurement which exist because the observed units of the design have been repeatedly measured. No particular structure in $\mathbf{\Sigma }$ is assumed since it would only make sense if we would have precise knowledge about the mean structure, which is not the case in this article. For example, it does not make sense to assume an autocorrelated structure if the mean has not been established. The model in Equation (2.4)(2.4) $$\mathit{Y}=\mathbf{\text{ABC}}+\mathit{E}$$(2.4) is often called Growth Curve model, GMANOVA and recently bilinear regression model (see von Rosen 2018, for basic results and references).

As an example of A we can have $$\mathit{A}=\left(\begin{array}{ccc}1& t_1& t_1^2\\ 1& t_2& t_2^2\\ 1& t_3& t_3^2\\ 1& t_4& t_4^2\\ 1& t_5& t_5^2\end{array}\right)$$ meaning that there is a second order trend which is observed over five years.

When expressing θ_ij as a linear function of an unknown parameter vector $\beta ,$ as it was done in Equation (2.3)(2.3) $$\mathit{y}=\mathit{\beta }\prime \mathit{C}+ϵ$$(2.3) , we usually cannot say that the model is completely true for all subareas. Instead we believe that there is some random variation around a true imaginary model, i.e. we introduce a variance component. Thus, instead of Equation (2.4)(2.4) $$\mathit{Y}=\mathbf{\text{ABC}}+\mathit{E}$$(2.4) , the following model is set up: (2.5) $$\mathit{Y}=\mathbf{\text{ABC}}+\mathbf{\text{UZ}}+\mathit{E}$$(2.5) where $\mathit{U}\sim N_{p,l}(\mathbf{0},{\mathbf{\Sigma }}_u,{\mathit{I}}_l)$ is independent of E and Z is related to C because, in the model, usually via C, the observations are linked to the subareas and Z helps to describe the random variation (errors) in the subareas. Therefore, formally, $\mathcal{C}(\mathit{Z}\prime )\subseteq \mathcal{C}(\mathit{C}\prime ),$ where $\mathcal{C}(•)$ denotes the column vector space. Additionally we will standardize these random effects by assuming $\mathit{Z}{\mathit{V}}^{-1}\mathit{Z}\prime ={\mathit{I}}_l,$ meaning that it has been adjusted for a different number of subareas within an area, taking into an account the uncertainty in the survey estimates. This has two important implications. One is that it simplifies the mathematical treatment of the model and the second is that it becomes easier to compare predicted random effects.

Now observed covariables will be introduced in the model. There exist different types of covariables but we are not going to distinguish between different types. For example one type of covariables are those which are accounted for in the survey design (e.g. variables defining different strata) and another type can be covariables observed in registries and a third type can be baseline data which were available when a survey study with repeated measurements on units was initiated. The effect of the covariate will be modeled by a term ${\mathit{B}}_1{\mathit{C}}_1$ where ${\mathit{B}}_1$ stands for the unknown effect and ${\mathit{C}}_1$ collects the observed covariables leading to the model (2.6) $$\mathit{Y}=\mathbf{\text{ABC}}+{\mathit{B}}_1{\mathit{C}}_1+\mathbf{\text{UZ}}+\mathit{E}$$(2.6)

Moreover, suppose that there are a number of latent variables. It can be large clusters of variables, e.g. weather variables, socio-econometric variables, chemical characteristics of soil and water, etc. A typical phenomena of these clusters are that one knows that there should be an effect on the survey estimates but it is difficult to express this via some functional relationship. Instead we will model the relationship via rank restrictions on parameter matrices and will consider the following model: (2.7) $$\mathit{Y}=\mathbf{\text{ABC}}+{\mathit{B}}_1{\mathit{C}}_1+\mathbf{\Psi }\mathit{F}+\mathbf{\text{UZ}}+\mathit{E}$$(2.7) where F: t × n is known and $\mathbf{\Psi }$ is unknown of rank $r(\mathbf{\Psi })=r\le \min (p,t).$ The idea is that all clustered variables appear in F and these variables are governed by r unobserved latent variables (processes) which are incorporated in the model via the rank restrictions on $\mathbf{\Psi }.$ It is important to note that p has to be larger than r meaning that modeling the number of underlying latent variables depends on how many times the survey has been repeated. Thus, we really see an advantage of performing repeated surveys when latent variables are supposed to exist, which indeed we believe is a common phenomena.

3. Detailed model specification

Let Y consist of the direct estimates of the subareas where each column corresponds to a subarea including p estimators from the p repeated measurements. When performing likelihood inference the direct estimates are used but in the presentation we will not distinguish between estimators and estimates.

Definition 3.1.

Let $$\mathit{Y}=\mathbf{\text{ABC}}+{\mathit{B}}_1{\mathit{C}}_1+\mathbf{\Psi }\mathit{F}+\mathbf{\text{UZ}}+\mathit{E}$$ where A: p × q, C: k × n, ${\mathit{C}}_1$: $k_1\times n,$ F: t × n, are all known matrices, $\mathcal{C}(\mathit{F}\prime )\subseteq \mathcal{C}(\mathit{C}\prime ),r(\mathbf{\Psi })=r<\min (p,t),$ Z: l × n, $p\le l<n,\mathit{Z}{\mathit{V}}^{-1}\mathit{Z}\prime ={\mathit{I}}_l,\mathit{U}\sim N_{p,l}(\mathbf{0},{\mathbf{\Sigma }}_u,{\mathit{I}}_l),\mathit{E}\sim N_{p,n}(\mathbf{0},{\mathbf{\Sigma }}_e,\mathit{V})$ where U and E are independently distributed. The parameters ${\mathbf{\Sigma }}_u$ and ${\mathbf{\Sigma }}_e$ are supposed to be unknown and positive definite whereas V is known and determined from the survey.

When in Definition 3.1 ${\mathit{B}}_1{\mathit{C}}_1=\mathbf{0},$ $\mathbf{\Psi }\mathit{F}=\mathbf{0}$ and $\mathbf{\text{UZ}}=\mathbf{0}$ then we have the classical growth curve model (GMANOVA) (see Potthoff and Roy 1964; von Rosen 2018). When $\mathbf{\Psi }\mathit{F}=\mathbf{0}$ and $\mathbf{\text{UZ}}=\mathbf{0}$ then the growth curve model with covariate information appears which sometimes is identified as a mixture of the GMANOVA and MANOVA models. Some references to this model, as well as more general models (see Chinchilli and Elswick 1985; Verbyla and Venables 1988; von Rosen 1989; and Bai and Shi 2007). If ${\mathit{B}}_1{\mathit{C}}_1=\mathbf{0}$ and $\mathbf{\Psi }\mathit{F}=\mathbf{0}$ we say that we have the growth curve model with random effects (see Ip et al. 2007) whereas if only $\mathbf{\Psi }\mathit{F}=\mathbf{0}$ holds (see Yokoyama and Fujikoshi 1992; Yokoyama 1995) where similar models are considered and where references to earlier works can be found. However, usually in these works one puts structures on the covariance matrices which will not be discussed in this article.

4. Estimation of parameters

Let Q be a matrix of basis vectors such that $\mathcal{C}(\mathit{Q})=\mathcal{C}({\mathit{V}}^{-1/2}\mathit{C}\prime :{\mathit{V}}^{-1/2}{\mathit{C}}_1^{\prime }:{\mathit{V}}^{-1/2}\mathit{Z}\prime ),$ where ${\mathit{V}}^{1/2}$ is a symmetric square root of V in $\mathit{E}\sim N_{p,n}(\mathbf{0},\mathbf{\Sigma },\mathit{V}).$ Further, let ${\mathit{Q}}^o$ be any matrix of full rank satisfying $\mathcal{C}({\mathit{Q}}^o)=\mathcal{C}{({\mathit{C}}^{\prime }:{\mathit{C}}_1^{\prime }:{\mathit{Z}}^{\prime })}^{\perp }.$ We assume $\mathit{Q}\prime \mathit{Q}={\mathit{I}}_v,$ where $v=r(\mathit{C}\prime :{\mathit{C}}_1^{\prime }:\mathit{Z}\prime ),$ v > l.

A one-one transformation of the model in Definition 3.1, using Q, yields (4.8) $$\begin{array}{c}\mathit{Y}{\mathit{V}}^{-1/2}\mathit{Q}=\mathbf{\text{ABC}}{\mathit{V}}^{-1/2}\mathit{Q}+{\mathit{B}}_1{\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}+\mathbf{\Psi }\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q}\\ +\mathbf{\text{UZ}}{\mathit{V}}^{-1/2}\mathit{Q}+\mathit{E}{\mathit{V}}^{-1/2}\mathit{Q},\end{array}$$(4.8) (4.9) $$\mathit{Y}{\mathit{V}}^{-1/2}{\mathit{Q}}^o=\mathit{E}{\mathit{V}}^{-1/2}{\mathit{Q}}^o$$(4.9)

The matrix $\mathit{G}=\mathit{Q}\prime {\mathit{V}}^{-1/2}\mathit{Z}\prime \mathit{Z}{\mathit{V}}^{-1/2}\mathit{Q}$ of size v × v is idempotent since $\mathcal{C}({\mathit{V}}^{-1/2}\mathit{Z}\prime )\subseteq \mathcal{C}(\mathit{Q})$ and $\mathit{Z}{\mathit{V}}^{-1}\mathit{Z}\prime =\mathit{I},$ and therefore $$\mathit{G}=\mathbf{\Gamma }\left(\begin{array}{cc}{\mathit{I}}_l& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right)\mathbf{\Gamma }\prime$$ where $\mathbf{\Gamma }:$ v × v is a known orthogonal matrix. The identity in Equation (4.8)(4.8) $$\begin{array}{c}\mathit{Y}{\mathit{V}}^{-1/2}\mathit{Q}=\mathbf{\text{ABC}}{\mathit{V}}^{-1/2}\mathit{Q}+{\mathit{B}}_1{\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}+\mathbf{\Psi }\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q}\\ +\mathbf{\text{UZ}}{\mathit{V}}^{-1/2}\mathit{Q}+\mathit{E}{\mathit{V}}^{-1/2}\mathit{Q},\end{array}$$(4.8) is post-multiplied by $\mathbf{\Gamma }$ leading to the model (4.10) $$\begin{array}{c}\mathit{Y}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }=\mathbf{\text{ABC}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+{\mathit{B}}_1{\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+\mathbf{\Psi }\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }\\ +\mathbf{\text{UZ}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+\mathit{E}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }.\end{array}$$(4.10)

However, since $\mathbf{\Gamma }\prime \mathit{Q}\prime {\mathit{V}}^{-1/2}\mathit{Z}\prime \mathit{Z}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }=\mathit{I},$ $$D\left[\mathbf{\text{UZ}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }\right]=\left(\begin{array}{cc}{\mathit{I}}_l& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right)\otimes {\mathbf{\Sigma }}_u$$ and this leads to that the model in Equation (4.10)(4.10) $$\begin{array}{c}\mathit{Y}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }=\mathbf{\text{ABC}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+{\mathit{B}}_1{\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+\mathbf{\Psi }\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }\\ +\mathbf{\text{UZ}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+\mathit{E}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }.\end{array}$$(4.10) will be split into two models. Let $\mathbf{\Gamma }=({\mathbf{\Gamma }}_1:{\mathbf{\Gamma }}_2):$ $v\times l,v\times (v-l).$ Then we have three models which will be used when finding estimators:

1. $\mathit{Y}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_1=\mathbf{\text{ABC}}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_1+{\mathit{B}}_1{\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_1+\mathbf{\Psi }\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_1\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\mathbf{\text{UZ}}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_1+\mathit{E}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_1,\hspace{1em}\hspace{1em}\mathbf{\text{UZ}}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_1\sim N_{p,l}(\mathbf{0},{\mathbf{\Sigma }}_u,{\mathit{I}}_l),\hspace{1em}\mathit{E}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_1\sim N_{p,l}(\mathbf{0},{\mathbf{\Sigma }}_e,{\mathit{I}}_l),\hspace{1em}\hspace{1em}\hspace{1em}$(4.11)

2. $\mathit{Y}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_2=\mathbf{\text{ABC}}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_2+{\mathit{B}}_1{\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_2+\mathbf{\Psi }\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_2+ \mathit{E}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_2,\hspace{1em}\mathit{E}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_2\sim N_{p,v-l}(\mathbf{0},{\mathbf{\Sigma }}_e,{\mathit{I}}_{v-l}),$(4.12)

3. $\mathit{Y}{\mathit{V}}^{-1/2}{\mathit{Q}}^o=\mathit{E}{\mathit{V}}^{-1/2}{\mathit{Q}}^o,\hspace{1em}\mathit{E}{\mathit{V}}^{-1/2}{\mathit{Q}}^o\sim N_{p,n-v}(\mathbf{0},{\mathbf{\Sigma }}_e,{\mathit{I}}_{n-v}).$(4.13)

The idea is to utilize Equations (4.12)(4.10) $$\begin{array}{c}\mathit{Y}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }=\mathbf{\text{ABC}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+{\mathit{B}}_1{\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+\mathbf{\Psi }\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }\\ +\mathbf{\text{UZ}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+\mathit{E}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }.\end{array}$$(4.10) and (4.13)(4.10) $$\begin{array}{c}\mathit{Y}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }=\mathbf{\text{ABC}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+{\mathit{B}}_1{\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+\mathbf{\Psi }\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }\\ +\mathbf{\text{UZ}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+\mathit{E}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }.\end{array}$$(4.10) when estimating B, ${\mathit{B}}_1,\mathbf{\Psi }$ and ${\mathbf{\Sigma }}_e.$ To estimate ${\mathbf{\Sigma }}_u$ these estimators are inserted in Equation (4.11)(4.10) $$\begin{array}{c}\mathit{Y}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }=\mathbf{\text{ABC}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+{\mathit{B}}_1{\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+\mathbf{\Psi }\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }\\ +\mathbf{\text{UZ}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+\mathit{E}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }.\end{array}$$(4.10) . Suppose that the estimators $\widehat{\mathit{B}},{\widehat{\mathit{B}}}_1,$ $\widehat{\mathbf{\Psi }}$ and ${\widehat{\mathbf{\Sigma }}}_e$ have been obtained, and let $$\begin{array}{c}{\mathit{Y}}_0=\mathit{X}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_1-\mathit{A}\widehat{\mathit{B}}\mathit{C}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_1-{\widehat{\mathit{B}}}_1{\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_1-\widehat{\mathbf{\Psi }}\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q},\\ \mathbf{\Omega }={\mathbf{\Sigma }}_u+{\mathbf{\Sigma }}_e.\end{array}$$

Under the assumption of no randomness in $\widehat{\mathit{B}},$ ${\widehat{\mathit{B}}}_1$ and $\widehat{\mathbf{\Psi }},$ $${\mathit{Y}}_0=\tilde{\mathit{E}},\hspace{1em}\tilde{\mathit{E}}\sim N_{p,l}(0,\mathbf{\Omega },{\mathit{I}}_l)$$ Thus, if $p\le l,\widehat{\mathbf{\Omega }}=\frac{1}{l}{\mathit{Y}}_0{\mathit{Y}}_0^{\prime }$ , implying $${\widehat{\mathbf{\Sigma }}}_u=\frac{1}{l}{\mathit{Y}}_0{\mathit{Y}}_0^{\prime }-{\widehat{\mathbf{\Sigma }}}_e.$$

Here it is assumed that the difference is positive definite. If ${\widehat{\mathbf{\Sigma }}}_u$ is not positive definite, which should occur with some positive probability, this indicates that data does not support the model presented in Definition 3.1 or the estimation procedure has some deficiencies. It can be noted that the problem with negative variance components in mixed linear models has a long history and an early reference is Nelder (1954). A more recent contribution to the discussion of the problem is provided by Molenberghs and Verbeke (2011). If we would start with the model in Equation (2.5)(2.5) $$\mathit{Y}=\mathbf{\text{ABC}}+\mathbf{\text{UZ}}+\mathit{E}$$(2.5) and only consider $D[\mathit{Y}]=\mathbf{\text{ZZ}}\prime \otimes {\mathbf{\Sigma }}_u+\mathit{V}\otimes {\mathbf{\Sigma }}_e$ we could avoid to assume ${\mathbf{\Sigma }}_u$ to be positive definite. However, in our model derivation, via $\mathit{U}\sim N_{p,l}(\mathbf{0},{\mathbf{\Sigma }}_u,{\mathit{I}}_l),$ ${\mathbf{\Sigma }}_u$ has to be positive definite. Thus, if obtaining a non-positive definite ${\widehat{\mathbf{\Sigma }}}_u$ we have to either reformulate the model or change the estimator. One way to modify ${\widehat{\mathbf{\Sigma }}}_u$ is to work with the eigenvalues of ${\widehat{\mathbf{\Sigma }}}_u.$ If there are a few “small” (by absolute value) non-positive eigenvalues these eigenvalues can be replaced by small positive quantities. This should of course take place with caution and common sense but as far as we know there is no general recipe of how to handle a non-positive definite ${\widehat{\mathbf{\Sigma }}}_u.$

To estimate $\widehat{\mathit{B}},{\widehat{\mathit{B}}}_1,$ $\widehat{\mathbf{\Psi }}$ and ${\widehat{\mathbf{\Sigma }}}_e$ Equations (4.12)(4.10) $$\begin{array}{c}\mathit{Y}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }=\mathbf{\text{ABC}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+{\mathit{B}}_1{\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+\mathbf{\Psi }\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }\\ +\mathbf{\text{UZ}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+\mathit{E}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }.\end{array}$$(4.10) and (4.13)(4.10) $$\begin{array}{c}\mathit{Y}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }=\mathbf{\text{ABC}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+{\mathit{B}}_1{\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+\mathbf{\Psi }\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }\\ +\mathbf{\text{UZ}}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }+\mathit{E}{\mathit{V}}^{-1/2}\mathit{Q}\mathbf{\Gamma }.\end{array}$$(4.10) are merged: (4.14) $$\begin{array}{c}\mathit{Y}{\mathit{V}}^{-1/2}(\mathit{Q}{\mathbf{\Gamma }}_2:{\mathit{Q}}^o)=\mathbf{\text{AB}}(\mathit{C}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_2:\mathbf{0})+{\mathit{B}}_1({\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_2:\mathbf{0})\\ +\mathbf{\Psi }(\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_2:\mathbf{0})+\tilde{\mathit{E}},\end{array}$$(4.14) where $\tilde{\mathit{E}}\sim N_{p,n-l}(\mathbf{0},{\mathbf{\Sigma }}_e,{\mathit{I}}_{n-l}).$ Put $$\begin{array}{c}\mathit{X}=\mathit{Y}{\mathit{V}}^{-1/2}(\mathit{Q}{\mathbf{\Gamma }}_2:{\mathit{Q}}^o):p\times (n-l),\\ {\mathit{D}}_1=(\mathit{C}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_2:\mathbf{0}):k\times (n-l),\\ {\mathit{D}}_2=({\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_2:\mathbf{0}):k_1\times (n-l),\\ {\mathit{D}}_3=(\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q}{\mathbf{\Gamma }}_2:\mathbf{0}):t\times (n-l).\end{array}$$

Equation (4.14)(4.8) $$\begin{array}{c}\mathit{Y}{\mathit{V}}^{-1/2}\mathit{Q}=\mathbf{\text{ABC}}{\mathit{V}}^{-1/2}\mathit{Q}+{\mathit{B}}_1{\mathit{C}}_1{\mathit{V}}^{-1/2}\mathit{Q}+\mathbf{\Psi }\mathit{F}{\mathit{V}}^{-1/2}\mathit{Q}\\ +\mathbf{\text{UZ}}{\mathit{V}}^{-1/2}\mathit{Q}+\mathit{E}{\mathit{V}}^{-1/2}\mathit{Q},\end{array}$$(4.8) is identical to (4.15) $$\mathit{X}=\mathbf{\text{AB}}{\mathit{D}}_1+{\mathit{B}}_1{\mathit{D}}_2+\mathbf{\Psi }{\mathit{D}}_3+\tilde{\mathit{E}}$$(4.15) which is an extended Growth Curve model (GMANOVA + MANOVA) with a reduced rank regression component. Furthermore, the likelihood function which corresponds to the model in Equation (4.15)(4.15) $$\mathit{X}=\mathbf{\text{AB}}{\mathit{D}}_1+{\mathit{B}}_1{\mathit{D}}_2+\mathbf{\Psi }{\mathit{D}}_3+\tilde{\mathit{E}}$$(4.15) equals (4.16) $$\begin{array}{c}L(\mathit{B},{\mathit{B}}_1,\mathbf{\Psi },{\mathbf{\Sigma }}_e)={(2\pi )}^{-\frac{1}{2}p(n-l)}|{\mathbf{\Sigma }}_e{|}^{-\frac{1}{2}(n-l)}\\ \hspace{1em}\hspace{1em}\times \hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{1}{2}\text{tr}\left\{{\mathbf{\Sigma }}_e^{-1}\right(\mathit{X}-\text{AB}{\mathit{D}}_1-{\mathit{B}}_1{\mathit{D}}_2-\mathbf{\Psi }{\mathit{D}}_3)()\prime \}\right),\end{array}$$(4.16) where we have used the convention that instead of $(\mathit{H})(\mathit{H})\prime$ it is written $(\mathit{H})()\prime$ for any arbitrary matrix expression H. Our first observation is that the likelihood in Equation (4.16)(4.16) $$\begin{array}{c}L(\mathit{B},{\mathit{B}}_1,\mathbf{\Psi },{\mathbf{\Sigma }}_e)={(2\pi )}^{-\frac{1}{2}p(n-l)}|{\mathbf{\Sigma }}_e{|}^{-\frac{1}{2}(n-l)}\\ \hspace{1em}\hspace{1em}\times \hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{1}{2}\text{tr}\left\{{\mathbf{\Sigma }}_e^{-1}\right(\mathit{X}-\text{AB}{\mathit{D}}_1-{\mathit{B}}_1{\mathit{D}}_2-\mathbf{\Psi }{\mathit{D}}_3)()\prime \}\right),\end{array}$$(4.16) is smaller or equal to (4.17) $$\begin{array}{c}{(2\pi )}^{-\frac{1}{2}p(n-l)}|{\mathbf{\Sigma }}_e{|}^{-\frac{1}{2}(n-l)}\\ \times \hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{1}{2}\text{tr}\left\{{\mathbf{\Sigma }}_e^{-1}\right((\mathit{X}-\mathbf{\text{AB}}{\mathit{D}}_1)(\mathit{I}-{\mathit{P}}_{D_2^{\prime }})-\mathbf{\Psi }{\mathit{D}}_3(\mathit{I}-{\mathit{P}}_{D_2^{\prime }})()\prime \left)\right\}\right)\hspace{1em}\hspace{1em}\end{array}$$(4.17) with equality if and only if $${\mathit{B}}_1{\mathit{D}}_2=\mathit{X}{\mathit{P}}_{D_2^{\prime }}-\mathbf{\text{AB}}{\mathit{D}}_1{\mathit{P}}_{D_2^{\prime }}-\mathbf{\Psi }{\mathit{D}}_3{\mathit{P}}_{D_2^{\prime }}$$ which, under some full rank conditions on ${\mathit{D}}_2,$ determines ${\mathit{B}}_1$ as a function of B and $\mathbf{\Psi }.$ Thus, if B and $\mathbf{\Psi }$ can be estimated the covariate effect described by ${\mathit{B}}_1{\mathit{D}}_2$ can be estimated. The density in Equation (4.17)(4.17) $$\begin{array}{c}{(2\pi )}^{-\frac{1}{2}p(n-l)}|{\mathbf{\Sigma }}_e{|}^{-\frac{1}{2}(n-l)}\\ \times \hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{1}{2}\text{tr}\left\{{\mathbf{\Sigma }}_e^{-1}\right((\mathit{X}-\mathbf{\text{AB}}{\mathit{D}}_1)(\mathit{I}-{\mathit{P}}_{D_2^{\prime }})-\mathbf{\Psi }{\mathit{D}}_3(\mathit{I}-{\mathit{P}}_{D_2^{\prime }})()\prime \left)\right\}\right)\hspace{1em}\hspace{1em}\end{array}$$(4.17) corresponds to the model (4.18) $$\begin{array}{c}\mathit{X}(\mathit{I}-{\mathit{P}}_{D_2^{\prime }})=\mathbf{\text{AB}}{\mathit{D}}_1(\mathit{I}-{\mathit{P}}_{D_2^{\prime }})+\mathbf{\Psi }{\mathit{D}}_3(\mathit{I}-{\mathit{P}}_{D_2^{\prime }})+\tilde{\mathit{E}},\\ \tilde{\mathit{E}}\sim N_{p,n-l}(\mathbf{0},{\mathbf{\Sigma }}_e,{\mathit{I}}_{n-l})\end{array}$$(4.18) which is a Growth Curve model with latent variables (reduced rank effect). The model has been considered (Reinsel and Velu 1998; von Rosen and von Rosen 2017), where also other references can be found. Now it is shortly described how B and $\mathbf{\Psi }$ can be estimated. For notational conveniences the model in Equation (4.18)(4.18) $$\begin{array}{c}\mathit{X}(\mathit{I}-{\mathit{P}}_{D_2^{\prime }})=\mathbf{\text{AB}}{\mathit{D}}_1(\mathit{I}-{\mathit{P}}_{D_2^{\prime }})+\mathbf{\Psi }{\mathit{D}}_3(\mathit{I}-{\mathit{P}}_{D_2^{\prime }})+\tilde{\mathit{E}},\\ \tilde{\mathit{E}}\sim N_{p,n-l}(\mathbf{0},{\mathbf{\Sigma }}_e,{\mathit{I}}_{n-l})\end{array}$$(4.18) is written (4.19) $${\mathit{X}}_1=\mathbf{\text{AB}}{\tilde{\mathit{D}}}_1+\mathbf{\Psi }{\tilde{\mathit{D}}}_3+\tilde{\mathit{E}},\hspace{1em}\tilde{\mathit{E}}\sim N_{p,n-l}(\mathbf{0},{\mathbf{\Sigma }}_e,{\mathit{I}}_{n-l})$$(4.19)

Note that $\mathcal{C}({\tilde{\mathit{D}}}_3^{\prime })\subseteq \mathcal{C}({\tilde{\mathit{D}}}_1^{\prime })$ because $\mathcal{C}(\mathit{F}\prime )\subseteq \mathcal{C}(\mathit{C}\prime )$ which is essential for being able to obtain explicit estimators.

Let $$\mathit{W}=({\mathit{X}}_1-\mathbf{\text{AB}}{\tilde{\mathit{D}}}_1-\mathbf{\Psi }{\tilde{\mathit{D}}}_3)()\prime$$

and then the likelihood function for model Equation (4.19)(4.19) $${\mathit{X}}_1=\mathbf{\text{AB}}{\tilde{\mathit{D}}}_1+\mathbf{\Psi }{\tilde{\mathit{D}}}_3+\tilde{\mathit{E}},\hspace{1em}\tilde{\mathit{E}}\sim N_{p,n-l}(\mathbf{0},{\mathbf{\Sigma }}_e,{\mathit{I}}_{n-l})$$(4.19) equals (4.20) $$L(\mathit{B},\mathbf{\Psi },{\mathbf{\Sigma }}_e)={(2\pi )}^{-(n-l)p/2}|{\mathbf{\Sigma }}_e{|}^{-(n-l)/2}\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{1}{2}\text{tr}\left\{{\mathbf{\Sigma }}_e^{-1}\mathit{W}\right\}\right)$$(4.20)

The upper bound of the likelihood function in Equation (4.20)$$\mathit{W}=({\mathit{X}}_1-\mathbf{\text{AB}}{\tilde{\mathit{D}}}_1-\mathbf{\Psi }{\tilde{\mathit{D}}}_3)()\prime$$ is well known (see Srivastava and Khatri 1979, Theorem 1.10.4): (4.21) $$L(\mathit{B},\mathbf{\Psi },{\mathbf{\Sigma }}_e)\le {(2\pi )}^{-(n-l)p/2}{e}^{\{-\frac{1}{2}(n-l\left)p\right\}}|\frac{1}{n-l}\mathit{W}{|}^{-(n-l)/2}$$(4.21) and in Equation (4.21)(4.21) $$L(\mathit{B},\mathbf{\Psi },{\mathbf{\Sigma }}_e)\le {(2\pi )}^{-(n-l)p/2}{e}^{\{-\frac{1}{2}(n-l\left)p\right\}}|\frac{1}{n-l}\mathit{W}{|}^{-(n-l)/2}$$(4.21) equality holds if and only if $(n-l){\mathbf{\Sigma }}_e=\mathit{W}.$ Hence, if B and $\mathbf{\Psi }$ are estimated then ${\mathbf{\Sigma }}_e$ is also estimated. Maximizing the log-likelihood function is equivalent to minimizing the determinant $\left|\mathit{W}\right|$ which now will take place.

It is worth noting that it is used, since $r(\mathbf{\Psi })=r,$ that the matrix $\mathbf{\Psi }$ can be factored into two terms, i.e. $\mathbf{\Psi }={\mathbf{\Psi }}_1{\mathbf{\Psi }}_2,$ where ${\mathbf{\Psi }}_1$ and ${\mathbf{\Psi }}_2$ are of size p × r and size r × t, respectively. Depending on the knowledge about the model ${\mathbf{\Psi }}_1$ and ${\mathbf{\Psi }}_2$ can be interpreted but there are also cases where the factorization has no clear interpretation.

In the subsequent, let ${\mathit{P}}_Q=\mathit{Q}{({\mathit{Q}}^{\prime }\mathit{Q})}^{-}\mathit{Q}\prime$ and ${\mathit{P}}_{Q,V}=\mathit{Q}{({\mathit{Q}}^{\prime }{\mathit{V}}^{-1}\mathit{Q})}^{-}\mathit{Q}\prime {\mathit{V}}^{-1}$ with ${(•)}^{-}$ denoting an arbitrary generalized inverse.

Following the approach presented in (Kollo and von Rosen 2005, Chapter 4) we have (4.22) $$\left|\mathit{W}\right|\ge |{\mathit{S}}_1+{\mathit{T}}_1(\mathit{X}{\mathit{P}}_{{\tilde{D}}_1^{\prime }}-{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3\left)\right(\mathit{X}{\mathit{P}}_{{\tilde{D}}_1^{\prime }}-{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3)\prime \mathit{T}{\prime }_1|$$(4.22) where (4.23) $${\mathit{S}}_1={\mathit{X}}_1({\mathit{I}}_{n-l}-{\mathit{P}}_{{\tilde{D}}_1^{\prime }}){\mathit{X}}_1^{\prime },$$(4.23) (4.24) $${\mathit{T}}_1={\mathit{I}}_p-{\mathit{P}}_{A,S_1}$$(4.24) and equality in Equation (4.22)(4.22) $$\left|\mathit{W}\right|\ge |{\mathit{S}}_1+{\mathit{T}}_1(\mathit{X}{\mathit{P}}_{{\tilde{D}}_1^{\prime }}-{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3\left)\right(\mathit{X}{\mathit{P}}_{{\tilde{D}}_1^{\prime }}-{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3)\prime \mathit{T}{\prime }_1|$$(4.22) holds if and only if $\text{AB}{\tilde{\mathit{D}}}_1=\mathit{P}{}_{A,S_1}(\mathit{X}{\mathit{P}}_{{\tilde{D}}_1^{\prime }}-{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3).$ Throughout the article it is supposed that ${\mathit{S}}_1$ is positive definite which holds with probability 1. Thus, B as a function of $\mathbf{\Psi }={\mathbf{\Psi }}_1{\mathbf{\Psi }}_2$ can be obtained. Moreover, (4.25) $$\begin{array}{c}|{\mathit{S}}_1+{\mathit{T}}_1(\mathbf{\text{XP}}{}_{{\tilde{D}}_1^{\prime }}-{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3\left)\right(\mathbf{\text{XP}}{}_{{\tilde{D}}_1^{\prime }}-{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3)\prime \mathit{T}{\prime }_1|\\ \ge |{\mathit{S}}_2+{\mathit{T}}_2{\mathit{T}}_1\mathbf{\text{XP}}{}_{{\tilde{D}}_3^{\prime }}\mathit{X}\prime \mathit{T}{\prime }_1\mathit{T}{\prime }_2|,\end{array}$$(4.25) where (4.26) $$\begin{array}{c}{\mathit{S}}_2={\mathit{S}}_1+{\mathit{T}}_1\mathit{X}(\mathit{P}{}_{{\tilde{D}}_1^{\prime }}-\mathit{P}{}_{{\tilde{D}}_3^{\prime }})\mathit{X}\prime \mathit{T}{\prime }_1,\\ {\mathit{T}}_2={\mathit{I}}_p-\mathit{P}{}_{T_1{\Psi }_1,S_2},\end{array}$$(4.26)

and the equality in Equation (4.25)(4.25) $$\begin{array}{c}|{\mathit{S}}_1+{\mathit{T}}_1(\mathbf{\text{XP}}{}_{{\tilde{D}}_1^{\prime }}-{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3\left)\right(\mathbf{\text{XP}}{}_{{\tilde{D}}_1^{\prime }}-{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3)\prime \mathit{T}{\prime }_1|\\ \ge |{\mathit{S}}_2+{\mathit{T}}_2{\mathit{T}}_1\mathbf{\text{XP}}{}_{{\tilde{D}}_3^{\prime }}\mathit{X}\prime \mathit{T}{\prime }_1\mathit{T}{\prime }_2|,\end{array}$$(4.25) holds if and only if (4.27) $${\mathit{T}}_1{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3=\mathit{P}{}_{T_1{\Psi }_1,S_2}\text{XP}{}_{{\tilde{D}}_3^{\prime }}$$(4.27)

Given ${\mathbf{\Psi }}_1,$ the linear system of equations in Equation (4.27)(4.27) $${\mathit{T}}_1{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3=\mathit{P}{}_{T_1{\Psi }_1,S_2}\text{XP}{}_{{\tilde{D}}_3^{\prime }}$$(4.27) is consistent and hence ${\mathbf{\Psi }}_2$ can always be estimated as a function of ${\mathbf{\Psi }}_1:$ (4.28) $${\mathbf{\Psi }}_2={\left({\mathbf{\Psi }}_1^{\prime }{\mathit{T}}_1^{\prime }{\mathit{S}}_2^{-1}{\mathit{T}}_1{\mathbf{\Psi }}_1\right)}^{-1}\mathbf{\Psi }{\prime }_1\mathit{T}{\prime }_1{\mathit{S}}_2^{-1}\mathit{X}{\tilde{\mathit{D}}}_3^{\prime }{\left({\tilde{\mathit{D}}}_3{\tilde{\mathit{D}}}_3^{\prime }\right)}^{-1},$$(4.28) under the assumption that $r({\tilde{\mathit{D}}}_3)=k_1$ and $\mathcal{C}(\mathit{A})\cap \mathcal{C}({\mathbf{\Psi }}_1)=\left\{0\right\}.$ The only parameter in the right-hand side of Equation (4.25)(4.25) $$\begin{array}{c}|{\mathit{S}}_1+{\mathit{T}}_1(\mathbf{\text{XP}}{}_{{\tilde{D}}_1^{\prime }}-{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3\left)\right(\mathbf{\text{XP}}{}_{{\tilde{D}}_1^{\prime }}-{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3)\prime \mathit{T}{\prime }_1|\\ \ge |{\mathit{S}}_2+{\mathit{T}}_2{\mathit{T}}_1\mathbf{\text{XP}}{}_{{\tilde{D}}_3^{\prime }}\mathit{X}\prime \mathit{T}{\prime }_1\mathit{T}{\prime }_2|,\end{array}$$(4.25) which is left to estimate is ${\mathbf{\Psi }}_1.$ The right-hand side of Equation (4.25)(4.25) $$\begin{array}{c}|{\mathit{S}}_1+{\mathit{T}}_1(\mathbf{\text{XP}}{}_{{\tilde{D}}_1^{\prime }}-{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3\left)\right(\mathbf{\text{XP}}{}_{{\tilde{D}}_1^{\prime }}-{\mathbf{\Psi }}_1{\mathbf{\Psi }}_2{\tilde{\mathit{D}}}_3)\prime \mathit{T}{\prime }_1|\\ \ge |{\mathit{S}}_2+{\mathit{T}}_2{\mathit{T}}_1\mathbf{\text{XP}}{}_{{\tilde{D}}_3^{\prime }}\mathit{X}\prime \mathit{T}{\prime }_1\mathit{T}{\prime }_2|,\end{array}$$(4.25) can be factored as $$|{\mathit{S}}_2+{\mathit{T}}_2{\mathit{T}}_1\mathbf{\text{XP}}{}_{{\tilde{D}}_3^{\prime }}\mathit{X}\prime \mathit{T}{\prime }_1\mathit{T}{\prime }_2|=\left|{\mathit{S}}_2\right|\left|\mathit{R}\right||\mathit{M}\prime \mathbf{\text{UM}}|,$$ where $$\begin{array}{c}\mathit{R}={\mathit{I}}_n+\mathit{P}{}_{{\tilde{\mathit{D}}}_3^{\prime }}\mathit{X}\prime \mathit{T}{\prime }_1{\mathit{S}}_2^{-1}{\mathit{T}}_1\mathbf{\text{XP}}{}_{{\tilde{D}}_3^{\prime }},\\ \mathit{M}=\mathit{H}\prime {\mathbf{\Psi }}_1{\left({\mathbf{\Psi }}_1^{\prime }{\mathbf{\text{HH}}}^{\prime }{\mathbf{\Psi }}_1\right)}^{-1/2},\\ \mathit{U}={\mathit{I}}_{p-r\left(\mathit{A}\right)}-\mathit{H}\prime \mathbf{\text{XP}}{}_{{\tilde{D}}_3^{\prime }}{\mathit{R}}^{-1}\mathit{P}{}_{{\tilde{D}}_3^{\prime }}\mathit{X}\prime \mathit{H},\end{array}$$ with $\mathit{H}:p\times r(\mathit{A})$ is such that $\mathit{T}{\prime }_1{\mathit{S}}_2^{-1}{\mathit{T}}_1=\mathbf{\text{HH}}\prime ,$ $\mathit{M}:(p-r(\mathit{A}))\times r$ is semi-orthogonal of rank r and U is a positive definite matrix. The square root in M is supposed to be symmetric. Due to the Poincaré separation theorem (see Rao 1979) $$\left|{\mathit{S}}_2\right|\left|\mathit{R}\right||\mathit{M}\prime \mathbf{\text{UM}}|\ge \left|{\mathit{S}}_2\right|\left|\mathit{R}\right|\prod _{i=1}^r{\lambda }_{p-r\left(\mathit{A}\right)-r+i},$$ where ${\lambda }_1\ge \cdots \ge {\lambda }_{p-r(\mathit{A})}$ are the eigenvalues of U which do not depend on any unknown parameter. The lower bound is then attained if M consists of the corresponding eigenvectors ${\mathit{v}}_1,\dots ,{\mathit{v}}_r.$ Let $\tilde{\mathit{M}}=({\mathit{v}}_1,\dots ,{\mathit{v}}_r).$ Therefore, a maximum likelihood estimator of ${\mathbf{\Psi }}_1$ is found if we can find ${\mathbf{\Psi }}_1$ satisfying the following equality: (4.29) $$\mathit{H}\prime {\mathbf{\Psi }}_1{\left({\mathbf{\Psi }}_1^{\prime }{\mathbf{\text{HH}}}^{\prime }{\mathbf{\Psi }}_1\right)}^{-1/2}=\tilde{\mathit{M}}.$$(4.29)

Since $\tilde{\mathit{M}}\prime \tilde{\mathit{M}}={\mathit{I}}_r,$ the following choice of ${\mathbf{\Psi }}_1$ is appropriate: (4.30) $${\widehat{\mathbf{\Psi }}}_1=\mathit{H}{\left({\mathit{H}}^{\prime }\mathit{H}\right)}^{-1}\tilde{\mathit{M}}.$$(4.30)

There remains an unresolved problem of presenting all solutions to Equation (4.29)(4.29) $$\mathit{H}\prime {\mathbf{\Psi }}_1{\left({\mathbf{\Psi }}_1^{\prime }{\mathbf{\text{HH}}}^{\prime }{\mathbf{\Psi }}_1\right)}^{-1/2}=\tilde{\mathit{M}}.$$(4.29) , but this is not considered here.

Together with Equation (4.28)$$\left|{\mathit{S}}_2\right|\left|\mathit{R}\right||\mathit{M}\prime \mathbf{\text{UM}}|\ge \left|{\mathit{S}}_2\right|\left|\mathit{R}\right|\prod _{i=1}^r{\lambda }_{p-r\left(\mathit{A}\right)-r+i},$$ this yields that $\widehat{\Psi }{\tilde{\mathit{D}}}_3={\widehat{\mathbf{\Psi }}}_1\widehat{\mathbf{\Psi }}{\prime }_1\mathit{T}{\prime }_1{\mathit{S}}_2^{-1}\mathit{X}{\mathit{P}}_{\tilde{D}{\prime }_3}.$ Now, the obtained results will be stated in the following theorem.

Theorem 4.1.

Let the direct estimates of a survey follow the model given in Definition 3.1. Moreover, all matrices in the statements have earlier been defined in this section.

1. ${\widehat{\mathbf{\Psi }}}_1=\mathit{H}{({\mathit{H}}^{\prime }\mathit{H})}^{-1}{\tilde{\mathit{M}}}_1;$

2. If $\mathcal{C}(\mathit{A})\cap \mathcal{C}(\mathbf{\Psi })=\left\{\mathbf{0}\right\},$ $\widehat{\mathbf{\Psi }}{\tilde{\mathit{D}}}_3={\widehat{\mathbf{\Psi }}}_1{\widehat{\mathbf{\Psi }}}_1^{\prime }{\mathit{T}}_1^{\prime }{\mathit{S}}_2^{-1}{\mathit{X}}_1{\mathit{P}}_{{\tilde{D}}_3^{\prime }};$

   If additionally $r({\tilde{\mathit{D}}}_3^{\prime })=t,$ $\widehat{\mathbf{\Psi }}={\widehat{\mathbf{\Psi }}}_1{\widehat{\mathbf{\Psi }}}_1^{\prime }{\mathit{T}}_1^{\prime }{\mathit{S}}_2^{-1}{\mathit{X}}_1{\tilde{\mathit{D}}}_3^{\prime }{({\tilde{\mathit{D}}}_3{\tilde{\mathit{D}}}_3^{\prime })}^{-1};$

3. If $\mathcal{C}(\mathit{A})\cap \mathcal{C}(\mathbf{\Psi })=\left\{\mathbf{0}\right\},$ $\mathit{A}\widehat{\mathit{B}}{\tilde{\mathit{D}}}_1=\mathit{A}{({\mathit{A}}^{\prime }{\mathit{S}}_1^{-1}\mathit{A})}^{-}\mathit{A}\prime {\mathit{S}}_1^{-1}({\mathit{X}}_1{\mathit{P}}_{{\tilde{D}}_1^{\prime }}-\widehat{\mathbf{\Psi }}{\tilde{\mathit{D}}}_3);$

   If additionally $r({\tilde{\mathit{D}}}_1)=k$ and $r(\mathit{A})=q$ $$\widehat{\mathit{B}}={\left({\mathit{A}}^{\prime }{\mathit{S}}_1^{-1}\mathit{A}\right)}^{-1}\mathit{A}\prime {\mathit{S}}_1^{-1}({\mathit{Y}}_1{\tilde{\mathit{D}}}_1^{\prime }{\left({\tilde{\mathit{D}}}_1{\tilde{\mathit{D}}}_1^{\prime }\right)}^{-1}-\widehat{\mathbf{\Psi }}{\tilde{\mathit{D}}}_3{\tilde{\mathit{D}}}_1^{\prime }{\left({\tilde{\mathit{D}}}_1{\tilde{\mathit{D}}}_1^{\prime }\right)}^{-1});$$

4. $\widehat{\mathit{B}}{\mathit{D}}_2=\mathit{X}{\mathit{P}}_{D_2^{\prime }}-\mathit{A}{\widehat{\mathit{B}}}_1{\mathit{D}}_1{\mathit{P}}_{D_2^{\prime }}-\widehat{\mathbf{\Psi }}{\mathit{D}}_3{\mathit{P}}_{D_2^{\prime }};$

   If additionally $r({\mathit{D}}_2)=k_2$ $${\widehat{\mathit{B}}}_1=\mathit{X}{\mathit{D}}_2^{\prime }{\left({\mathit{D}}_2{\mathit{D}}_2^{\prime }\right)}^{-1}-\mathit{A}\widehat{\mathit{B}}{\mathit{D}}_1{\mathit{D}}_2^{\prime }{\left({\mathit{D}}_2{\mathit{D}}_2^{\prime }\right)}^{-1}-\widehat{\mathbf{\Psi }}{\mathit{D}}_3{\mathit{D}}_2^{\prime }{\left({\mathit{D}}_2{\mathit{D}}_2^{\prime }\right)}^{-1};$$

5. $(n-l){\widehat{\mathbf{\Sigma }}}_e={\mathit{S}}_2+{\widehat{\mathit{T}}}_2{\mathit{T}}_1{\mathit{Y}}_1{\mathit{P}}_{{\tilde{D}}_3^{\prime }}{\mathit{Y}}_1^{\prime }{\mathit{T}}_1^{\prime }{\widehat{\mathit{T}}}_2^{\prime },$ where $${\widehat{\mathit{T}}}_2={\mathit{I}}_p-{\mathit{T}}_1{\widehat{\mathbf{\Psi }}}_1{\left({\widehat{\mathbf{\Psi }}}_1^{\prime }{\mathit{T}}_1^{\prime }{\mathit{S}}_2^{-1}{\mathit{T}}_1{\widehat{\mathbf{\Psi }}}_1\right)}^{-}{\widehat{\mathbf{\Psi }}}_1^{\prime }{\mathit{T}}_1^{\prime }{\mathit{S}}_2^{-1};$$

6. Suppose $\frac{1}{l}{\mathit{Y}}_0{\mathit{Y}}_0^{\prime }-{\widehat{\mathbf{\Sigma }}}_e$ is positive definite. Then $${\widehat{\mathbf{\Sigma }}}_u=\frac{1}{l}{\mathit{Y}}_0{\mathit{Y}}_0^{\prime }-{\widehat{\mathbf{\Sigma }}}_e.$$

5. Prediction

Small area estimation is very often about prediction of unobserved units (subareas). For our model, due to normality assumptions, it is relatively straight forward to construct predictors, in particular if the explicit estimators of the previous section are utilized. The strategy will be to first predict U which will be based on the observed data and thereafter the unobserved units are predicted.

For prediction of U the conditional mean $E[\mathit{U}|\mathit{Y}]$ is crucial and thus the joint distribution of U and Y will be derived. According to Definition 3.1, U and E are independently distributed $(\text{vec}\mathit{U},\text{vec}\mathit{Y})$ is normally distributed, i.e. $$\left(\begin{array}{c}\text{vec}\mathit{U}\\ \text{vec}\mathit{Y}\end{array}\right)\sim N_{p(n+l)}\left(E\right[•],D[•\left]\right),$$ where $$\begin{array}{c}E[•]=(\mathbf{0},\text{vec}\prime (\mathbf{\text{ABC}}+{\mathit{B}}_1{\mathit{C}}_1+\mathbf{\Psi }\mathit{F}\left)\right)\prime ,\\ D[•]=\left(\begin{array}{cc}\mathit{I}\otimes {\mathbf{\Sigma }}_u& \mathit{Z}\otimes {\mathbf{\Sigma }}_u\\ \mathit{Z}\prime \otimes {\mathbf{\Sigma }}_u& \mathit{Z}\prime \mathit{Z}\otimes {\mathbf{\Sigma }}_u+\mathit{V}\otimes {\mathbf{\Sigma }}_e\end{array}\right).\end{array}$$

Thus, $$\begin{array}{c}E\left[\text{vec}\mathit{U}\right|\text{vec}\mathit{Y}]=(\mathit{Z}\otimes {\mathbf{\Sigma }}_u\left){({\mathit{Z}}^{\prime }\mathit{Z}\otimes {\mathbf{\Sigma }}_u+\mathit{V}\otimes {\mathbf{\Sigma }}_e)}^{-1}\text{vec}\right(\mathit{Y}-\mathbf{\text{ABC}}-{\mathit{B}}_1{\mathit{C}}_1-\mathbf{\Psi }\mathit{F})\end{array}$$ and the next proposition can be stated.

Proposition 5.1.

Let the model be given in Definition 3.1 and use the notation introduced earlier in Section 4.

1. The predicted value $\widehat{\mathit{U}}$ is given by $$\text{vec}\widehat{\mathit{U}}=(\mathit{Z}\otimes {\widehat{\mathbf{\Sigma }}}_u){({\mathit{Z}}^{\prime }\mathit{Z}\otimes {\widehat{\mathbf{\Sigma }}}_u+\mathit{V}\otimes {\widehat{\mathbf{\Sigma }}}_e)}^{-1}\text{vec}(\mathit{Y}-\mathit{A}\widehat{\mathit{B}}\mathit{C}-{\widehat{\mathit{B}}}_1{\mathit{C}}_1-\widehat{\mathbf{\Psi }}\mathit{F}),$$ where ${\widehat{\mathbf{\Sigma }}}_u,{\widehat{\mathbf{\Sigma }}}_e,$ $\mathit{A}\widehat{\mathit{B}}\mathit{C},$ ${\widehat{\mathit{B}}}_1{\mathit{C}}_1$ and $\widehat{\mathbf{\Psi }}\mathit{F}$ are presented in Theorem 4.1 and Y consists of subarea estimates.

2. Let ${\mathit{Y}}_o$ consist of the unobserved direct estimates, corresponding to the non-sampled units. The corresponding covariables are available which are denoted ${\mathit{C}}_o,{\mathit{C}}_{1o},{\mathit{F}}_o$ and ${\mathit{Z}}_o.$ Then $$\text{vec}{\widehat{\mathit{Y}}}_o=\text{vec}(\mathit{A}\widehat{\mathit{B}}{\mathit{C}}_o+{\widehat{\mathit{B}}}_1{\mathit{C}}_{1o}+\widehat{\mathbf{\Psi }}{\mathit{F}}_o)+({\mathit{Z}}_o^{\prime }\otimes \mathit{I})\text{vec}\widehat{\mathit{U}},$$ where $\text{vec}\widehat{\mathit{U}}$ is given in (i).

Additional information

Funding

This research has been supported by The Swedish Foundation for Humanities and Social Sciences (P14-0641:1) and The Swedish Natural Research Council (2017-03003).