How to cluster standard errors when fitting a model using maximum likelihood estimation in R?

Question

I am estimating several parameters using the mle() function with the L-BFGS-B method in R to data from an experiment in which participants had to make multiple choices. My goal is to obtain population estimates, i.e., pool all the data and find the set of parameters that describes aggregate choices best.

Since each participant provides some of the data, I want to control for heterogeneity in choice error between participants, i.e. control for the notion that person A's choices may be are more erratic than those of person B. As far as I know, I should pursue to cluster the standard errors (SEs) by participant. Whereas I have found packages on how to cluster SEs when using maximum likelihood estimation (MLE) to obtain the coefficients of (non-) linear models, e.g.

• Clustered (grouped) standard errors MLE in R

• https://cran.r-project.org/web/packages/fixest/vignettes/fixest_walkthrough.html

I have not found anything that is applicable to (the output of) mle(). My question has been asked before on the internet, https://stat.ethz.ch/pipermail/r-help/2014-July/376336.html but has remained unanswered so far.

Below is the basic model that I estimate and a small portion of the data. You will notice that the array "subject.id" is not utilised at the moment. My question is how I can incorporate it for clustered SEs.

### Functions
LL <- function(n,a,b,s)
{
  V = (v(z1,n)-v(z2,n))*w(p,a,b) + v(z2,n) 
  res = zce - v.inv(V,n)
  ll = dnorm(res, 0, s,log=T)
  return(-sum(ll))
}

# with:
u <- function(x,n) 
{
  ifelse(n!=1,util <- x^(1-n)/(1-n), util <- log(x))
  return(util)
}
u.inv <- function(x,n)
{
  ifelse(n !=1, inv.util <- ((1-n)*(x))^(1/(1-n)), inv.util <- exp(x))
  return(inv.util)
}

v = function(x,n){return(1/(u(maxz,n)-u(minz,n))*(u(x,n)-u(minz,n)))}
v.inv = function(x,n){return(u.inv(x*(u(maxz,n)-u(minz,n))+u(minz,n),n))}

w <- function(p,a,b){return(exp(-b*(-log(p))^(1-a)))}

maxz = 135
minz = 0

### Data 
z1 <- c(0.1111111, 0.1037037, 0.1222222, 0.1111111, 0.1074074, 0.1666667, 0.1333333, 0.2000000, 0.1333333, 0.1074074,
        0.1037037, 0.1111111, 0.1333333, 0.2000000, 0.1222222, 0.1111111, 0.1666667, 0.1333333, 0.1111111, 0.1333333,
        0.1111111, 0.1666667, 0.1074074, 0.1333333, 0.1222222, 0.2000000, 0.1037037)

z2 <- c(0.08888889, 0.06666667, 0.07777778, 0.00000000, 0.03333333, 0.09259259, 0.09629630, 0.08888889, 0.06666667,
        0.03333333, 0.06666667, 0.08888889, 0.06666667, 0.08888889, 0.07777778, 0.00000000, 0.09259259, 0.09629630,
        0.00000000, 0.09629630, 0.08888889, 0.09259259, 0.03333333, 0.06666667, 0.07777778, 0.08888889, 0.06666667)

p <-  c(0.5, 0.9, 0.5, 0.9, 0.9, 0.1, 0.1, 0.1, 0.5, 0.9, 0.9, 0.5, 0.5, 0.1, 0.5, 0.9, 0.1, 0.1, 0.9, 0.1, 0.5, 0.1, 0.9, 0.5, 0.5, 0.1, 0.9)

zce <- c(0.11055556, 0.10277778, 0.11000000, 0.10833333, 0.10185185, 0.11666667, 0.13240741, 0.14166667, 0.13166667,
         0.07222222, 0.08796296, 0.09944444, 0.09500000,0.10833333, 0.09444444, 0.05277778, 0.10925926, 0.11759259,
         0.05833333, 0.10277778, 0.09277778, 0.10925926, 0.06111111, 0.08833333, 0.09222222, 0.12500000, 0.09166667)

subject.id <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 4, 4, 4, 4, 4, 4, 4, 4)

### mle()
fit <- mle(LL,
           start = list(n = 0.1,a=0.1,b=0.1,s=0.1),
           method = "L-BFGS-B",
           lower = list(n=-Inf,a = -Inf, b = 0.0001, s=0.0001),
           upper = list(n=0.9999,a = 0.9999, b = Inf, s=Inf),
           control = list(maxit = 500, ndeps = rep(0.000001,4)),
           nobs=length(z1))

Some papers within the field I am working in have employed a seemingly simple method to account for individual heterogeneity in errors, but colleagues of mine have pointed out potential issues with this method (which induced me to post this question here). Next, I will explain this method and the critique I have heard regarding it and would be grateful for any further comments.

In order to account for individual heterogeneity in errors, previous literature has estimated an error variance s for each participant individually, like this:

### New LL
LL.id <- function(n,a,b,s1,s2,s4)
{
  V = (v(z1,n)-v(z2,n))*w(p,a,b) + v(z2,n) 
  res = zce - v.inv(V,n)

  sigma <- as.numeric(mget(paste("s",subject.id,sep="")))

  ll = dnorm(res, 0, sigma,log=T)
  return(-sum(ll))
}


### New mle()
fit.id <- mle(LL.id,
           start = list(n = 0.1,a=0.1,b=0.1,s1=0.1,s2=0.1,s4=0.1),
           method = "L-BFGS-B",
           lower = list(n=-Inf,a = -Inf, b = 0.0001, s1=0.0001,s2=0.0001,s4=0.0001),
           upper = list(n=0.9999,a = 0.9999, b = Inf, s1=Inf,s2=Inf,s4=Inf),
           control = list(maxit = 500, ndeps = rep(0.000001,6)),
           nobs=length(z1))

Potential issues with this method:

• Using this method means there is a "huge" amount of free parameters that must be estimated, resulting in a model that is overfitted. Moreover, it becomes less likely that the global optimum is found by the optimisation algorithm, as there are simply too many dimensions.

• It is undesirable that for each set of additional observations (i.e. adding the responses of an additional participant) there is an additional free parameter that must be estimated (this is, apparently, somehow related to the incidental parameters problem?).

To summarise, my questions boil down to this: is it possible to cluster SEs by participant for the mle that I am running?, if yes, how?, and will it be a better way to estimate the parameters compared to the alternative method that I presented?

I welcome any type of advice!

Update 1

For my actual estimation I have a few more parameters and a lot more observations (up to ±5000). As Achim rightly points out, it takes the estfun.mle some time with more data to complete calculation, which is why I have added parallel computing using mclapply from the "parallel" package, which works for mac and linux. For windows I recommend https://www.r-bloggers.com/implementing-mclapply-on-windows-a-primer-on-embarrassingly-parallel-computation-on-multicore-systems-with-r/:

library(parallel)
library(numDeriv)
library(stats4)

estfun.mle <- function(x, ...) {
  form <- formals(x@minuslogl)
  if(names(form)[length(form)] != "i" | !is.null(form[[length(form)]])) {
    stop("cannot compute gradient contributions, last argument of minuslogl() must be i = NULL")
  }
  func <- function(par, i = NULL) -do.call(x@minuslogl, c(as.list(par), list(i = i)))
  func2 <- function(i){grad(func, coef(x), i = i)}
  save1 <- mclapply(1:nobs(x), func2, mc.cores = numcores)
  save2 <- do.call(rbind,save1)
  return(save2)
}

Answer 1

Clustered sandwich standard errors can be computed with the vcovCL() function from the sandwich package, see vignette("sandwich-CL", package = "sandwich"). The ingredients that are needed are:

• A bread() method that extracts the "bread" matrix in the sandwich computation, typically n times the inverse of the Hessian (observed or expected).
• An estfun() method that extracts the "estimating functions" or "gradient contributions", i.e., the score function evaluated at the maximum likelihood estimates of the parameters separately for each observation in the sample.

The former can be easily computed via nobs(fit) * vcov(fit) as nobs() and vcov() methods are provided for "mle" objects. However, the objects do not provide scores or gradients. So either this information must be provided separately, ideally using analytical scores/gradients. A fallback solution would be to use grad() from the numDeriv package to approximate the gradients numerically. This typically works well but is relatively slow (compared to vectorized analytical solutions).

However, there is one problem: The LL() function you pass to mle() only evaluates the sum of the log-likelihood contributions for all observations. But we need the possibility to obtain the results observation by observation. So one trick would be the following:

LL <- function(n, a, b, s, i = NULL)
{
  V <- (v(z1, n) - v(z2, n)) * w(p, a, b) + v(z2, n) 
  res <- zce - v.inv(V,n)
  ll <- dnorm(res, 0, s, log = TRUE)
  if(is.null(i)) i <- 1:length(z1)
  -sum(ll[i])
}

This is essentially the same log-likelihood you used but I added the argument i = NULL. This enables to get the likelihood only for observation i = 1 or i = 2 etc. Then mle() can be called exactly as before:

fit <- mle(LL,
           start = list(n = 0.1, a = 0.1, b = 0.1, s = 0.1),
           method = "L-BFGS-B",
           lower = list(n = -Inf, a = -Inf, b = 0.0001, s = 0.0001),
           upper = list(n = 0.9999,a = 0.9999, b = Inf, s = Inf),
           control = list(maxit = 500, ndeps = rep(0.000001, 4)),
           nobs = length(z1))

The bread() can be extracted by:

bread.mle <- function(x, ...) {
  stats4::nobs(x) * stats4::vcov(x)
}

And the estfun() by:

estfun.mle <- function(x, ...) {
  form <- formals(x@minuslogl)
  if(names(form)[length(form)] != "i" | !is.null(form[[length(form)]])) {
    stop("cannot compute gradient contributions, last argument of minuslogl() must be i = NULL")
  }
  func <- function(par, i = NULL) -do.call(x@minuslogl, c(as.list(par), list(i = i)))
  t(sapply(1:stats4::nobs(x), function(i) numDeriv::grad(func, stats4::coef(x), i = i)))
}

This code is somewhat ugly and condensed but it does the following:

• First, it is checked whether the log-likelihood underlying the "mle" fit has i = NULL as the last argument. If yes, it is assumed that this can be used to select the observation. If no, the function stops.
• Second, we set up the positive log-likelihood func(). This is necessary because the numDeriv package wants a single argument with a vector of parameters whereas mle() wants a list of scalar parameters. This is interfaced suitable in func(par, i = NULL).
• Third, numDeriv::grad() is called for this func() at the estimated parameters and for i = 1, ..., i = nobs(fit).

But with this we can now get the clustered covariance matrix estimate:

library("lmtest")
vcovCL(fit, subject.id)
##               n             a            b             s
## n  0.0848076882 -0.0453881312 -0.064477117  5.871361e-04
## a -0.0453881312  0.0327678769  0.054799305 -5.799149e-04
## b -0.0644771167  0.0547993052  0.097596087 -1.082398e-03
## s  0.0005871361 -0.0005799149 -0.001082398  1.239231e-05

We can also see that except for parameter n this increases the standard errors, e.g., by using the coeftest() function from lmtest for computing Wald tests for all coefficients.

library("lmtest")
coeftest(fit)
## z test of coefficients:
## 
##    Estimate Std. Error z value  Pr(>|z|)    
## n 0.1653341  0.4859535  0.3402    0.7337    
## a 0.6525431  0.1465998  4.4512 8.540e-06 ***
## b 0.7872708  0.1728246  4.5553 5.231e-06 ***
## s 0.0147600  0.0020084  7.3490 1.997e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(fit, vcov = vcovCL, cluster = subject.id)
## z test of coefficients:
## 
##    Estimate Std. Error z value  Pr(>|z|)    
## n 0.1653341  0.2912176  0.5677 0.5702156    
## a 0.6525431  0.1810190  3.6048 0.0003124 ***
## b 0.7872708  0.3124037  2.5200 0.0117340 *  
## s 0.0147600  0.0035203  4.1928 2.755e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wald confidence intervals can be computed with coefci() instead of coeftest().

(P.S.: Rather than computing the parameter estimates under an independence working model and then adjusting the standard errors for clustering, you can also try to model the clustering explicitly in the likelihood. I don't understand enough about your model to understand whether your LL.id() does this correctly or not.)