Optimal estimation of a mean from non-independent data

Question

I have the following model:

$Y_1=\beta+\varepsilon_1+\varepsilon_2$

$Y_2=\beta+\varepsilon_3+\varepsilon_4$

$Y_3=\beta+\varepsilon_1+\varepsilon_4+\varepsilon_5$

$Y_4=\beta+\varepsilon_2+\varepsilon_3+\varepsilon_5$

$\varepsilon_i\thicksim \text{iid } \mathcal{N}(0,\sigma^2), \forall i$

I would like to obtain the best (unbiased and with minimum variance) estimator of $\beta$. That is, I would like to know $\hat{\beta}=f(Y_1,Y_2,Y_3,Y_4)$. How should I obtain it?

I will really appreciate your help.

Answer 1

As you've described the problem, ${\boldsymbol Y} = \{Y_1, Y_2, Y_3, Y_4\}$ will have a multivariate normal distribution with mean ${\boldsymbol \mu} = (\beta, \beta, \beta, \beta)'$ and covariance matrix

$$ \Sigma = \sigma^2 \left( \begin{array}{cccc} 2 & 0 & 1 & 1 \\ 0 & 2 & 1 & 1 \\ 1 & 1 & 3 & 1 \\ 1 & 1 & 1 & 3 \\ \end{array} \right) $$

Normally this type of covariance structure model would require some kind of software like MPLUS but I believe it may be simple enough to "trick" lme into fitting a model like this but it is simple enough to "build-your-own".

I'm not sure about getting the unbiased minimum variance estimator (although I'm sure the ordinary sample mean would be competitive), but I can describe how to get the maximum likelihood estimator (MLE), which is desirable for the reasons mentioned by Michael Chernick. The log-likelihood for a single observation of ${\boldsymbol Y}$ is

$$L(\beta, \sigma^2) = \log \left( \frac{1}{(2\pi)^2 |\Sigma|^{1/2}} \right) -\frac{1}{2} ({\boldsymbol Y}-{\boldsymbol \mu})' \Sigma^{-1} ({\boldsymbol Y}-{\boldsymbol \mu}) $$

which is only a function of $\beta$ and $\sigma^2$ since ${\boldsymbol \mu}$ only depends on $\beta$ and $\Sigma$ only depends on $\sigma^2$. We sum over the observations and optimize the resulting function as a function of $\beta, \sigma^2$ to get the MLE. I'll use the dmnorm() function from the R package mnormt to do this and give a rather crudely programmed example:

set.seed(1234) 
N <- 100 
s = matrix(0,4,4)
s[1,]=c(2,0,1,1)
s[2,]=c(0,2,1,1)
s[3,]=c(1,1,3,1)
s[4,]=c(1,1,1,3)

# generate data where true values are beta=1, sigma^2 = 3. 
y <- list()
for(i in 1:N) y[[i]] <- rmnorm(1,mean=c(1,1,1,1),varcov=3*s)

# P[1] is beta, P[2] is sigma squared
L <- function(P)
{
   # crude barrier to prevent sigma squared being negative
   if( P[2] <= 0 ) return(Inf) 

   like <- 0 
   for(i in 1:N) 
   {
      like <- like + dmnorm(y[[i]], mean=rep(P[1],4), varcov=P[2]*s, log=TRUE)
   }
   return(-like)
}
# chose arbitrary starting values of beta=1,sigma^2=1 for the optimization
optim(c(1,1),L)$par
[1] 0.9109401 3.0786393

You can get approximate confidence intervals either by bootstrapping or using the fisher information, which will require derivatives of the log-likelihood either numerically (which is returned by optim()) or analytically, which you may find this thread helpful for.

Answer 2

A general result on gaussian families which should be in your lecture notes says that $\hat\beta$ is an unbiased affine transform of the vector $Y=(Y_k)_{1\leqslant k\leqslant 4}$. Unbiasedness for every value of $\beta$ imposes that this transform must be linear. Since the coefficient of $\beta$ in each $Y_k$ is $1$, one sees that $\hat\beta=\langle x, Y\rangle=\sum\limits_{k=1}^4x_kY_k$ for some $x=(x_k)_{1\leqslant k\leqslant 4}$ such that $\langle x, 1\rangle=\sum\limits_{k=1}^4x_k=1$.

At this point, Lagrange multiplier's method readily yields the value of $x=(x_k)_{1\leqslant k\leqslant 4}$, hence, of $\hat\beta$, but, in the present case, symmetry considerations offer a nice alternative proof.

To see this, note that the symmetry $\varepsilon_1\leftrightarrow\varepsilon_3$, $\varepsilon_2\leftrightarrow\varepsilon_4$, exchanges $Y_1$ and $Y_2$ and exchanges $Y_3$ and $Y_4$. Since the distribution of $Y$ is invariant by this operation, this yields $x_1=x_2$ and $x_3=x_4$. Hence $x_1=x_2=\frac12(1-t)$ and $x_3=x_4=\frac12t$ for some $t$.

For every $x=(x_k)_{1\leqslant k\leqslant 4}$, the variance of $\langle x, Y\rangle=\sum\limits_{k=1}^4x_kY_k$ is $$\sigma^2\cdot((x_1+x_3)^2+(x_1+x_4)^2+(x_2+x_4)^2+(x_2+x_3)^2+(x_3+x_4)^2), $$ and, when $x=(x_k)_{1\leqslant k\leqslant 4}$ is as above, the sum in the parenthesis is $\frac14+\frac14+\frac14+\frac14+t^2$, which is minimum for $t=0$. Finally, all this proves that $$ \hat\beta=\tfrac12(Y_1+Y_2). $$

Answer 3

The $Y_i$s are dependent because some share the same noise terms. However you still have a well defined likelihood function for beta. The theory of maximum likelihood still applies. Determine the likelihood and obtain the mle for beta.

A weighted average of the Yis would be unbiased and you could find the weights that minimize the bias amoung the linear estimators. The mle could be biased but I am reasonably sure it will be consistent, asymptotically normal and achieve the Cramer-Rao lower bound. This should be checked.

For the best linear unbiased estimator of beta let that be $B= a_1Y_1+a_2Y_2+a_3Y_3+ a_4Y_4$.

Then $${\rm Var}(B) =a_1^2 {\rm Var}(Y_1)+a_2^2{\rm Var}(Y_2) + a_3^2 {\rm Var}(Y_3) +a_4^2 {\rm Var}(Y_4) + 2a_1 a_2 {\rm Cov}(Y_1,Y_2)+2a_1 a_3 {\rm Cov}(Y_1, Y_3)+2a_1a_4 {\rm Cov}(Y_1,Y_4)+2a_2a_3 {\rm Cov}(Y_2,Y_3)+2a_2a_4 {\rm Cov}(Y_2,Y_4) +2a_3a_4 {\rm Cov}(Y_3, Y_4)$$.

Now by the definition of the $Y_i$s, ${\rm Cov}(Y_1,Y_2)=0$. All the other covariance terms = sigma square because the Yi pairs each contain one common epsilon.

So $${\rm Var}(B)= \sigma^2(2a_1^2 +2a_2^2+3 a_3^2+3a_4^2 + 2a_1 a_3 +2a_1 a_4 +2 a_2 a_3+2a_2 a_4 +2 a_3 a_4)$$

Now for $B$ to be unbiased we have the constraint $a_1+a_2+a_3+a_4=1$ So the solution does require the use of Lagrange multipliers.