Gram–Schmidt–Fisher scoring algorithm for parameter orthogonalization in MLE

Abstract

The estimation of parameters is a key component in statistical modelling and inference. However, parametrization of certain likelihood functions could lead to highly correlated estimates, causing numerical problems, mathematical complexities and difficulty in estimation or erroneous interpretation and subsequently inference. In statistical estimation, the concept of orthogonalization is familiar as a simplifying technique that allows parameters to be estimated independently and thus separates information from each other. We introduce a Fisher scoring iterative process that incorporates the Gram–Schmidt orthogonalization technique for maximum likelihood estimation. A finite mixture model for correlated binary data is used to illustrate the implementation of the method with discussion of application to oesophageal cancer data.

Keywords:

• correlated data
• Fisher scoring algorithm
• Gram–Schmidt orthogonalization
• maximum likelihood estimation
• orthogonal parameters

Public Interest Statement

In statistical estimation, the concept of orthogonalization is simply a technique that allows parameters to be estimated independently, through the rotation of the parameter space. The paper introduces an iterative algorithm that utilizes the Fisher scoring process and incorporates the Gram–Schmidt orthogonalization technique for maximum likelihood estimation. The algorithm that we propose provides an exceedingly convenient method for multi-parameter statistical problems first to reduce numerical difficulties and second to improve accuracy in parameter estimation.

1. Introduction

The problem of estimating parameters is one of the key stages in fitting a statistical model to a set of data. Typically, the model will contain some parameters which are not of interest in themselves, but whose values affect the inferences that we make about the parameters of direct interest. Parametrization used when working with certain distributions could lead to high correlations of the maximum likelihood estimates. Maximization algorithms converge rapidly if the initial estimates are good, the likelihood function is well approximated by a quadratic in the neighbourhood of the parameter space and the information matrix is well conditioned, which means that the parameter estimates are not strongly intercorrelated. High intercorrelation of parameters causes numerical problems, and difficulty or erroneous interpretation in the parameter estimation. In statistical estimation, the concept of orthogonalization is familiar as a simplifying technique that allows parameters to be estimated independently, through rotation of the parameter space. Computationally, all that is required is that the original parameters should be computable from the new ones, and vice-versa. An injective transformation will therefore be desirable, though not necessary.

Ross (ross) gave a comprehensive discussion of techniques that can be used to reduce correlation in particular problems. Among them are sequential and nested maximizations. Kalbliesch and Sprott (kalbliesch) discuss methods aimed at eliminating nuisance parameters from the likelihood function so that inference can be made about only the parameters of interest. Amari (amari, amari1) derived from a theoretical point of view, a general criteria for the existence and construction of orthogonal parameters for statistical models using differential geometry. Cox and Reid (cox) gave a more general procedure and covered the advantages in maximum likelihood estimation. The approach, though popular in theoretical statistics, is computationally impractical in many application situations. Willmot (willmot) discussed orthogonal parametrization but only for a two-parameter family of discrete distributions. His method, however, does not allow for higher parametric problems. Hurlimann (hurlimann) discusses the existence of orthogonal parameterization to the mean, characterized by a partial differential equation involving the mean, the variance and cummulant generating functions. Godambe (godambe) deals with the problem of nuisance parameters, within a semi-parametric set-up which includes the class of distributions associated with generalized linear models. Their technique uses the optimum orthogonal estimating functions of Godambe and Thompson (godambe1). Bonney and Kissling (bonney) applied the Gram–Schmidt orthogonalization (Clayton, clayton) to multi-normal variates and presented an application in genetics. Kwagyan, Apprey and Bonney (kwagyan1) presented the idea of Gram–Schmidt parameter orthogonalization in genetic analysis. A major consequence of parameter orthogonalization is that the maximum likelihood estimate of one parameter is not affected or changes slowly with the misspecification of the other. The most direct interpreted consequence is that the maximum likelihood estimates of orthogonal coordinates are asymptotically independent.

The remainder of this paper is organized as follows. In Section 2, the disposition model for clustered binary data proposed by Bonney (bonney1, bonney2) and further investigated by Kwagyan (kwagyan) is introduced and used as a motivation. In Section 3, an overview of maximum likelihood estimation and parameter orthogonalization is then presented. Section 4 develops the proposed Gram–Schmidt parameter orthogonalization and iterative scheme for maximum likelihood estimation. In section Section , fitting algorithm of the disposition model is discussed with application to oesophageal cancer data in Chinese families. Finally, Section wraps up with discussion of the methodology and outlines directions for further research.

2. Model for correlated binary data

Consider a sample of N groups, each of size $n_i,$$i=1,\dots ,N;$ and let $Y_i={(Y_{i1},\dots ,Y_{{in}_i})}^T$ denote the vector of binary outcomes for the ith group. It is postulated that the measures of the outcomes within a group are possibly correlated. Further suppose the jth subject in the ith group has $1\times p$ individual-specific covariates $X_{ij}=(x_{ij1},\cdots ,x_{ijp})$ and let the ith group have group-specific covariates $Z_i=(z_1,\dots ,z_q)$. Let ${\mathit{\delta }}_{ij}$ denote the conditional probability of $Y_{ij}=1$ given $Y_{i{j}^{\prime }}=1$; that is$$\begin{array}{c}\hfill {\mathit{\delta }}_{ij}=P(Y_{ij}=1|Y_{i{j}^{\prime }}=1),j\ne j^{\prime },j,j^{\prime }=1,2,\dots ,n.\end{array}$$

We call ${\mathit{\delta }}_{ij}$ the individual disposition which is simply interpreted as the probability of the outcome on one unit given another unit from the same cluster has the attribute. This in essence is an indication of aggregation. Assume further that within the group, a pair of observations satisfy the relation$$\frac{P(Y_{ij}=1,Y_{i{j}^{\prime }}=1)}{P(Y_{ij}=1)P(Y_{i{j}^{\prime }}=1)}=\frac{1}{{\mathit{\alpha }}_i},{\mathit{\alpha }}_i>0,j\ne j^{\prime },j,j^{\prime }=1,2,\cdots ,n.$$

where ${\mathit{\alpha }}_i$ is common for all pairs. Clearly, ${\mathit{\alpha }}_i=1$ implies independence of the observations. Thus, ${\mathit{\alpha }}_i$ is a measure of the departure from independence and is called the relative disposition. With the above definitions, Bonney (bonney1, bonney2), and from further investigation through a latent mixture formulation by Kwagyan (kwagyan), has shown that the joint distribution for the N groups can be based on(1) $$\begin{array}{c}\hfill L(\mathrm{\Theta })=\prod _{i=1}^N\left\{(1-{\mathit{\alpha }}_i)\prod _{j=1}^{n_i}(1-y_{ij})+{\mathit{\alpha }}_i\prod _{j=1}^{n_i}{\mathit{\delta }}_{ij}^{yij}{\left(1-{\mathit{\delta }}_{ij}\right)}^{1-y_{ij}}\right\}\end{array}$$(1)

The ${\mathit{\alpha }}_i$ and ${\mathit{\delta }}_{ij}$ are generally modelled as$${\mathit{\alpha }}_i=\frac{1+e^{[-\{M(Z_i)+D(Z_i)\}]}}{1+e^{-M(Z_i)}},{\mathit{\delta }}_{ij}=\frac{1}{1+e^{-[M(Z_i)+D(Z_i)+W(X_{ij})]}}$$

where $M(Z_i)$ is function of mean effects, $D(Z_i)$ is a measure of within-group dependence and $W(X_{ij})$ a function describing the effect of the individual covariates. Typically, $M(Z_i),D(Z_i)$ and $W(X_{ij})$ are modelled, respectively, as$$\begin{array}{cc}\hfill M(Z_i)& ={\mathit{\lambda }}_0+{\mathit{\lambda }}_1{Z}_1+\dots +{\mathit{\lambda }}_p{Z}_p\hfill \\ \hfill D(Z_i)& ={\mathit{\gamma }}_0+{\mathit{\gamma }}_1{Z}_1+\dots +{\mathit{\gamma }}_q{Z}_q\hfill \\ \hfill W(X_{ij})& ={\mathit{\beta }}_1{X}_1+{\mathit{\beta }}_2{X}_2+\dots +{\mathit{\beta }}_r{X}_r\hfill \end{array}$$

where$${\mathrm{\Theta }}^{\prime }=\{\mathrm{\Lambda },\mathrm{\Gamma },\mathit{\beta }\}$$

and$$\mathrm{\Lambda }={\mathit{\gamma }}_0,\dots ,{\mathit{\gamma }}_q,\mathrm{\Gamma }={\mathit{\lambda }}_0,\dots ,{\mathit{\lambda }}_q\},\mathit{\beta }=\{{\mathit{\beta }}_1,\dots ,{\mathit{\beta }}_r\}$$

are unknown parameters to be estimated. When data consist of a few large clusters or truncated or size-biased samples, Bonney (bonney3) has shown that there is little or no information to separate the effects of M(Z) and D(Z). The estimates of the parameters tend to be highly correlated. Similar problems occur in other approaches as well. The popular estimating equations approach, GEE (Liang & Zeger, liang), suffers the same problems when there are only a few large clusters (see discussion of Prentice, prentice).

3. Maximum likelihood estimation & parameter orthogonalization

There are different statistical methods for estimating parameters, but the approach most commonly used is that based on maximum likelihood estimation. Maximum likelihood estimates of parameters are those values of the parameters that make the likelihood function a maximum. Let $\mathbf{y}=(y_1,y_2,\dots ,y_n)$ be observed data from a population with a probability distribution $P(y;\mathit{\theta })$ indexed by the vector of unknown parameters $\mathit{\theta }=({\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_p)$. The likelihood function for $\mathit{\theta }$ is the joint distribution$$L(\mathit{\theta };y)=\prod _{i=1}^{n}P(y_i;\mathit{\theta })$$

According to the likelihood principle, $\widehat{\mathit{\theta }}$ is regarded as the value of $\mathit{\theta }$ which is most strongly supported by the observed data. In practice, $\widehat{\mathit{\theta }}$ is obtained by direct maximization of Log $L(\mathit{\theta })$ or by solving the set of equations from the score function, $U(\mathit{\theta }),$where$$U(\mathit{\theta })=\frac{\mathit{\partial }\text{log}L}{\mathit{\partial }\mathit{\theta }}=0$$

for which the Hessian matrix $H(\mathit{\theta })=(\frac{{\mathit{\partial }}^2\text{log}L}{\mathit{\partial }\mathit{\theta }\mathit{\partial }{\mathit{\theta }}^T})$ is negative definite. The matrix $I(\mathit{\theta })=E[-H(\mathit{\theta })]=E(-\frac{{\mathit{\partial }}^{2}logL}{\mathit{\partial }\mathit{\theta }\mathit{\partial }{\mathit{\theta }}^T})$ is called Fisher (expected) Information matrix for $\mathit{\theta }$ and its inverse, $\mathrm{\Omega }(\mathit{\theta })={[I(\mathit{\theta })]}^{-1}$ gives the asymptotic variances and covariances of the maximum likelihood estimates.

Suppose $\mathit{\theta }$ is a p parameter vector that is partitioned into two vectors ${\mathit{\theta }}_1$ and ${\mathit{\theta }}_2$ of lengths $p_1$ and $p_2$, respectively, $p_1+p_2=p$ . Cox and Reid (cox) define ${\mathit{\theta }}_1$ to be orthogonal to ${\mathit{\theta }}_2$ if the elements of the information matrix satisfy$$E\left(\frac{\mathit{\partial }l}{\mathit{\partial }{\mathit{\theta }}_s}\frac{\mathit{\partial }l}{\mathit{\partial }{\mathit{\theta }}_t}\right)=E\left(-{\displaystyle \frac{{\mathit{\partial }}^{2}l}{\mathit{\partial }{\mathit{\theta }}_s\mathit{\partial }{\mathit{\theta }}_t}}\right)=0,s=1,\dots ,p_1;t=p_1+1,\dots ,p$$

If this holds for all $\mathit{\theta }$ in the parameter space, it is sometimes called global orthogonality. If it holds at only one parameter value say ${\mathit{\theta }}^0$, then ${\mathit{\theta }}_1$ and ${\mathit{\theta }}_2$ are said to be locally orthogonal at ${\mathit{\theta }}^0$.

We now discuss Cox and Reid (cox) approach for the construction of orthogonal parameters. Suppose that initially the likelihood is specified in terms of $(\mathit{\theta },\mathit{\gamma })$ , $\mathit{\gamma }=({\mathit{\gamma }}_1,\dots ,{\mathit{\gamma }}_p)$ and let $l(\mathit{\theta },\mathit{\gamma })$ be the log-likelihood function. Assume $\mathit{\theta }$ is the parameter of interest and $\mathit{\gamma }$ the set of nuisance parameters. We seek a transformation from $(\mathit{\theta },\mathit{\gamma })$ to a new set of orthogonal parameters $(\mathit{\theta },\mathit{\lambda })$, ${\mathit{\lambda }=(\mathit{\lambda }}_1,\dots ,{\mathit{\lambda }}_p)$. It is easiest to think of the original parameters $\mathit{\gamma }$ as some function of the new parameters, $\mathit{\lambda }$ , that is$$\mathit{\gamma }=\mathit{\gamma }(\mathit{\theta },\mathit{\lambda })$$

Then, the log-likelihood function can be expressed as$$l^{\ast }(\mathit{\theta },\mathit{\lambda })=l(\mathit{\theta },{\mathit{\gamma }}_1(\mathit{\theta },\mathit{\lambda }),\dots ,{\mathit{\gamma }}_p(\mathit{\theta },\mathit{\lambda })),$$

where now $l^{\ast }$ refers to the log-likelihood in the new parametrization. Taking derivatives of this equation with respect to the new parameters, we have by the chain rule;$$\frac{\mathit{\partial }{l}^{\ast }}{\mathit{\partial }\mathit{\theta }}=\frac{\mathit{\partial }l}{\mathit{\partial }\mathit{\theta }}+\sum _{j=1}^p\frac{\mathit{\partial }l}{\mathit{\partial }{\mathit{\gamma }}_j}\frac{\mathit{\partial }{\mathit{\gamma }}_j}{\mathit{\partial }\mathit{\theta }}$$

and(2) $$\begin{array}{c}\hfill \frac{{\mathit{\partial }}^2{l}^{\ast }}{\mathit{\partial }\mathit{\theta }\mathit{\partial }{\mathit{\lambda }}_k}=\sum _l^p{\displaystyle \frac{{\mathit{\partial }}^{2}l}{\mathit{\partial }\mathit{\theta }\mathit{\partial }{\mathit{\gamma }}_l}}{\displaystyle \frac{\mathit{\partial }{\mathit{\gamma }}_l}{\mathit{\partial }{\mathit{\lambda }}_k}}+\sum _j^p\sum _l^p{\displaystyle \frac{{\mathit{\partial }}^{2}l}{\mathit{\partial }{\mathit{\gamma }}_j\mathit{\partial }{\mathit{\gamma }}_l}}{\displaystyle \frac{\mathit{\partial }{\mathit{\gamma }}_l}{\mathit{\partial }{\mathit{\theta }}_k}}{\displaystyle \frac{\mathit{\partial }{\mathit{\gamma }}_j}{\mathit{\partial }\mathit{\theta }}}+\sum _l^p{\displaystyle \frac{\mathit{\partial }l}{\mathit{\partial }{\mathit{\gamma }}_l}}{\displaystyle \frac{{\mathit{\partial }}^2{\mathit{\gamma }}_l}{\mathit{\partial }\mathit{\theta }\mathit{\partial }{\mathit{\lambda }}_k}}\end{array}$$(2)

The derivatives of $\mathit{\gamma }$ are not functions of the data, and hence are constant with respect to expectation. Therefore, after taken expectations with respect to the distribution of the data indexed by the parameters $\{\mathit{\theta },\mathit{\gamma }\}$ ,the last term in Equation 3.1 is zero. If the parameters are orthogonal, then again from Equation 3.1,$$\sum _{l=1}^p\frac{\mathit{\partial }{\mathit{\gamma }}_l}{\mathit{\partial }{\mathit{\lambda }}_k}\left(E\left\{{\displaystyle \frac{{\mathit{\partial }}^{2}l}{\mathit{\partial }\mathit{\theta }\mathit{\partial }{\mathit{\gamma }}_l}}\right\}+\sum _{j=1}^{p}E\left\{\frac{{\mathit{\partial }}^{2}l}{\mathit{\partial }{\mathit{\gamma }}_l\mathit{\partial }{\mathit{\gamma }}_j}\right\}\frac{\mathit{\partial }{\mathit{\gamma }}_j}{\mathit{\partial }\mathit{\theta }}\right)=0.$$

so that the orthogonality equations are$$\sum _{l=1}^p{\displaystyle \frac{\mathit{\partial }{\mathit{\gamma }}_l}{\mathit{\partial }{\mathit{\lambda }}_k}}\left(E\left\{\frac{{\mathit{\partial }}^{2}l}{\mathit{\partial }\mathit{\theta }\mathit{\partial }{\mathit{\gamma }}_l}\right\}+\sum _{j=1}^{p}E\left\{\frac{{\mathit{\partial }}^{2}l}{\mathit{\partial }{\mathit{\gamma }}_l\mathit{\partial }{\mathit{\gamma }}_j}\right\}{\displaystyle \frac{\mathit{\partial }{\mathit{\gamma }}_j}{\mathit{\partial }\mathit{\theta }}}\right)=0,k=1,\cdots ,p.$$

We require the transformation from $(\mathit{\theta },\mathit{\gamma })$ to $(\mathit{\theta },\mathit{\lambda })$ to have a non-zero Jacobian; hence,$$\sum _{j=1}^{p}E\left\{\frac{{\mathit{\partial }}^{2}l}{\mathit{\partial }{\mathit{\gamma }}_l\mathit{\partial }{\mathit{\gamma }}_k}\right\}{\displaystyle \frac{\mathit{\partial }{\mathit{\gamma }}_j}{\mathit{\partial }\mathit{\theta }}}=-E\left\{{\displaystyle \frac{{\mathit{\partial }}^{2}l}{\mathit{\partial }\mathit{\theta }\mathit{\partial }{\mathit{\gamma }}_k}}\right\},k=1,\dots ,p$$

This is a system of p differential equations which must be solved for $\mathit{\gamma }(\mathit{\theta },\mathit{\lambda })$. In fact, since $\mathit{\lambda }$ does not enter explicitly into the equations, the solution for ${\mathit{\gamma }}_j$ can contain an arbitrary function of $\mathit{\lambda }$ as the integrating constant. It is noted in the discussion of Cox and Reid that although it is sometimes theoretically possible to find a differential equation, simple explicit solutions of the differential equation were not feasible for the some models. There are also cases where explicit solutions are possible, but the original nuisance parameters could not be explicitly expressed in terms of the orthogonal parameter (Hills, hills). In general, global orthogonalization can also not be achieved by this approach.

4. Gram–Schmidt parameter orthogonalization

The Gram–Schmidt orthogonalization process (Clayton, clayton) is equivalent to the linear transformation $\mathrm{\Theta }=({\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_p)$ to ${\mathrm{\Theta }}^{\ast }=({\mathit{\theta }}_1^{\ast },\dots ,{\mathit{\theta }}_p^{\ast })$ defined as:(3) $$\begin{array}{cc}\hfill {\mathit{\theta }}_1^{\ast }& ={\mathit{\theta }}_1\hfill \\ \hfill {\mathit{\theta }}_2^{\ast }& ={\mathit{\theta }}_2-b_{21}{\mathit{\theta }}_1\hfill \\ \hfill & -------\hfill \\ \hfill {\mathit{\theta }}_j^{\ast }& ={\mathit{\theta }}_j-\sum _{k=1}^{j-1}{b}_{j(j-k)}{\mathit{\theta }}_{j-k},j=2,\dots ,p\hfill \end{array}$$(3)

In linear regression set-up, ${\mathit{\theta }}_j^{\ast }$ is ${\mathit{\theta }}_j$ adjusted for ${\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_p,$ and $b_{jk},$$j=2,\dots ,p;$$k=1,\dots ,(j-1)$ are the multiple regression coefficients. Writing this using matrix notation, we have$${\mathrm{\Theta }}^{\ast }=\mathbf{B}\mathrm{\Theta }$$

where $\mathbf{B}$, the transformation matrix, is lower triangular with ones along the diagonal. The transformation matrix is chosen so that ${\mathit{\theta }}_1^{\ast },{\mathit{\theta }}_2^{\ast },\dots ,{\mathit{\theta }}_p^{\ast }$ are mutually uncorrelated. The Jacobian of the transformation is unity. Suppose $\mathrm{\Omega }$ is the covariance matrix of $\mathrm{\Theta }$, then the covariance matrix of${\mathrm{\Theta }}^{\ast },$${\mathrm{\Omega }}^{\ast }=B\mathrm{\Omega }{B}^T$ is a diagonal matrix. We recommend that the parameters be ordered in terms of interest. This ensures that the parameters of most interest are least affected by round-off errors.

4.1. Evaluation of the transformation matrix

To evaluate the transformation matrix, $\mathbf{B}$, let $\mathrm{\Sigma }=(c_{ij}),$$i,j=1,\dots ,p,$ where $c_{ij}$$=cov({\mathit{\theta }}_i,{\mathit{\theta }}_j)={\left[E(-\frac{{\mathit{\partial }}^{2}logL}{\mathit{\partial }{\mathit{\theta }}_i\mathit{\partial }{\mathit{\theta }}_j})\right]}^{-1},$ are elements of the covariance matrix. The orthogonality relation, $c_{ij}^{\ast }=cov({\mathit{\theta }}_i^{\ast },{\mathit{\theta }}_j^{\ast })=0$$(i\ne j)$ implies that for $i<j$,$$\begin{array}{cc}\hfill 0& =c_{ij}^{\ast }=cov({\mathit{\theta }}_i^{\ast },{\mathit{\theta }}_j^{\ast })\hfill \\ \hfill & =cov({\mathit{\theta }}_i-\sum _{k^{\prime }=1}^{i-1}{b}_{i(i-k^{\prime })}{\mathit{\theta }}_{i-k^{\prime }},{\mathit{\theta }}_j-\sum _{k=1}^{j-1}{b}_{j(j-k)}{\mathit{\theta }}_{j-k})\hfill \\ \hfill & =cov({\mathit{\theta }}_i,{\mathit{\theta }}_j)-\sum _{k=1}^{j-1}{b}_{j(j-k)}cov({\mathit{\theta }}_i,{\mathit{\theta }}_{j-k})\hfill \\ \hfill & -\sum _{k^{\prime }=1}^{i-1}{b}_{i(i-k^{\prime })}\left\{cov({\mathit{\theta }}_{j-k^{\prime }},{\mathit{\theta }}_j)+\sum _{k=1}^{j-1}{b}_{j(j-k)}cov({\mathit{\theta }}_{i-k^{\prime }},{\mathit{\theta }}_{j-k})\right\}\hfill \\ \hfill & =c_{ij}-\sum _{k=1}^{j-1}{b}_{j(j-k)}{c}_{i,j-k}-\sum _{k^{\prime }=1}^{i-1}{b}_{i(i-k^{\prime })}\left\{c_{j-k^{\prime },j}+\sum _{k=1}^{j-1}{b}_{j(j-k)}{c}_{i-k^{\prime },j-k}\right\}\hfill \\ \hfill & =c_{ij}-\sum _{k=1}^{j-1}{b}_{j(j-k)}{c}_{i,j-k}-\sum _{k^{\prime }=1}^{i-1}{b}_{i(i-k^{\prime })}\left\{c_{j-k^{\prime },j}^{\ast }\right\}\hfill \\ \hfill & =c_{ij}-\sum _{k=1}^{j-1}{b}_{j(j-k)}{c}_{i,j-k},\text{since}{c}_{j-k^{\prime },j}^{\ast }=0\hfill \end{array}$$

Thus, the system of linear equations to be solved for entries of the transformation matrix based on the covariance matrix is:(4) $$\begin{array}{c}\hfill c_{ij}=\sum _{k=1}^{j-1}{b}_{j(j-k)}{c}_{i(j-k)},i=1,\dots (j-1),j=2,\dots ,p.\end{array}$$(4)

Consequently, solving the system of linear equations (4.2), we obtain the elements of the transformation matrix as(5) $$\begin{array}{c}\hfill b_{jk}=\sum _{r=1}^{j-1}{c}_{rj}{I}_{r(j-k)},j=2,\dots ,p;k=1,\dots ,(j-1).\end{array}$$(5)

where $I_{r(j-k)}$ are entries of the information matrix.The observed information matrix would be a good approximation of the expected information, if there is difficulty evaluating it.

For illustration, when $p=4$, we have Gram–Schmidt orthogonalization process in matrix notation as,$$\left(\begin{array}{c}{\mathit{\theta }}_1^{\ast }\\ {\mathit{\theta }}_2^{\ast }\\ {\mathit{\theta }}_3^{\ast }\\ {\mathit{\theta }}_4^{\ast }\end{array}\right)=\left(\begin{array}{cccc}1& 0& 0& 0\\ -b_{21}& 1& 0& 0\\ -b_{31}& -b_{32}& 1& 0\\ -b_{41}& -b_{42}& -b_{43}& 1\end{array}\right)\left(\begin{array}{c}{\mathit{\theta }}_1\\ {\mathit{\theta }}_2\\ {\mathit{\theta }}_3\\ {\mathit{\theta }}_4\end{array}\right)$$

The orthogonalization relationship of the parameters $({\mathit{\theta }}_1^{\ast },{\mathit{\theta }}_2^{\ast },{\mathit{\theta }}_3^{\ast },{\mathit{\theta }}_4^{\ast })$ implies$$\begin{array}{cc}\hfill 0& =c_{12}-b_{21}{c}_{11}\hfill \\ \hfill 0& =c_{13}-b_{32}{c}_{12}-b_{31}{c}_{11}\hfill \\ \hfill 0& =c_{23}-b_{32}{c}_{22}-b_{31}{c}_{12}\hfill \\ \hfill 0& =c_{14}-b_{43}{c}_{13}-b_{42}{c}_{12}-b_{41}{c}_{11}\hfill \\ \hfill 0& =c_{24}-b_{43}{c}_{23}-b_{42}{c}_{22}-b_{41}{c}_{12}\hfill \\ \hfill 0& =c_{34}-b_{43}{c}_{33}-b_{42}{c}_{23}-b_{41}{c}_{13}\hfill \end{array}$$

And so the system of linear equations to be solved for the elements of the transformation matrix is:

Let Q be the matrix of coefficients for the system of equations (4.4) Above; then, we note that Q is a patterned block diagonal matrix. Furthermore, let $D_r(r=1,2,3),$ be the block diagonal of Q; then, in this illustration,$$Q=\left(\begin{array}{ccc}{D}_1& 0& 0\\ 0& D_2& 0\\ 0& 0& D_3\end{array}\right)$$

where$$D_1=(c_{11}),D_2=\left(\begin{array}{cc}{c}_{11}& c_{12}\\ c_{12}& c_{22}\end{array}\right),D_3=\left(\begin{array}{ccc}{c}_{11}& c_{12}& c_{13}\\ c_{12}& c_{22}& c_{23}\\ c_{13}& c_{23}& c_{33}\end{array}\right)$$

Thus, $D_r$ is the covariance matrix of the first r parameters. It can easily be shown that $D_r$ is symmetric and positive definite and so Q is non-singular. A unique solution thus exists for the system of linear Equation (4.3) and in general for Equation (4.2). And so, the elements of the transformation matrix, B, can be obtained as

or

where $I_{ij}$ are the entries of the information matrix .

Having obtained ${\mathrm{\Theta }}^{\ast }=({\mathit{\theta }}_1^{\ast },\dots ,{\mathit{\theta }}_p^{\ast })$ the original parameters, $\mathrm{\Theta }=({\mathit{\theta }}_1,\dots ,{\mathit{\theta }}_p)$ can be obtained recursively as$$\begin{array}{cc}\hfill {\mathit{\theta }}_1& ={\mathit{\theta }}_1^{\ast }\hfill \\ \hfill {\mathit{\theta }}_j& ={\mathit{\theta }}_j^{\ast }+\sum _{k=1}^{j-1}{b}_{j(j-k)}{\mathit{\theta }}_{j-k},j=2,\dots ,p\hfill \end{array}$$

or writing this using matrix notation, we have$$\mathrm{\Theta }={\mathbf{B}}^{-1}{\mathrm{\Theta }}^{\ast }$$

4.2. Block orthogonalization

Let ${\mathrm{\Theta }}^T=\{\mathrm{\Lambda },\mathrm{\Gamma },\mathit{\beta }\}$ be a set of parameters to be estimated where$$\mathrm{\Lambda }={({\mathit{\lambda }}_0,{\mathit{\lambda }}_1,\dots ,{\mathit{\lambda }}_q)}^T;\mathrm{\Gamma }={({\mathit{\gamma }}_0,{\mathit{\gamma }}_1,\dots ,{\mathit{\gamma }}_p)}^T;\mathit{\beta }=({\mathit{\beta }}_1,{\mathit{\beta }}_{1,}\dots ,{\mathit{\beta }}_r{)}^T$$

Suppose further that the set ${\mathrm{\Theta }}^T=\{\mathrm{\Lambda },\mathrm{\Gamma },\mathit{\beta }\}$ is correlated. Then, we wish to find a new set${\mathrm{\Theta }}^{\ast T}=\{{\mathrm{\Lambda }}^{\ast },{\mathrm{\Gamma }}^{\ast },{\mathit{\beta }}^{\ast }\}$ through a linear transformation such that the vector of parameters in ${\mathrm{\Theta }}^{\ast }$ is mutually uncorrelated. We allow for correlation, if any, within each set of parameters in ${\mathrm{\Omega }}^{\ast }$.

Further suppose the vector $\mathit{\beta }$ is the set of parameters of interest. Then, the Gram–Schmidt orthogonalization process computes the new set of parameters in terms of the original parameters as follows:$$\begin{array}{cc}\hfill {\mathit{\beta }}^{\ast }& =\mathit{\beta }\hfill \\ \hfill {\mathrm{\Lambda }}^{\ast }& =\mathrm{\Lambda }-B_{21}\mathit{\beta }\hfill \\ \hfill {\mathrm{\Gamma }}^{\ast }& =\mathrm{\Gamma }-B_{31}\mathrm{\Lambda }-B_{32}\mathit{\beta }\hfill \end{array}$$

where $B_{21},$$B_{31}$ and $B_{32}$ are matrices with dimensions $q\times r$, $p\times q$ and $p\times r$, respectively.

In matrix notation, we write$$\left(\begin{array}{c}{\mathit{\beta }}^{\ast }\\ {\mathrm{\Lambda }}^{\ast }\\ {\mathrm{\Gamma }}^{\ast }\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ -B_{21}& 1& 0\\ -B_{31}& -B_{32}& 1\end{array}\right)\left(\begin{array}{c}\mathit{\beta }\\ \mathrm{\Lambda }\\ \mathrm{\Gamma }\end{array}\right)$$

or$${\mathrm{\Theta }}^{\ast }=\mathbf{B}\mathrm{\Theta }$$

where $\mathbf{B}$ is a lower triangular block matrix whose diagonal unit matrix is the transformation matrix chosen such that ${\mathrm{\Lambda }}^{\ast },{\mathrm{\Gamma }}^{\ast },{\mathit{\beta }}^{\ast }$ are mutually uncorrelated and where ${\mathrm{\Theta }}^{\ast T}\mathrm{\Theta }$ is a block diagonal matrix. The Jacobian of the transformation is unity. The only way in which this procedure could break down would be if one of the vectors of $\mathrm{\Theta }$ is identically zero. From the method of formation of ${\mathrm{\Theta }}^{\ast }$, it is clear that it is a linear combination of the vectors in $\mathrm{\Theta }$. If $\mathrm{\Theta }$ is identically zero, this means that $\mathrm{\Theta }$ are linearly dependent, thus contradicting the assumption of independence. Extensions to more than three sets of vectors of parameters readily follow.

4.3. Gram–Schmidt–Fisher scoring algorithm

We introduce a modification of the Fisher scoring algorithm incorporating the Gram–Schmidt process to obtain an information matrix which is diagonal and thus ensuring the approximate (near) orthogonality and consequently the stability of the estimates of the new parameters. Since the transformation is linear and injective, the original parameters are readily obtained.

Let $l(\mathit{\theta };\mathbf{y})$ be the log-likelihood in the original Parametrization; then, Fisher scoring algorithm is given by the iterative routine(6) $$\begin{array}{c}\hfill {\mathit{\theta }}_{m+1}={\mathit{\theta }}_m+{[{\mathbf{I}(\mathit{\theta }}_m)]}^{-1}{\mathbf{U}(\mathit{\theta }}_m)\end{array}$$(6)

where $\mathbf{I}(\mathit{\theta })$ is the expected information matrix and $\mathbf{U}(\mathit{\theta })$ is the score function.

Suppose $\mathit{\theta }$ is transformed to ${\mathit{\theta }}^{\ast }$ through the Gram–Schmidt orthogonalization process. Let $\mathbf{B}$$=(b_{ij})$ be the matrix of transformation, defined as in Equation (3.4). Then, since $\mathbf{B}$ is square and non-singular, we can write from Equation (4.5).(7) $$\begin{array}{cc}\hfill {\mathbf{B}\mathit{\theta }}_{m+1}& ={\mathbf{B}\mathit{\theta }}_m+\mathbf{B}{[\mathbf{I}({\mathit{\theta }}_m)]}^{-1}{[\mathbf{B}}^T{(\mathbf{B}}^T{)}^{-1}{]\mathbf{U}(\mathit{\theta }}_m)\hfill \\ \hfill {\mathbf{B}\mathit{\theta }}_{m+1}& ={\mathbf{B}\mathit{\theta }}_m+\{\mathbf{B}{[\mathbf{I}({\mathit{\theta }}_m)]}^{-1}{\mathbf{B}}^T\}[{({\mathbf{B}}^T)}^{-1}\mathbf{U}({\mathit{\theta }}_m)]\hfill \\ \hfill {\mathit{\theta }}_{m+1}^{\ast }& ={\mathit{\theta }}_m^{\ast }+{[\mathbf{I}({\mathit{\theta }}_m^{\ast })]}^{-1}{\mathbf{U}(\mathit{\theta }}_m^{\ast })\hfill \end{array}$$(7)

where$$\begin{array}{cc}\hfill {\mathit{\theta }}_m^{\ast }& ={\mathbf{B}\mathit{\theta }}_m\hfill \\ \hfill \mathbf{I}({\mathit{\theta }}_m^{\ast })& =\mathbf{B}{[\mathbf{I}({\mathit{\theta }}_m)]}^{-1}{\mathbf{B}}^T\text{is asymptotically diagonal}\hfill \\ \hfill \mathbf{U}({\mathit{\theta }}_m^{\ast })& ={({\mathbf{B}}^T)}^{-1}\mathbf{U}({\mathit{\theta }}_m)\hfill \end{array}$$

We shall call Equation (4.6) the Gram–Schmidt–Fisher scoring algorithm.

4.3.1. Algorithmic process

The proposed Gram–Schmit–Fisher scoring algorithm is a 2-stage iterative process that oscillates between Equations (4.3) and (4.6) until convergence and is described iteratively as follows:

1	=	Start with an initial estimate of the original parameter, $\mathit{\theta }$, denoted by ${\mathit{\theta }}^{(0)}$.
2	=	Estimate $\mathbf{B}$ at ${\mathit{\theta }}^{(0)}$ from Equation 4.3 that is $b_{jk}={\sum }_{r=1}^{j-1}{c}_{rj}{I}_{r(j-k)},j=2,\dots ,p;k=1,\dots ,(j-1)$.
3	=	Determine the initial estimate of the orthogonal parameterization, ${\mathit{\theta }}^{\ast (0)}$ from Equation 4.1, i.e. ${\mathit{\theta }}_1^{\ast (0)}={\mathit{\theta }}_1^{(0)}$ , ${\mathit{\theta }}_j^{\ast (0)}={\mathit{\theta }}_j^{(0)}-{\sum }_{k=1}^{j-1}{b}_{j(j-k)}{\mathit{\theta }}_{j-k}^{(0)},j=2,\dots ,p$.
4	=	Update ${\mathit{\theta }}^{\ast }$ to obtain a new value for ${\mathit{\theta }}^{\ast (1)}$, that is ${\mathit{\theta }}^{\ast (1)}={\mathit{\theta }}^{\ast (0)}+{\mathit{\theta }}^{\ast }$.
5	=	Update $\mathit{\theta }$ to obtain a new value for ${\mathit{\theta }}^{(1)}$, recursively as ${\mathit{\theta }}^{(1)}={\mathbf{B}}^{-1}{\mathit{\theta }}^{\ast (1)}$.
6	=	Repeat steps 2 through 5 using ${\mathit{\theta }}^{\ast (1)}$ to obtain ${\mathit{\theta }}^{\ast (2)}$.
7	=	Stop when $|{\mathit{\theta }}^{\ast (n-1)}-{\mathit{\theta }}^{\ast (n)}|<\mathit{\epsilon }$.

4.4. Example

We consider a sample problem which concerns inference about the difference between two exponential means. Let $Y_1$ and $Y_2$ be the independent exponential random variables with means $\mathit{\varphi }$ and $(\mathit{\varphi }+\mathit{\psi })$, respectively. Then, the joint distribution is given by the function$$f(y_1,y_2|\mathit{\varphi },\mathit{\psi })=\frac{1}{\mathit{\varphi }(\mathit{\varphi }+\mathit{\psi })}exp\left(-\left[\frac{y_1}{\mathit{\varphi }}+\frac{y_2}{\mathit{\varphi }+\mathit{\psi }}\right]\right)$$

The score vector is$$\mathbf{U}(\mathit{\varphi },\mathit{\psi })=\left(\begin{array}{c}\frac{\mathit{\partial }l}{\mathit{\partial }\mathit{\varphi }}\\ \frac{\mathit{\partial }l}{\mathit{\partial }\mathit{\psi }}\end{array}\right)=\left(\begin{array}{c}-\frac{n}{\mathit{\varphi }}-\frac{n}{\mathit{\varphi }+\mathit{\psi }}+\frac{{\sum }_{i=1}^n{y}_{1i}}{{\mathit{\varphi }}^2}+\frac{{\sum }_{i=1}^n{y}_{2i}}{{(\mathit{\varphi }+\mathit{\psi })}^2}\hfill \\ -\frac{n}{\mathit{\varphi }+\mathit{\psi }}+\frac{{\sum }_{i=1}^n{y}_{2i}}{{(\mathit{\varphi }+\mathit{\psi })}^2}\hfill \end{array}\right)$$

The information matrix is$$I_{\mathit{\varphi }\mathit{\psi }}=\left(\begin{array}{cc}{i}_{\mathit{\varphi }\mathit{\varphi }}\hfill & i_{\mathit{\varphi }\mathit{\psi }}\hfill \\ i_{\mathit{\varphi }\mathit{\psi }}\hfill & i_{\mathit{\psi }\mathit{\psi }}\hfill \end{array}\right)=\left(\begin{array}{cc}\frac{n}{{(\mathit{\varphi }+\mathit{\psi })}^2}& \frac{n}{{(\mathit{\varphi }+\mathit{\psi })}^2}\\ \frac{n}{{(\mathit{\varphi }+\mathit{\psi })}^2}& \frac{n}{{\mathit{\varphi }}^2}+\frac{n}{{(\mathit{\varphi }+\mathit{\psi })}^2}\end{array}\right)$$

and its inverse, the variance–covariance matrix, is$$I^{\mathit{\varphi }\mathit{\psi }}=\left(\begin{array}{cc}{i}^{\mathit{\varphi }\mathit{\varphi }}\hfill & i^{\mathit{\psi }\mathit{\varphi }}\hfill \\ i^{\mathit{\psi }\mathit{\varphi }}\hfill & i^{{}_{{}^{\mathit{\psi }\mathit{\psi }}}}\hfill \end{array}\right)=\left(\begin{array}{cc}\frac{1}{n}\left(2{\mathit{\varphi }}^2+2\mathit{\varphi }\mathit{\psi }+{\mathit{\psi }}^2\right)& -\frac{1}{n}{\mathit{\varphi }}^2\\ -\frac{1}{n}{\mathit{\varphi }}^2& \frac{1}{n}{\mathit{\varphi }}^2\end{array}\right)$$

4.4.1. Cox and Reid approach

Orthogonal parametrization, following Cox–Reid method, requires solving the differential equation$$\left\{\frac{1}{(\mathit{\varphi }+\mathit{\psi })}+\frac{1}{{\mathit{\varphi }}^2}\right\}\frac{\mathit{\partial }\mathit{\varphi }}{\mathit{\partial }\mathit{\psi }}=-\frac{1}{\mathit{\varphi }+\mathit{\psi }}$$

This can be solved by separation of variables, leading to $a(\mathit{\lambda })=\mathit{\varphi }(\mathit{\psi }+\mathit{\varphi })/(\mathit{\psi }+2\mathit{\varphi }),$ where $a(\mathit{\lambda })$ is an arbitrary function of $\mathit{\lambda }.$ Cox–Reid suggest setting $a(\mathit{\lambda })=e^{\mathit{\lambda }}$ as a suitable choice. Clearly, this does not produce a unique solution for $\mathit{\varphi }$, regardless of the parametrization of $a(\mathit{\lambda }).$ Thus, different reparametrizations may lead to different modified likelihoods, so that a Cox–Reid estimator may not exist or there may be many of them.

4.4.2. Gram–Schmidt approach

The proposed Gram–Schmidt approach requires seeking a transformation $(\mathit{\psi },\mathit{\varphi })$ to $(\mathit{\psi },\mathit{\lambda }(\mathit{\psi },\mathit{\varphi }))$ such that $\mathit{\lambda }=\mathit{\varphi }-b_{21}\mathit{\psi }$ and solving the equation$$i^{\mathit{\varphi }\mathit{\psi }}=b_{21}{i}^{\mathit{\psi }\mathit{\psi }}$$

From the variance–covariance matrix, we obtain$$b_{21}(\mathit{\varphi },\mathit{\psi })=-\frac{{\mathit{\varphi }}^2}{\left(2{\mathit{\varphi }}^2+2\mathit{\varphi }\mathit{\psi }+{\mathit{\psi }}^2\right)}⟹\mathit{\lambda }=\mathit{\varphi }+\left[\frac{{\mathit{\varphi }}^2}{\left(2{\mathit{\varphi }}^2+2\mathit{\varphi }\mathit{\psi }+{\mathit{\psi }}^2\right)}\right]\mathit{\psi }$$

Just like the Cox–Ried approach, it is impractical to obtain or express $\mathit{\varphi }$ uniquely in terms of $\mathit{\psi }$ and $\mathit{\lambda }$ and subsequently formulate a modified approximate orthogonal joint distribution function. Therefore, one could proceed iteratively, using the proposed Gram–Schmidt–Fisher scoring algorithm.

5. Fitting the disposition model

We will present procedures for estimating the parameters in the correlated model Equation (2.1). The calculations are standard, and so we will only outline the results. We write the likelihood of the joint distribution of the ith group as$$L_i(Y_i;\mathit{\theta })=(1-{\mathit{\alpha }}_i)\prod _{j=1}^{n_i}(1-y_{ij})+{\mathit{\alpha }}_i{L}_{0i}(\mathit{\theta }|y_i)$$

where$$\begin{array}{cc}\hfill L_{0i}(\mathit{\theta }|y_i)& =\prod _{j=1}^{n_i}{\mathit{\delta }}_{ij}^{y_{ij}}{(1-{\mathit{\delta }}_{ij})}^{1-y_{ij}}\hfill \\ \hfill {\mathit{\delta }}_{ij}(\mathit{\theta })& ={\displaystyle \frac{1}{1+exp\{-[M({\mathbf{Z}}_i)+D({\mathbf{Z}}_i)+W({\mathbf{X}}_{ij})]\}}}\hfill \\ \hfill {\mathit{\alpha }}_i(\mathit{\theta })& ={\displaystyle \frac{1+exp[-(M({\mathbf{Z}}_i)+D({\mathbf{Z}}_i)]}{1+exp[-M({\mathbf{Z}}_i)]}}\hfill \end{array}$$

Let$${\mathrm{{\rm Y}}}_{ij}(\mathit{\theta })=M({\mathbf{Z}}_i)+D({\mathbf{Z}}_i)+W({\mathbf{X}}_{ij})$$

and define the following$${\mathrm{{\rm Y}}}_{ij}^{(1)}=\frac{\mathit{\partial }}{\mathit{\partial }\mathit{\theta }}{\mathrm{{\rm Y}}}_{ij}(\mathit{\theta })\text{and}{\mathrm{{\rm Y}}}_{ij}^{(2)}=\frac{{\mathit{\partial }}^2}{\mathit{\partial }{\mathit{\theta }}^T\mathit{\partial }\mathit{\theta }}{\mathrm{{\rm Y}}}_{ij}(\mathit{\theta })$$

The contribution of the i-th group to the log-likelihood is the term$$log{L}_i(\mathit{\theta }|y_i)=log\left[\left\{1-\mathit{\alpha }(\mathit{\theta })\right\}\prod (1-y_{ij})+\mathit{\alpha }(\mathit{\theta })L_{0i}(\mathit{\theta }|y_i)\right]$$

and the contribution to the score function is the term$${\mathbf{U}}_{{}_i}=\frac{\mathit{\partial }}{\mathit{\partial }\mathit{\theta }}log{L}_i(\mathit{\theta })=A_i(\mathit{\theta }){\mathit{\alpha }}_i^{\ast }+D_i(\mathit{\theta }){\mathbf{U}}_{0i}(\mathit{\theta }|\mathbf{y}),$$

where$$\begin{array}{cc}\hfill {\mathit{\alpha }}_i^{\ast }(\mathit{\theta }|y)& =\frac{\mathit{\partial }}{\mathit{\partial }\mathit{\theta }}log\mathit{\alpha }(\mathit{\theta })={\mathit{\delta }}_{i0}(1-{\mathit{\alpha }}_i)\frac{\mathit{\partial }}{\mathit{\partial }\mathit{\theta }}{M}_i(\mathbf{Z})-(1-{\mathit{\delta }}_{i0})\frac{\mathit{\partial }}{\mathit{\partial }\mathit{\theta }}{D}_i(\mathbf{Z})\hfill \\ \hfill {\mathit{\delta }}_{i0}& ={\displaystyle \frac{1}{1+exp\{-[M({\mathbf{Z}}_i)+D({\mathbf{Z}}_i)]\}}}\hfill \\ \hfill A_i(\mathit{\theta })& ={\mathit{\alpha }}_i(\mathit{\theta })\left[L_{0i}(\mathit{\theta }|y)-\prod (1-y_{ij})\right]/L_i,\hfill \\ \hfill D_i(\mathit{\theta })& =\frac{{\mathit{\alpha }}_i{L}_{0_i}}{L_{{}_i}}\hfill \\ \hfill {\mathbf{U}}_{0i}(\mathit{\theta }|y_i)& =\frac{\mathit{\partial }}{\mathit{\partial }\mathit{\theta }}log{L}_{0i}(\mathit{\theta }|y)=\sum _{j=1}^{n_i}(y_{ij}-{\mathit{\delta }}_{ij}){\mathrm{{\rm Y}}}_{ij}^{(1)};\hfill \end{array}$$

Setting $\mathbf{U}(\mathit{\theta })={\sum }_{i=1}^N{\mathbf{U}}_i=\mathbf{0}$, we obtain the score equations. Closed-form solutions are not possible and so the equations are solved by iterative procedures to obtain the maximum likelihood estimates of the parameters.

The contribution of the i-th group to the Hessian matrix is the term$${\mathbf{H}}_i(\mathit{\theta })=\frac{{\mathit{\partial }}^2{l}_i(\mathit{\theta })}{\mathit{\partial }\mathit{\theta }\mathit{\partial }{\mathit{\theta }}^t}=\sum _{j=1}^{n_i}[(y_{ij}-{\mathit{\delta }}_{ij}){\mathrm{{\rm Y}}}_{ij}^{(2)}-{\mathit{\delta }}_{ij}(1-{\mathit{\delta }}_{ij}){\mathrm{{\rm Y}}}_{ij}^{(1)}{\mathrm{{\rm Y}}}_{ij}^{(1)T}].$$

Estimates of the parameter vector can be obtained by the Newton–Raphson iteration routine, which is given by the updating the formula$${\mathit{\theta }}_{s+1}={\mathit{\theta }}_s-{[\mathbf{H}({\mathit{\theta }}_s)]}^{-1}\mathbf{U}({\mathit{\theta }}_s)$$

where ${\mathit{\theta }}_s$ is the estimate of the sth iteration.

The contribution of the i-th group to the Fisher information matrix is the term$${\mathbf{I}}_i(\mathit{\theta })={\mathit{\alpha }}_i{\mathbf{I}}_{0i}+A_i{\mathit{\alpha }}_i^{\ast }{\mathit{\alpha }}^{\ast T}+D_i[{\mathit{\alpha }}^{\ast }{\mathbf{U}}_{0i}^T+{\mathbf{U}}_{0i}{\mathit{\alpha }}^{\ast T}]-D_i(1-{\mathit{\alpha }}_i){\mathbf{U}}_{0i}{\mathbf{U}}_{0i}^T$$

where$${\mathbf{I}}_{0_i}=\sum _{j=1}^{n_i}{\mathit{\delta }}_{ij}(1-{\mathit{\delta }}_{ij}){\mathrm{{\rm Y}}}_{ij}^{(1)}{\mathrm{{\rm Y}}}_{ij}^{(1)T}$$

and where $A_i,$$D_i$ and ${\mathbf{U}}_{0i}$ are evaluated at $\mathbf{y}=(y_1,y_2,\dots ,y_n)=\mathbf{0}$. The estimates can be alternatively obtained by the Fisher scoring method which is given by$${\mathit{\theta }}_{s+1}={\mathit{\theta }}_s+{[\mathbf{I}({\mathit{\theta }}_s)]}^{-1}\mathbf{U}({\mathit{\theta }}_s)$$

The asymptotic variance–covariance matrix of the parameter estimates, $\mathbf{C}(\mathit{\theta })={[\mathbf{I}(\mathit{\theta })]}^{-1}$, is the inverse of the information matrix. Thus, the transformation matrix for use of the proposed Gram–Schmidt–Fisher scoring algorithm can be obtained as described from Equation 4.4. Specifically, the transformation matrix, a lower triangular matrix with ones along the diagonal, is $\mathbf{B}=(b_{jk})$, where $b_{jk}={\sum }_{r=1}^{j-1}{c}_{rj}{I}_{r(j-k)},$$j=2,\dots ,p;$$k=1,\dots ,(j-1)$. and where $(c_{rj})$ and $(I_{r(j-k)})$ are entries of the variance–covariance matrix ,$\mathbf{C}(\mathit{\theta }),$ and the information matrix,$[\mathbf{I}(\mathit{\theta })],$respectively.

The parameter estimates can subsequently be obtained by the Gram–Schmidt–Fisher scoring algorithm given by$${\mathit{\theta }}_{s+1}^{\ast }={\mathit{\theta }}_s^{\ast }+{[\mathbf{I}({\mathit{\theta }}_s^{\ast })]}^{-1}\mathbf{U}({\mathit{\theta }}_s^{\ast })$$

where$${\mathit{\theta }}_s^{\ast }={\mathbf{B}\mathit{\theta }}_s;\mathbf{U}({\mathit{\theta }}_s^{\ast })={({\mathbf{B}}^T{)}^{-1}\mathbf{U}({\mathit{\theta }}_s);\mathbf{I}({\mathit{\theta }}_s^{\ast })=\mathbf{B}[\mathbf{I}({\mathit{\theta }}_s)]}^{-1}{\mathbf{B}}^T$$

5.1. Application to oesophageal cancer in Chinese families

This application involves the study of oesophageal cancer in 2951 nuclear families collected in Yangcheng County, Shanxi Province in China (Kwagyan, kwagyan). In this study, we consider as a group the nuclear family unit. The outcome variable is whether an individual is affected with oesophageal cancer or not. The objective of the study is to assess the presence and aggregation of oesophageal cancer in these families. Table summarizes the distribution of number of affected individuals by family size. Of the 2951 families, 1580 (53%) had no affected individuals. The combined total number of individuals from the studied families was 14310, with mean $\pm$ sd age of 48.26 $\pm$ 18.17 years. Males comprise 56.4 and 25.4% indicated drinking alcohol.

Table 1. Distribution of family size by number of affecteds

	Number of Affecteds (a)
Family size (s)	0	1	2	3	4	5	6	7	Total
3	435	151	34	3					623
4	536	215	52	16	6				819
5	335	203	81	34	8	0			659
6	159	155	59	27	1	4	0		412
7	74	95	46	12	5	3	0	1	232
8	30	54	25	11	2	3	1	0	129
9	6	21	10	3	1	1	0	0	43
10	4	8	6	2	1	1	1	0	23
11	0	3	1	1	0	1	1	0	8
12	0	1	0	1	0	0	0	0	2
13	1	0	0	0	0	0	0	0	1
Total	1580	906	314	110	24	13	3	1	2951

The respondent within a family has correlated outcomes which are influenced in part or wholly by the group variables as well as the variables on the individual respondent. The main objective here is to assess the presence of familial aggregation of oesophageal cancer adjusting for measured risk factors: sex, age and alcohol consumption. Here, the model for the regression analysis of disposition is parametrization as$$M(Z)=\mathit{\lambda },D(Z)=\mathit{\gamma }\text{and}W(X)={\mathit{\beta }}_1\ast \text{sex}+{\mathit{\beta }}_2\ast \text{alcohol}+{\mathit{\beta }}_3\ast \text{age}$$

In this application, the vector of parameters $\mathit{\beta }=({\mathit{\beta }}_1,{\mathit{\beta }}_2,{\mathit{\beta }}_3)$ is the set of parameters of interest. Computations are performed using the computer program we developed CORRDAT (Bonney, Kwagyan, Slater, & Slifker, bonney5), which was linked with the likelihood optimization software MULTIMAX (Bonney, Kwagyan, & Apprey, bonney4). Computations can also be accomplished in MATLAB. Table gives estimates of the correlation matrices. The correlations between the original parameters are quite high compared to those of the orthogonal parameters. In particular, the correlation between $\mathit{\lambda }$ and $\mathit{\gamma }$ is $-0.179$; the correlation between $\mathit{\lambda }$ and ${\mathit{\beta }}_1$ is $-0.414$; and that between $\mathit{\lambda }$ and ${\mathit{\beta }}_3$ is high, $-0.841$. The correlations between the orthogonal parameters are near zero or nonexistent. The correlation between ${\mathit{\lambda }}^{\ast }$ is and ${\mathit{\gamma }}^{\ast }$ is $-0.026$; the correlation between ${\mathit{\lambda }}^{\ast }$ and ${\mathit{\beta }}_1^{\ast }$ is now low, 0.004 and that between ${\mathit{\lambda }}^{\ast }$ and ${\mathit{\beta }}_3^{\ast }$ is also low, 0.0188.

As expected, the orthogonal likelihood converged in fewer iterations than the non-orthogonal one—the orthogonal likelihood converged after 17 iterations; the non-orthogonal one converged after 56 iterations. Estimates of the parameters are summarized in Table . The maximum likelihood estimate of the relative dependence parameter, $\widehat{\mathit{\gamma }}=0.577$ with a corresponding asymptotic 95% confidence interval of (0.529, 0.625) suggests that there is significant familial aggregation of oesophageal cancer in the families sampled. Sex and age have a positive significance, while alcohol was negative. The results suggest that males are at a higher risk of getting oesophageal cancer than females and also that it is more prevalent in older people. The negative effect of alcohol seems to suggest that it has the propensity to lower the risk of oesophageal cancer in the Chinese population studied, perhaps if drank in moderation. In summary, we conclude that oesophageal cancer aggregates in the families sampled.

Table 2. Correlation matrix from original and orthogonal parametrization. Fisher scoring algorithm used for original parametrization and Gram–Schmidt–Fisher scoring algorithm used of orthogonal parametrization

Original parametrization			Orthogonal parametrization
	$\mathit{\lambda }$	$\mathit{\gamma }$		${\mathit{\lambda }}^{\ast }$	${\mathit{\gamma }}^{\ast }$
$\mathit{\lambda }$	1		${\mathit{\lambda }}^{\ast }$	1
$\mathit{\gamma }$	$-0.179$	1	${\mathit{\gamma }}^{\ast }$	0.026	1
${\mathit{\beta }}_1$	$-0.414$	$-0.061$	${\mathit{\beta }}_1^{\ast }$	0.004	0.000
${\mathit{\beta }}_2$	$-0.086$	$-0.036$	${\mathit{\beta }}_2^{\ast }$	$-0.016$	0.001
${\mathit{\beta }}_3$	$-0.841$	0.075	${\mathit{\beta }}_3^{\ast }$	$-0.188$	0.024

Table 3. Regression analysis of disposition to oesophageal cancer. Fisher scoring algorithm used for original parametrization and Gram–Schmidt–Fisher scoring algorithm for orthogonal parametrization

Original parametrization		Orthogonal parametrization
Parameter	Estimates (SE)	Parameter	Estimates (SE)
$\mathit{\lambda }$	$-4.262(0.103)$	${\mathit{\lambda }}^{\ast }$	$-1.873(0.028)$
$\mathit{\gamma }$	0.577 (0.024)	${\mathit{\gamma }}^{\ast }$	$-0.008(0.038)$
${\mathit{\beta }}_1$	0.911 (0.056)	${\mathit{\beta }}_1^{\ast }$	0.915 (0.054)
${\mathit{\beta }}_2\dots$	$-1.131(0.166)$	${\mathit{\beta }}_2^{\ast }$	$-1.128(0.171)$
${\mathit{\beta }}_3$	0.039 (0.002)	${\mathit{\beta }}_3^{\ast }$	0.038 (0.002)
log (Likelihood)	$-5371.294$	log (Likelihood)	$-5377.876$
No. of iterations	56	No. of iterations	27

6. Conclusions

Parameter orthogonalization is used as an aid in computation, approximation, interpretation and elimination or reduction of the effect of nuisance parameters. The resulting algorithm that we have proposed based on the Gram–Schmidt orthogonalization process is computationally feasible. Unlike the approach of Cox–Reid which requires solving a system of differential equations to obtain orthogonality, the proposed method only requires solving a system of linear equations. Our approach has some similarities with the conjugate gradient algorithm, but is distinctly different from it. The conjugate gradient method is an optimization routine, whereas the method we have proposed is principally an approximate orthogonalization technique combined with a maximization algorithm to aid in parameter estimation. Amari (amari1) claims that global orthogonalization is possible only in a special case. Our method allows for global orthogonalization, if the information matrix can be found or estimated. The approach we have proposed is based on exact calculations. The transformation is surjective, that is the original parameters can be sufficiently obtained. Clearly, the process can be cumbersome if a large number of parameters must be accurately estimated because of the inversion of the information matrix at every stage to obtain the covariance matrix and subsequently the transformation matrix. The advantages, however, are that convergence is generally rapid in iteration times, sure and accurate. Block parametrization of parameters would be desirable if interest is in sets of parameters and where the dimension of the parameter space is large. The work we have presented in this article is the first attempt to consider the use of Gram–Schmidt process for estimation of the parameters. It is possible, however, that closer scrutiny, practical considerations, numerical studies for numerical stability and properties, and simulation studies to evaluate and compare convergence times would suggest modifications and refinements to the methods we have discussed.

In conclusion, the algorithm that we have proposed provides an exceedingly convenient method for multi-parameter statistical problems first to reduce numerical difficulties and second improve accuracy in parameter estimation. The method would be efficient for use of function minimization in both linear and non-linear maximum likelihood estimations and particularly useful for small parameter space estimation problems.

Additional information

Funding

The work was supported in part from NIH/NCATS [grant number UL1TR000101 previously UL1RR031975] and NIH/NIMHHD [grant number G12MD007597].

Notes on contributors

John Kwagyan

Dr John Kwagyan is a mathematician and a public health and medical researcher. His research interests include mathematical and statistical modelling of correlated data, clinical trials, survival models, big data analytics, statistical genetics and pharmacokinetic & pharmacodynamics modeling. He is currently the director of Biostatistics, Epidemiology and Research Design, of the Georgetown-Howard University Center for Clinical and Translational Science, Howard University College of Medicine and adjunct professor in the Department of Mathematics, Howard University, Washington, DC, USA.