Multilevel probit models

Stata allows you to fit multilevel mixed-effects probit models with meprobit. A multilevel mixed-effects probit model is an example of a multilevel mixed-effects generalized linear model (GLM). You can fit the latter in Stata using meglm.

Let's fit a crossed-effects probit model. A crossed-effects model is a multilevel model in which the levels of random effects are not nested. We investigate the extent to which two salamander populations, whiteside and roughbutt, cross-breed. We label whiteside males wsm, whiteside females wsf, roughbutt males rbm, and roughbutt females rbf. Our dependent variable y is coded 1 if there was a successful mating and 0 otherwise. Let's fit our model:

. webuse salamander

. meprobit y wsm##wsf || _all: R.male || female:

note: crossed random-effects model specified; option intmethod(laplace) 
      implied.

Fitting fixed-effects model:

Iteration 0:   log likelihood = -223.01026  
Iteration 1:   log likelihood = -222.78736  
Iteration 2:   log likelihood = -222.78735  

Refining starting values:

Grid node 0:   log likelihood = -216.49485

Fitting full model:

Iteration 0:   log likelihood = -216.49485  (not concave)
Iteration 1:   log likelihood = -214.34477  (not concave)
Iteration 2:   log likelihood = -212.30675  
Iteration 3:   log likelihood = -209.21068  (not concave)
Iteration 4:   log likelihood = -209.18271  
Iteration 5:   log likelihood = -208.21198  
Iteration 6:   log likelihood = -208.11598  (not concave)
Iteration 7:   log likelihood = -208.11296  
Iteration 8:   log likelihood = -208.11185  
Iteration 9:   log likelihood = -208.11182  
Iteration 10:  log likelihood = -208.11182

Mixed-effects probit regression                 Number of obs     =        360

        Grouping information




 
     No. of       Observations per group

 Group variable
 
     groups    Minimum    Average    Maximum

 
 
           _all
 
          1        360      360.0        360

         female
 
         60          6        6.0          6




Integration method: laplace

                                                Wald chi2(3)      =      42.97
Log likelihood = -208.11182                     Prob > chi2       =     0.0000



           y 
 
 Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
 
 
 
       1.wsm 
 
  -.4121819   .2580721    -1.60   0.110     -.917994    .0936301
       1.wsf 
 
  -1.720302   .3131195    -5.49   0.000    -2.334005   -1.106599
              
 
     wsm#wsf 
 
        1 1  
 
   2.121092   .3492928     6.07   0.000      1.43649    2.805693
              
 
       _cons 
 
   .5950767   .2214112     2.69   0.007     .1611186    1.029035
 
 
 
_all>male    
 
   var(_cons)
 
   .3867363   .1795286                      .1556951    .9606274
 
 
 
female       
 
   var(_cons)
 
   .4464042   .1980406                      .1871136    1.065004


LR test vs. probit model: chi2(2) = 29.35                 Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.

Our model has two random-effects equations, separated by ||. We use the _all notation that identifies all the observations as one big group. We use the R. notation to tell Stata to treat male as an indicator variable.

The output table includes the fixed-effect portion of our model and the estimated variance components. The estimates of the random intercepts suggest that the heterogeneity among the female salamanders is larger than the heterogeneity among the male salamanders.

If we wish, we can constrain the two random intercepts to be equal.

. constraint 1 _b[/var(_cons[_all>male])] = _b[/var(_cons[female])]

. meprobit y wsm##wsf || _all: R.male || female:, constraint(1) nolog

note: crossed random-effects model specified; option intmethod(laplace)
      implied.

Mixed-effects probit regression                 Number of obs     =        360

        Grouping information




 
     No. of       Observations per group

 Group variable
 
     groups    Minimum    Average    Maximum

 
 
           _all
 
          1        360      360.0        360

         female
 
         60          6        6.0          6




Integration method: laplace

                                                Wald chi2(3)      =      40.12
Log likelihood = -208.14476                     Prob > chi2       =     0.0000
 ( 1)  [/]var(_cons[_all>male]) - [/]var(_cons[female]) = 0



           y 
 
 Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
 
 
 
       1.wsm 
 
  -.4130233   .2945617    -1.40   0.161    -.9903537    .1643071
       1.wsf 
 
  -1.720578   .3294304    -5.22   0.000     -2.36625   -1.074907
              
 
     wsm#wsf 
 
        1 1  
 
  -1.720578   .3294304    -5.22   0.000     -2.36625   -1.074907
              
 
       _cons 
 
   .5963019   .2514219     2.37   0.018      .103524     1.08908
 
 
 
_all>male    
 
   var(_cons)
 
   .4156124   .1500219                       .204849    .8432244
 
 
 
female       
 
   var(_cons)
 
   .4156124   .1500219                       .204849    .8432244

You can also fit our model using a logit model (see melogit) or a complementary log-log model (see mecloglog).

Also watch A tour of multilevel GLMs.