Кластерные (сгруппированные) стандартные ошибки Максимальное правдоподобие R

Question

Я делаю следующую оценку максимального правдоподобия, используя функцию mle2 из пакета bbmle:

llik.probit2<- function(aL,beta, Ks, Kw, Bs, Bw, dta){
  Y <- as.matrix(dta$qualified)
  sce1 <- as.matrix(dta$eds)
  wce1 <- as.matrix(dta$edw)
  sce1_obs <- (as.matrix(dta$eds_obs))
  wce1_obs <- (as.matrix(dta$edw_obs))
  obs <- as.matrix(dta$CombinedObservable1)
  c <- as.matrix(dta$const)
  phi <- pnorm(ifelse(Y == 0, -1, 1) * (-aL*c + beta*obs + (Bs+Bs*Ks-aL)*sce1 -beta*Ks*sce1_obs + (Bw+Bw*Kw-aL)*wce1 - beta*Kw*wce1_obs), log.p = TRUE)
  -sum(phi)
}

starting.pointmle2 <- list(aL=1.4, beta=0.3, Ks=0.5, Kw=0.5, Bs=0.5, Bw=0.5)

result1 <- mle2(llik.probit2, start = starting.pointmle2, data=list(dta=Mydata), skip.hessian=FALSE)

Мне нужно оценить кластерные стандартные ошибки этой модели, кластеризованные по переменной c_id. Я пытаюсь реализовать сэндвич- оценщик, но не могу получить часть "мяса" из вывода mle2. Я также попробовал с nlm и optim (mle2 в основном является оберткой для этих методов). Любой совет, как это сделать? В общем, есть ли пакет, который предоставляет функции для вычисления кластерных стандартных ошибок для оценки MLE общей функции?