Maximum Likelihood Estimator for equicorrelation model

Question

Consider the equicorrelation model for multivariate normal. Let $X_1, \dots, X_n\sim \mathbf{N_p(\mu, \Sigma)}$, where $\mathbf{\Sigma}=\sigma^2((1-\rho)\mathbf{I_p}+\rho\mathbf{J_p})$ where $\mathbf{J_p}$ is a matrix of all $1$'s. Now how do we find the maximum likelihood estimator of $\sigma^2, \rho$?. Estimating $\mathbf\mu$ is easier, it comes out to be $\overline{X}$. Can someone help me out? If I take inverses and try to differentiate the likelihood it gets quite messy. Is there a slick method? Thanks a lot.

Answer 1

The useful formulas have been given in the answer by @Alecos Papadopoulos, which I will refer to. Note that $\mathbf{\Sigma}$ is positive definite iff $\rho_{\text{min}} < \rho < 1$ with $\rho_{\text{min}} := - 1/(p-1)$.

For the ML estimation you can concentrate the vector $\boldsymbol{\mu}$ and the scalar $\sigma^2$ out of the log-likelihood function by replacing them by their ML estimates, namely $\widehat{\boldsymbol{\mu}} = \bar{\mathbf{X}}$ and $$ \widehat{\sigma}^2 = \frac{1}{np} \, \sum_{i = 1}^n \left[\mathbf{X}_i - \bar{\mathbf{X}} \right]^\top \boldsymbol{\Sigma}^{\star -1} \left[\mathbf{X}_i - \bar{\mathbf{X}}\right] $$ where $\boldsymbol{\Sigma}^\star:= \sigma^{-2} \boldsymbol{\Sigma}$. Thus the estimation boils down to the one-dimensional optimisation of a function of $\rho$.

An interesting point is that for a vector $\mathbf{z}$ of length $p$ we have $$ \mathbf{z}^\top \boldsymbol{\Sigma}^{\star -1} \mathbf{z} = \frac{1}{1-\rho} \left[\mathbf{z}^\top \mathbf{z} - \nu (\mathbf{1}_p^\top \mathbf{z})^2 \right] $$ where $\mathbf{1}_p$ is a vector of $p$ ones and $\nu:= \rho\, / \, [1 + (p-1) \rho]$. So taking $\mathbf{z}_i:= \mathbf{X}_i - \bar{\mathbf{X}}$ we can use the sums $$ A := \sum_i \mathbf{z}_i^\top \mathbf{z}_i, \qquad B := \sum_i (\mathbf{1}_p^\top \mathbf{z}_i)^2 $$ which do not depend on the parameter. The concentrated log-likelihood obtained after replacing $\boldsymbol{\mu}$ and $\sigma^2$ by their estimate is $$ \ell_{\text{c}}(\rho) = -\frac{np}{2} \, \log(2 \pi) - \frac{n}{2}\,\log \widehat{\sigma}^2 - \frac{n}{2}\,\log|\boldsymbol{\Sigma}^\star| - \frac{np}{2}. $$ which is easily maximised.

##' MLE of a normal vector with compound symmetric covariance.
##' @title Compound symmetry
##' @param X A matrix with n row and p columns.
##' @return A list of ML Estimates.
MLE <- function(X) {

    n <- nrow(X)
    p <- ncol(X)
    p1  <- p - 1
    rho.min <- - 1 / p1

    ## substract column means. We could use 'scale' as well
    xbar <- apply(X, 2, mean)
    Xcent <- sweep(X, MARGIN = 2, STATS = xbar, FUN = "-")
    
    ## sum of 'n' crossprods and sum of 'n' squared sum of 'p'
    ## components
    A <- sum(apply(Xcent, MARGIN = 1, FUN = crossprod))
    B <- sum(apply(Xcent, MARGIN = 1, sum)^2)
    
    ## concentrated (or profiled) deviance = - 2 log L, to be
    ## minimised
    deviance  <- function(rho1) {
        nu <- rho1  / (1.0 + p1 * rho1)
        sigma2.hat <- (A - nu * B) / n / p / (1.0 - rho1)
        dev <- n * (p * log(sigma2.hat) + p1 * log(1.0 - rho1)
                    + log(1.0 + p1 * rho1) + p)
        attr(dev, "sigma2.hat") <- sigma2.hat
        dev
    }

    opt <- optimise(deviance,  interval = c(rho.min, 1.0))
    list(mu = xbar,
         sigma2 = attr(opt$objective, "sigma2.hat"),
         rho = opt$minimum,
         loglik = -opt$objective[1] / 2)
}

## Now try the function change/remove set.seed for other examples
set.seed(123)
p <- 10L; n <- 1000
rho.min <- -1 / (p - 1)

## draw 'rho', 'mu' and 'sigma2'
rho <- runif(1, min = rho.min, max = 1)
mu <- rnorm(p)
sigma2 <- rexp(1)

## build a p * p covariance matrix 'Sigma'
Sigma <- matrix(rho, nrow = p, ncol = p)
diag(Sigma) <- 1
Sigma <- sigma2 * Sigma

## build the matrix 'X'
G <- chol(Sigma)
X <- matrix(rnorm(n * p), nrow = n, ncol = p) %*% G
X <- sweep(X, MARGIN = 2, STATS = mu, FUN = "+")

## 
MLE(X)

Answer 2

One can verify that the inverse of this matrix equals

$$\mathbf{\Sigma}^{-1} =\frac {1}{\sigma^2(1-\rho)}\left(\mathbf{I_p}-\frac{\rho}{1+(p-1)\rho}\mathbf{J_p}\right)$$

Its determinant is

$$|\mathbf{\Sigma}| = \sigma^{2p}\cdot [1+(p-1)\rho]\cdot (1-\rho)^{p-1}$$

(see Tong 1990, p.104)

I don't think it gets simpler than that (and it is not very messy either).

Beware of asymptotics (which are always prominent in maximum likelihood estimation), since the matrix may become singular at the limit. And MLE with a sigular covariance matrix is a whole different ball-game. If I remember correctly, the existing theoretical results are acrimonious to follow in such a case.

Answer 3

I found this question in an exercise of Seber's Multivariate Observations (1984).

I don't have a proof but apparently the maximum likelihood estimators of $\sigma^2$ and $\rho$ are given by

$$\hat\sigma^2=\frac1p\sum_{j=1}^p s_{jj}$$

and $$\hat\rho=\frac1{\hat\sigma^2}\sum_{j\ne k}\frac{s_{jk}}{p(p-1)}\,,$$

where $s_{jk}$ is the $(j,k)$th element of $\boldsymbol S=\frac1{n-1}\sum\limits_{i=1}^n (\boldsymbol X_i-\boldsymbol{\overline X})(\boldsymbol X_i-\boldsymbol{\overline X})^T$.

The relevant extract from page 95 of the book:

The hint given in the exercise is to find $\boldsymbol \Sigma^{-1}$ (which exists for $\rho\in \left(-\frac1{p-1},1\right)$ as pointed out by @Yves) and $\det\boldsymbol \Sigma$, as given in the answer by @Alecos Papadopoulos.

Answer 4

There exists explicit (closed-form) solutions for the MLE's you are struggling to find. In fact, a little intuition, after the re-parametrization in step 1 shown below, reveals that the average of the diagonal elements of xx' is a reasonable estimator for the variance and the average of the off-diagonal elements of xx' is a reasonable estimator for the covariance. But to be specific, follow these steps. For simplicity, let the mean vector be equal to the zero vector.

1) Parameterize the problem so that the diagonal elements of the variance-covariance matrix are each "a" and the off-diagonal elements are each "b". 2) Use well known formulas for the inverse and determinant of the matrix constructed in step 1 to set-up the likelihood equation. 3) Re-parameterize the problem again by letting c = a - b and d = a + b*(p-1). 4) Find the MLE of c and d (very easy): you should obtain c = (x'x)/(n-1) - (x'j)^2/(n(n-1) and d = (x'j)^2/n 5) Use the invariance property of the MLE to show that a = ((n-1)*c + d)/n and b = (d-c)/p