I am trying to re-create some indicators from the World Bank, using the methodology described in this paper, and I need to do maximum likelihood estimation, preferably using R.

The aim is to get an overall (unobserved) indicator from many different (observed) scores. My data looks like this:

   Country   Score Source
   <chr>     <dbl> <chr> 
 1 Albania   0.6   EIU   
 2 Albania   0.797 GCS   
 3 Albania   0.833 WMO   
 4 Albania   0.756 PRS   
 5 Algeria   0.45  EIU   
 6 Algeria   0.739 GCS   
 7 Algeria   0.5   WMO   
 8 Algeria   0.634 PRS   
 9 Argentina 0.65  EIU   
10 Argentina 0.725 GCS   
11 Argentina 0.75  WMO   
12 Argentina 0.795 PRS

Based on the formulas and variables given in the paper, I am estimating the following model:

$y_j{_k} = \alpha_k + \beta_k(g_j + \varepsilon_j{_k})$

Where $y{_j}{_k}$ is the observed score of country $j$ from source $k$. This observed score is a linear function of an unobserved indicator, $g_j$, for each country $j$ and a disturbance term, $\varepsilon_j{_k}$, where $\alpha_k$ and $\beta_k$ are parameters which map the unobserved indicator in country $j$, $g_j$ into the observed data $y_j{_k}$. $g_j$ is assumed to be a normally-distributed random variable with mean zero and variance one. We also assume the error term is normally distributed, with zero mean and a variance that is the same across countries, but differs across indicators, i.e. $V[\varepsilon_j{_k}] = \sigma_k^{2}$. We also assume that the errors are independent across sources, i.e., $E[\varepsilon_j{_k}\varepsilon_j{_m}] = 0$ for source $m$ different from source $k$.

I want to estimate the parameters $\alpha_k$, $\beta_k$, and $\sigma_k^2$ for every data source $k$, using maximum likelihood. The paper I referenced says that the contribution to the log-likelihood of country $j$ is given by:

$ln L(\alpha, \beta, \sigma^2) \propto ln|\Omega| + (y_j - \alpha)'\Omega^{-1}(y_j-\alpha)$

Where $\alpha = \Big(\alpha_1, ..., \alpha_K{_j}\Big)$, $\beta = \Big(\beta_1, ..., \beta_K{_j}\Big)$, and $\sigma^2 = \Big(\sigma^2_1, ..., \sigma^2_K{_j}\Big)$, and $B$ and $\Sigma$ are diagonal matrices with $\alpha$ and $\sigma^2$ on the diagonal. Furthermore, the mean of the vector of observed data for each country $j$, $y_j$, is $\alpha$ and the variance is $\Omega = \beta\beta' + B\Sigma B'$. Summing over all countries $j$ and then maximizing over the unknown parameters should give the maximum-likelihood estimates of $\alpha_k$, $\beta_k$, and $\sigma^2_k$ for every data source $k$.

There are a couple of things that confuse me about this notation:

  • are $\alpha$, $\beta$, and $\sigma^2$ matrices or vectors? The paper gives them inside of big parens, which I take to be matrix notation. But they only have one dimension, so are they $1$ X $k$, or $k$ X $k$, with $\alpha ... \alpha_K{_j}$ as the diagonal?
  • why are we taking transposes of diagonal matrices (like $B$)?
  • why are there pipes around the log of omega, i.e., $ln|\Omega|$ instead of parens, i.e., $ln(\Omega)$

So, most of the resources I can find for calculating MLE in R are based on pretty simple examples, where you estimate only a few parameters. In this case, I want to estimate whole matrices of parameters, so I am a little unclear on how to do that. Moreover, we are taking the formula given and summing it across every country.

The best attempt I have so far looks like this (where dat is the dataframe with my data):

nll <- function(alpha1, alpha2, alpha3, alpha4, beta1, beta2, beta3, beta4, sigma1, sigma2, sigma3, sigma4){
  alpha <- diag(c(alpha1, alpha2, alpha3, alpha4))
  beta <- diag(c(beta1, beta2, beta3, beta4))
  sigma <- diag(c(sigma1, sigma2, sigma3, sigma4))

  LL <- 0
  for (country in unique(dat$Country)){
    ix <- dat$Country == country

    y <- dat$Score[ix]

    omega <- beta*t(beta) + alpha*sigma*t(alpha)

    LL <- LL + log(omega) + t(y - alpha)*omega^-1*(y - alpha)

stats4::mle(minuslogl = nll, start = list(alpha1=0.5, alpha2=0.5, alpha3=0.5, alpha4=0.5, 
                                          beta1=1, beta2=1, beta3=1, beta4=1, 
                                          sigma1=1, sigma2=1, sigma3=1, sigma4=1))

However, this throws an error when I try to run the estimation. I also dont think I should be hardcoding a separate alpha, beta, and sigma for each score.


1. "Are $\alpha$, $\beta$, and $\sigma^2$ matrices or vectors? The paper gives them inside of big parens, which I take to be matrix notation. But they only have one dimension, so are they $1 \times K$, or $K \times K$, with $\alpha_1,...,\alpha_{Kj}$ as the diagonal?"

A vector is just a matrix with a single column, so it doesn't really matter what you call $\alpha$, $\beta$, and $\sigma^2$, so long as the are used in conformable ways with other matrices or vectors. Sometimes authors write vectors on a single line in the body of text to save space, but they will usually write them as transposed to make it clear they are vectors (e.g. $X^\prime=(x_1,...,x_k)$, but this isn't always the case and some authors get a little sloppy and forget to write the transpose sign, writing a vector as just $X=(x_1,...,x_k)$ in the body text. Often the reader must use context clues and conformability checks to determine how they are treated. In this case, you should consider $\alpha$, $\beta$, and $\sigma^2$ as vectors (of dimensions $K_j\times 1$) because later in the text, the author states that the $B$ and $\Sigma$ are diagonal matrices with $\alpha$ and $\sigma^2$ on the diagonals (which implies $B$ and $\Sigma$ themselves are $K_j\times K_j$) and adds $\beta\beta^\prime$ to $B{\Sigma}B^\prime$ and the only way that would be conformable is if $\beta$ where a vector.

Note: your question asked if the vectors have a dimension of $K$. You should be careful to note that none of these have a $K$ dimension, but they have a $K_j$ dimension (the rows dimension). This is important because $K$ changes from country to country, so in the $j$th country, each of these vectors would be of length $K_j$. So in one country $\alpha$, $\beta$, and $\sigma^2$ might be of length 10 (i.e. $K_j=10$). In another, they might be of length 15 (i.e. $K_j=15$).

2. "Why are we taking transposes of diagonal matrices (like $B$)?"

Transposes are taken on matrices or vectors in order to make them conformable with other matrices and/or vectors to obtain a desired result. In the referenced paper, the transpose is taken on $B$ in order to square the $a_1, a_2,...,a_{K_{j}}$ terms on the diagonal. To see this, note that, as described in the text, $B$, looks like this:

\begin{eqnarray*} & \left[\begin{array}{cccc} a_{1} & 0 & 0 & 0\\ 0 & a_{2} & 0 & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \ldots & a_{K_{j}} \end{array}\right] \end{eqnarray*}

So, $BB^{\prime}$ will look like this:

\begin{eqnarray*} BB^{\prime} & = & \left[\begin{array}{cccc} a_{1} & 0 & 0 & 0\\ 0 & a_{2} & 0 & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \ldots & a_{K_{j}} \end{array}\right]\left[\begin{array}{cccc} a_{1} & 0 & 0 & 0\\ 0 & a_{2} & 0 & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \ldots & a_{K_{j}} \end{array}\right]=\left[\begin{array}{cccc} a_{1}^{2} & 0 & 0 & 0\\ 0 & a_{2}^{2} & 0 & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \ldots & a_{K_{j}}^{2} \end{array}\right] \end{eqnarray*}

At one point in text you see $B{\Sigma} B^\prime$. This is doing the same thing as above, only adding the $\sigma_j$ terms on the diagonal since $\Sigma$ is diagonal. In another place in the paper, the vector $\beta$ is transposed and premultiplied by $\beta$. This is done to create a $K_j\times K_j$ matrix with the squared $\beta_j$ terms on the diagonal and the $\beta_j \beta_k$ terms on the off diagonals where $i\ne k$ .

3. "why are there pipes around the log of omega, i.e., $ln|\Omega|$ instead of parens, i.e., $ln(\Omega)$"

In mathematical notation, two vertical bars around a matrix like $\Omega$ indicates that the determinant of the matrix argument should be taken. So, in this case $|\Omega|$ means that the determinant of $\Omega$. This is also sometimes written $det(\Omega)$.

Now that you know the answers to 1, 2, and 3 above, I'd suggest, you try to revisit your R code. If you still cannot get your R code to work, I would suggest you accept this answer (if it's helpful) and then post a totally separate question regarding any problems you run into with R. As an aside, I would recommend you consult with a statistician on this. From a very cursory review of the paper you referenced (I didn't read the whole paper, of course), it looks like the authors are building their model using linear mixed effects models. I think this is the case since I saw terms like $\Omega=\beta \beta^\prime + B\Sigma B$ on page 24 and this type of equation usually pops up in linear mixed effects models. This is important, because if this is the case, you don't need to re-invent the wheel with your R code and computer anything manually. Instead, you could probably invoke the appropriate mixed model commands (e.g. those found in the pages like lme4, nlme) after simply reading in your data, and the program will compute everything for you - so probably no need to write that big ugly function you seem to be working on.

  • $\begingroup$ Thanks for the feedback @StatsStudent. Your response to #3 and your suggestion to look into mixed models definitely help. I'm still confused about the transposes though. For part 1, if $\beta$ is a $1$X$K_j$ matrix, how can I take $\beta\beta'$? Wouldnt they be non conformable (ie, one is $1$X$K_j$ and one is $K_j$X$1$)? Also, for part 2, how is $BB'$ different from just $BB$ (without transposing the second $B$)? Aren't they identical? $\endgroup$ – Amadou Kone Jan 7 at 21:41
  • 1
    $\begingroup$ No problem. No, if you are considering $\beta$ a matrix, then it is $K_j\times 1$ not vice versa. I think you are better off calling it "vector" (although technically either is right, but usually what we call $N\times 1$ structures a vector). So then $\beta \beta^\prime$ is conformable. When you post-multiply a $K_j\times 1$ vector by a $1 \times K_j$ (transposed) vector the result is a $K_j\times K_j$ matrix. In part two, $BB=BB^\prime$ in this particular case since $B$ is diagonal. This is not generally always the case though. $\endgroup$ – StatsStudent Jan 7 at 21:47

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.