4
$\begingroup$

I'm doing the following maximum likelihood estimation using mle2 function from bbmle package:

library(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)

I need to estimate clustered standard errors of this model, clustering on variable c_id. I'm trying to implement the sandwich estimator but cannot retrieve 'meat' part from mle2 output. I have also tried with nlm, and optim (mle2 is basically wrapper for these methods) Any advice how to do it? In general is there any package which provides functions to calculate clustered standard errors for MLE estimation of a general function?

Here is is a small sample of the data file (MyData) in csv format:

const,qualified,eds,edw,eds_obs,edw_obs,CombinedObservable1,c_id
1,0,0,1,0,0,-0.6838316166,1
1,1,0,1,0,0,-0.1433684328,1
1,0,0,1,0,0.0758113685,0.0758113685,1
1,0,1,0,0.2084778637,0,0.2084778637,34
1,0,0,1,0,0,-0.1622519262,34
1,0,0,1,0,0,-0.5061082675,34
1,0,0,1,0,0,-0.6952748613,34
1,0,1,0,0.9883178986,0,0.9883178986,34
1,0,0,1,0,0,-0.5311805315,34
1,0,0,1,0,0,2.7546881325,34
1,1,1,0,-0.2174263974,0,-0.2174263974,34
1,0,0,1,0,0,-0.4397037288,1
1,0,0,1,0,0,-0.3189328097,1
1,0,0,1,0,0,-0.6132276964,12
1,0,0,1,0,0,0.2459941348,12
1,0,0,1,0,0.0589196966,0.0589196966,12
Improve this question
$\endgroup$
2
  • $\begingroup$ Do you have an idea how can we apply two-way clustering using with Robcov function? $\endgroup$
    – Onurhan Kurt
    Commented Sep 25, 2019 at 7:07
  • $\begingroup$ I faced a similar problem to which Achim Zeiles offered a solution that works for me: stats.stackexchange.com/questions/451853/… $\endgroup$
    – shofla
    Commented Mar 2, 2020 at 17:46

1 Answer 1

Reset to default
2
$\begingroup$

With out sample data it's hard to test what you are looking for given your code, but I would suggest using the rms package. It is a great package and has what you are looking for.

In particular,

lrm : Fit binary and proportional odds ordinal logistic regression models using maximum likelihood estimation or penalized maximum likelihood estimation

robcov: Uses the Huber-White method to adjust the variance-covariance matrix of a fit from maximum likelihood or least squares, to correct for heteroscedasticity and for correlated responses from cluster samples

Here is an example using your data provided with a very simple model:

Sample data

df <- structure(list(const = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L), qualified = c(0L, 1L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L), eds = c(0L, 0L, 
0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L), edw = c(1L, 
1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L), 
    eds_obs = c(0, 0, 0, 0.2084778637, 0, 0, 0, 0.9883178986, 
    0, 0, -0.2174263974, 0, 0, 0, 0, 0), edw_obs = c(0, 0, 0.0758113685, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0589196966), CombinedObservable1 = c(-0.6838316166, 
    -0.1433684328, 0.0758113685, 0.2084778637, -0.1622519262, 
    -0.5061082675, -0.6952748613, 0.9883178986, -0.5311805315, 
    2.7546881325, -0.2174263974, -0.4397037288, -0.3189328097, 
    -0.6132276964, 0.2459941348, 0.0589196966), c_id = c(1L, 
    1L, 1L, 34L, 34L, 34L, 34L, 34L, 34L, 34L, 34L, 1L, 1L, 12L, 
    12L, 12L)), .Names = c("const", "qualified", "eds", "edw", 
"eds_obs", "edw_obs", "CombinedObservable1", "c_id"), class = "data.frame", row.names = c(NA, 
-16L))

Regression w/o clustering:

library(rms)
d <- datadist(df)
options("datadist" = "d")
mod <- lrm(qualified ~ eds, data = df, x = TRUE, y = TRUE)

Logistic Regression Model

lrm(formula = qualified ~ eds, data = df, x = TRUE, y = TRUE)
                     Model Likelihood     Discrimination    Rank Discrim.    
                        Ratio Test            Indexes          Indexes       
Obs            16    LR chi2      1.19    R2       0.135    C       0.679    
 0             14    d.f.            1    g        0.582    Dxy     0.357    
 1              2    Pr(> chi2) 0.2760    gr       1.790    gamma   0.714    
max |deriv| 1e-09                         gp       0.083    tau-a   0.083    
                                          Brier    0.099                     

          Coef    S.E.   Wald Z Pr(>|Z|)
Intercept -2.4849 1.0408 -2.39  0.0170  
eds        1.7918 1.6073  1.11  0.2649

Cluster around c_id:

robcov(mod, df$c_id, method = c('huber'))

Logistic Regression Model

lrm(formula = qualified ~ eds, data = df, x = TRUE, y = TRUE)
                      Model Likelihood     Discrimination    Rank Discrim.    
                         Ratio Test            Indexes          Indexes       
Obs             16    LR chi2      1.19    R2       0.135    C       0.679    
 0              14    d.f.            1    g        0.582    Dxy     0.357    
 1               2    Pr(> chi2) 0.2760    gr       1.790    gamma   0.714    
Cluster on df$c_id                         gp       0.083    tau-a   0.083    
Clusters         3                         Brier    0.099                     
max |deriv|  1e-09                                                            

          Coef    S.E.   Wald Z Pr(>|Z|)
Intercept -2.4849 0.8250 -3.01  0.0026  
eds        1.7918 0.8250  2.17  0.0299
Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks, I will check it out. I have also included a small data sample (csv) $\endgroup$
    – Nidjsi
    Commented Feb 21, 2016 at 21:51
  • $\begingroup$ @Nidjsi I've updated the answer to provide a simple model with the package. Hope this helps. $\endgroup$
    – Amstell
    Commented Feb 21, 2016 at 23:31
  • $\begingroup$ Thanks Amstell. Unfortunately, I think this solution is not suitable for the problem I have. If you look at my specification of log-likelihood function, I'm actually estimating the parameters of the following probit model. $Prob(Y=1|sce1, wce1, sce1\_obs, wce1\_obs, obs)= \Phi(-aL + \beta *obs + (Bs + Bs* Ks -aL)*sce1 -\beta *Ks*sce1\_obs+(Bw+Bw *Kw-aL)*wce1 - \beta *Kw *wce1\_obs)$ and I'm interested in parameters $a_L, \beta, Bs, Ks, Bw, Kw$. $\endgroup$
    – Nidjsi
    Commented Feb 22, 2016 at 16:22
  • $\begingroup$ [Part 2]If I use one among many R's function to estimate probit (i.e. glm) I don't get the estimates of my parameters of interest, but I get the estimates $(\alpha_0, \alpha_1, \alpha_2, \alpha_3, \alpha_4, \alpha_5)$ of corresponding parameters of a different model $Prob(Y=1|sce1, wce1, sce1\_obs, wce1\_obs, obs)= \Phi(\alpha_0 + \alpha_1 *obs + \alpha_2*sce1 -\alpha_3*sce1\_obs+\alpha_4*wce1 - \alpha_5wce1\_obs )$. This is because in the case of solution you have proposed I cannot write the expression formula in function lrm appropriately. $\endgroup$
    – Nidjsi
    Commented Feb 22, 2016 at 16:27
  • $\begingroup$ [Part 3] Given that I'm not interested in $\alpha$'s but in other parameters, I esitmate the model by maximizing log-likelihood function, rather than using glm or lrm function. However I don't know how to use the output from that procedure to solve my problem (for example by using robcov function as you do in your solution) $\endgroup$
    – Nidjsi
    Commented Feb 22, 2016 at 16:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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