Cross level interaction - use of testparm after mixed vs manual calculation for statistical significance

[Emily Mann]

Hello readers of the Stata list,

I cannot find a previously posted answer to my question.

I have run 2 multilevel models (people nested within neighbourhoods); one without and one with an interaction.

mixed WEMWBStotal i.agecat09 i.income09 i.NHses09 if sex2==2 || habneigh2:, cov(unstr) var vce(r)
mixed WEMWBStotal i.agecat09 i.income09 i.NHses09 NHses09#i.income09 if sex2==2 || habneigh2:income09, cov(unstr) var vce(r)

I have used two methods to calculate statistical significance of the interaction:

1. a manual calculation using -2*log likelihood (model one) minus -2* loglikehood (model 2). I then take the difference in the degrees of freedom (22) and the value from the calculation and check a chi2 table to see if the value from the calculation is greater than the p=0.05 value. When I use this method the interaction is not statistically significant.

2. I then used the post estimation test testparm (in the hope of saving some manual labour!). So, I type testparm NHses09#i.income09 after after the mixed model. I notice that when I use this command there are two degrees of freedom less (20 verse 22) because the test doesn't seem to take the reference interaction into account (eg. 1.NHses09#1.income09). When I use this test, the result from test parm is statistically significant.

( 1) [WEMWBStotal]2.NHses09#2.income09 = 0
( 2) [WEMWBStotal]2.NHses09#3.income09 = 0
( 3) [WEMWBStotal]2.NHses09#4.income09 = 0
( 4) [WEMWBStotal]2.NHses09#5.income09 = 0
( 5) [WEMWBStotal]2.NHses09#9.income09 = 0
( 6) [WEMWBStotal]3.NHses09#2.income09 = 0
( 7) [WEMWBStotal]3.NHses09#3.income09 = 0
( 8) [WEMWBStotal]3.NHses09#4.income09 = 0
( 9) [WEMWBStotal]3.NHses09#5.income09 = 0
(10) [WEMWBStotal]3.NHses09#9.income09 = 0
(11) [WEMWBStotal]4.NHses09#2.income09 = 0
(12) [WEMWBStotal]4.NHses09#3.income09 = 0
(13) [WEMWBStotal]4.NHses09#4.income09 = 0
(14) [WEMWBStotal]4.NHses09#5.income09 = 0
(15) [WEMWBStotal]4.NHses09#9.income09 = 0
(16) [WEMWBStotal]5.NHses09#2.income09 = 0
(17) [WEMWBStotal]5.NHses09#3.income09 = 0
(18) [WEMWBStotal]5.NHses09#4.income09 = 0
(19) [WEMWBStotal]5.NHses09#5.income09 = 0
(20) [WEMWBStotal]5.NHses09#9.income09 = 0

chi2( 20) = 33.61
Prob > chi2 = 0.0289

Is there any way to have the reference categories included in the testparm command? If not, should I stick to the manual calculation?

My purpose for needing to make a statement about statistical significance of the interaction is so that I can write one line in a manuscript. While the significance of an interaction is not a major issue, as the interaction is meaningful in other ways, I must feel confident in the steps I'm undertaking to test for significance and justify what I'm writing to my supervisors.

I look forward to your feedback.
Emily

[Clyde Schechter]

Well, you don't show the details, so it is hard to be specific on where things are going wrong, but something is wrong with what you describe.

As best I can guess from seeing only the commands and not the full outputs, the likelihood-ratio calculation you are doing involves comparing a model that has both bottom level interactions of income09 with NHSES09 and a cross-level interaction of income09 with habneigh2 (aka a random slope), with a model that has neither. By contrast the -testparm- command deals only with the bottom level interactions. So your two tests are testing two different hypotheses, and there is no reason the df should be the same, nor the results the same.

That said, Stata is correct in not including the reference levels when it runs -testparm-: they are not represented in the model as variables and the removal of the other interaction terms only reduces the df by the 20 that you see. You should not want those reference categories included in the chi square test.

I'm guessing that your 22 df count for the likelihood ratio test is incorrect. I think it should be 21: 20 for the same interaction terms you show in the -testparm- output plus 1 for the random slope term.

Finally, there is no reason to hand-calculate the likelihood ratio test, and no reason to look the result up in a table. Use the -lrtest- command. See -help lrtest-.

Added: By the way, even when both are done precisely correctly and applied to the same hypothesis, they do not, in general, produce the same result. They are usually pretty close, but not identical. For an understanding of the difference between them, see https://stats.idre.ucla.edu/other/mu...andor-similar/.

[Dirk Enzmann]

Thanks, Clyde, for pointing to the UCLA page explaining the difference between LR, Wald, and Score tests!

Do you (or someone else) happen(s) to know how to how to perform the score tests for omitted variables as shown at the UCLA page?

[Emily Mann]

Thank you for your reply, Clyde.

The link clarifies much of the info I have read in a concise way. It's good to know that model fit results will be slightly different for each test. The UCLA Stata pages are excellent!

Yes, it would be fantastic to use the -lrtest- command! But this doesn't work with robust vce. The Stata documentation also makes note of this being due to taking clustering into account. In this study, neighbourhoods were sampled first. I would use the -lrtest- if there's another way.

. 
. lrtest income1 income2
LR test likely invalid for models with robust vce
r(498);

I thought there may be an issue with attempting to compare two models that are not nested .

So I'm confused as to what direction I should now take to establish if the cross-level interaction is significant or not.... I guess I have two options.

1. Use -testparm- to run a Wald test on the as I've done on:

mixed WEMWBStotal i.agecat09 i.income09 i.NHses09 NHses09#i.income09 if sex2==2 || habneigh2:income09, cov(unstr) var vce(r)

But it seems you are saying that this won't test cross level interactions? I was of the understanding that neighbourhood ses interaction with the individual income would be considered a cross level interaction.

2. Do the manual lrtest again, but run:

Model 1: mixed WEMWBStotal i.agecat09 i.income09 i.NHses09 if sex2==2 || habneigh2:income09, cov(unstr) var vce(r)
Model 2: mixed WEMWBStotal i.agecat09 i.income09 i.NHses09 NHses09#i.income09 if sex2==2 || habneigh2:income09, cov(unstr) var vce(r)

I have included calculations not shown in post 1 to show how I arrived at 22df, not 21 as you suggest:

The calculations (log likelihoods are from the regression output below):
. di (-2*-14100.44) - (-2*-14087.58)
25.72

df (from estat IC info below) = 38 (model 2 #df) - 16 (model 1 #df) = 22

. *income
. mixed WEMWBStotal i.agecat09 i.income09 i.NHses09 if sex2==2 || habneigh2:, cov(unstr) var vce(r)
Note: single-variable random-effects specification in habneigh2 equation; covariance structure set to identity

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log pseudolikelihood = -14101.266  
Iteration 1:   log pseudolikelihood = -14100.438  
Iteration 2:   log pseudolikelihood = -14100.437  

Computing standard errors:

Mixed-effects regression                        Number of obs     =      4,006
Group variable: habneigh2                       Number of groups  =        200

                                                Obs per group:
                                                              min =          3
                                                              avg =       20.0
                                                              max =         66

                                                Wald chi2(13)     =     220.93
Log pseudolikelihood = -14100.437               Prob > chi2       =     0.0000

                                                        (Std. Err. adjusted for 200 clusters in habneigh2)
----------------------------------------------------------------------------------------------------------
                                         |               Robust
                             WEMWBStotal |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-----------------------------------------+----------------------------------------------------------------
                                agecat09 |
                                  47-51  |   .2241137   .3694646     0.61   0.544    -.5000236     .948251
                                  52-56  |   1.568205   .3794711     4.13   0.000     .8244551    2.311954
                                  57-61  |   2.419533   .4224282     5.73   0.000     1.591589    3.247477
                                  62-70  |   4.070298   .3921503    10.38   0.000     3.301698    4.838898
                                         |
                                income09 |
                          $72800-129999  |  -1.011768   .3772156    -2.68   0.007    -1.751097    -.272439
                           $52000-72799  |  -2.383174   .4482803    -5.32   0.000    -3.261787    -1.50456
                           $26000-51599  |  -2.320487   .5078401    -4.57   0.000    -3.315835   -1.325139
                       Less than $25999  |  -4.509945   .5602949    -8.05   0.000    -5.608103   -3.411787
missing/Don't want to answer/Don't know  |  -1.590864     .45103    -3.53   0.000    -2.474866   -.7068611
                                         |
                                 NHses09 |
                                     Q2  |  -.0367962   .3561964    -0.10   0.918    -.7349283    .6613358
                                     Q3  |  -.5921428   .4517334    -1.31   0.190    -1.477524    .2932385
                                     Q4  |  -.7473764   .3921647    -1.91   0.057    -1.516005    .0212523
                           Q5(most dis)  |  -1.825169   .5509012    -3.31   0.001    -2.904916   -.7454227
                                         |
                                   _cons |   51.76537    .383879   134.85   0.000     51.01298    52.51776
----------------------------------------------------------------------------------------------------------

------------------------------------------------------------------------------
                             |               Robust           
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
habneigh2: Identity          |
                  var(_cons) |   .4120262   .3686788      .0713303    2.379992
-----------------------------+------------------------------------------------
               var(Residual) |   66.42289   1.850917      62.89245    70.15151
------------------------------------------------------------------------------

. estat ic

Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |      4,006         .  -14100.44      16    28232.87    28333.6
-----------------------------------------------------------------------------
               Note: N=Obs used in calculating BIC; see [R] BIC note.

. 
. mixed WEMWBStotal i.agecat09 i.income09 i.NHses09 NHses09#i.income09 if sex2==2 || habneigh2:income09, cov(unstr) var
>  vce(r)

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log pseudolikelihood = -14091.163  
Iteration 1:   log pseudolikelihood = -14087.624  (not concave)
Iteration 2:   log pseudolikelihood = -14087.621  (backed up)
Iteration 3:   log pseudolikelihood = -14087.579  
Iteration 4:   log pseudolikelihood = -14087.579  

Computing standard errors:

Mixed-effects regression                        Number of obs     =      4,006
Group variable: habneigh2                       Number of groups  =        200

                                                Obs per group:
                                                              min =          3
                                                              avg =       20.0
                                                              max =         66

                                                Wald chi2(33)     =     295.48
Log pseudolikelihood = -14087.579               Prob > chi2       =     0.0000

                                                           (Std. Err. adjusted for 200 clusters in habneigh2)
-------------------------------------------------------------------------------------------------------------
                                            |               Robust
                                WEMWBStotal |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------------------------------------+----------------------------------------------------------------
                                   agecat09 |
                                     47-51  |   .2617741   .3718799     0.70   0.481    -.4670971    .9906453
                                     52-56  |    1.54549   .3840106     4.02   0.000     .7928427    2.298137
                                     57-61  |   2.421821   .4198101     5.77   0.000     1.599009    3.244634
                                     62-70  |   4.048828   .3873469    10.45   0.000     3.289642    4.808014
                                            |
                                   income09 |
                             $72800-129999  |  -.6872012   .6021817    -1.14   0.254    -1.867456    .4930532
                              $52000-72799  |  -2.845769   .8373992    -3.40   0.001    -4.487042   -1.204497
                              $26000-51599  |  -1.782084   .8521309    -2.09   0.036    -3.452229   -.1119378
                          Less than $25999  |  -4.444242   1.486314    -2.99   0.003    -7.357363   -1.531121
   missing/Don't want to answer/Don't know  |  -1.458451   .7301205    -2.00   0.046    -2.889461   -.0274411
                                            |
                                    NHses09 |
                                        Q2  |   .1859142   .7209583     0.26   0.797    -1.227138    1.598966
                                        Q3  |  -.9465613   .8573326    -1.10   0.270    -2.626902    .7337798
                                        Q4  |  -1.238881   1.177055    -1.05   0.293    -3.545866    1.068104
                              Q5(most dis)  |   2.061356   1.349126     1.53   0.127    -.5828816    4.705594
                                            |
                           NHses09#income09 |
                          Q2#$72800-129999  |  -.8118072   .9394585    -0.86   0.388    -2.653112    1.029498
                           Q2#$52000-72799  |  -.4672019   1.161241    -0.40   0.687    -2.743192    1.808789
                           Q2#$26000-51599  |  -.1866604   1.251139    -0.15   0.881    -2.638849    2.265528
                       Q2#Less than $25999  |   .5880443   1.815593     0.32   0.746    -2.970453    4.146541
Q2#missing/Don't want to answer/Don't know  |   .3177766     1.0559     0.30   0.763    -1.751748    2.387302
                          Q3#$72800-129999  |    .504742   1.015798     0.50   0.619    -1.486186     2.49567
                           Q3#$52000-72799  |   1.008761   1.385297     0.73   0.466    -1.706371    3.723893
                           Q3#$26000-51599  |  -.8607832   1.374017    -0.63   0.531    -3.553807     1.83224
                       Q3#Less than $25999  |   1.026593   1.822654     0.56   0.573    -2.545743     4.59893
Q3#missing/Don't want to answer/Don't know  |   .9819355   1.176457     0.83   0.404    -1.323877    3.287748
                          Q4#$72800-129999  |  -.1457969   1.293291    -0.11   0.910      -2.6806    2.389006
                           Q4#$52000-72799  |   1.789663   1.384368     1.29   0.196    -.9236486    4.502975
                           Q4#$26000-51599  |   .3499844   1.602927     0.22   0.827    -2.791696    3.491665
                       Q4#Less than $25999  |   1.044397   1.973897     0.53   0.597     -2.82437    4.913165
Q4#missing/Don't want to answer/Don't know  |    .057093   1.571414     0.04   0.971    -3.022822    3.137008
                Q5(most dis)#$72800-129999  |   -3.57381   1.714341    -2.08   0.037    -6.933857   -.2137634
                 Q5(most dis)#$52000-72799  |  -2.111637   1.727934    -1.22   0.222    -5.498325    1.275051
                 Q5(most dis)#$26000-51599  |  -4.435635   1.703615    -2.60   0.009    -7.774658   -1.096612
             Q5(most dis)#Less than $25999  |  -4.887929   2.076589    -2.35   0.019    -8.957969   -.8178894
                              Q5(most dis) #|
   missing/Don't want to answer/Don't know  |  -4.612908   1.864828    -2.47   0.013    -8.267903   -.9579129
                                            |
                                      _cons |   51.62718   .4678249   110.36   0.000     50.71026     52.5441
-------------------------------------------------------------------------------------------------------------

------------------------------------------------------------------------------
                             |               Robust           
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
habneigh2: Unstructured      |
               var(income09) |   .0304079   .0183549       .009315    .0992638
                  var(_cons) |    .009237   .0005131      .0082842    .0102994
         cov(income09,_cons) |   .0167594   .0054833      .0060123    .0275065
-----------------------------+------------------------------------------------
               var(Residual) |   65.67563   1.809054      62.22396    69.31876
------------------------------------------------------------------------------

. estat ic

Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |        Obs  ll(null)  ll(model)      df         AIC        BIC
-------------+---------------------------------------------------------------
           . |      4,006         .  -14087.58      38    28251.16   28490.39
-----------------------------------------------------------------------------
               Note: N=Obs used in calculating BIC; see [R] BIC note.

. 
. testparm NHses09#i.income09

 ( 1)  [WEMWBStotal]2.NHses09#2.income09 = 0
 ( 2)  [WEMWBStotal]2.NHses09#3.income09 = 0
 ( 3)  [WEMWBStotal]2.NHses09#4.income09 = 0
 ( 4)  [WEMWBStotal]2.NHses09#5.income09 = 0
 ( 5)  [WEMWBStotal]2.NHses09#9.income09 = 0
 ( 6)  [WEMWBStotal]3.NHses09#2.income09 = 0
 ( 7)  [WEMWBStotal]3.NHses09#3.income09 = 0
 ( 8)  [WEMWBStotal]3.NHses09#4.income09 = 0
 ( 9)  [WEMWBStotal]3.NHses09#5.income09 = 0
 (10)  [WEMWBStotal]3.NHses09#9.income09 = 0
 (11)  [WEMWBStotal]4.NHses09#2.income09 = 0
 (12)  [WEMWBStotal]4.NHses09#3.income09 = 0
 (13)  [WEMWBStotal]4.NHses09#4.income09 = 0
 (14)  [WEMWBStotal]4.NHses09#5.income09 = 0
 (15)  [WEMWBStotal]4.NHses09#9.income09 = 0
 (16)  [WEMWBStotal]5.NHses09#2.income09 = 0
 (17)  [WEMWBStotal]5.NHses09#3.income09 = 0
 (18)  [WEMWBStotal]5.NHses09#4.income09 = 0
 (19)  [WEMWBStotal]5.NHses09#5.income09 = 0
 (20)  [WEMWBStotal]5.NHses09#9.income09 = 0

           chi2( 20) =   33.61
         Prob > chi2 =    0.0289

[Clyde Schechter]

OK, your df calculation was correct. I had overlooked the fact that when you remove the random slope (cross-level interaction) you also lose the covariance between the slope and intercept, so the difference really is 22 df as you say. That said, you are also correct that the likelihood-ratio test cannot be used with robust vce. But it is actually the likelihood ratio test that cannot be used; it's not just that the -lrtest- command is being uncooperative. -lrtest- won't run because the likelihood ratio test itself is not permissible here. So your hand calculated likelihood ratio test is invalid in any case. Sorry I hadn't paid attention to that vce(robust) before.

I do not know of any way to use -testparm- to include the cross-level interaction. -testparm- is giving you correct results for the bottom-level interactions, but that is the best it can do for you.

So we face the question of why you want to do this test in the first place. If testing this hypothesis is one of your research goals (and, I would have to say that it is a strange composite hypothesis, but a perfectly legal one), then I would say give up the robust VCE, rerun the models without it, and just do the likelihood ratio test comparing the model containing all of the interactions with the model containing none of them. Those are nested models, and if you use ordinary VCE instead of robust, -lrtest- will give you your answer. (I'm not sure what vce(robust) is getting you anyway. If your worry is within-group error correlation, vce(robust) doesn't fix that. If your worry is heteroscedasticity, -mixed- allows you to estimate group-specific heteroscedastic residuals. See the -by()- suboption of the -residuals()- option in the manual section for -mixed.)

If this is about model selection, I will spare you my rant on why I believe model selection should not be based on summary statistics or tests and just point out that, given that you are going to do that, AIC or BIC would probably be your best bet.

[Clyde Schechter]

Wait, I see you have an even bigger problem in your model. income09 is actually a categorical variable in the fixed-effects level of your model. So putting a random slope on income09 is not valid, because that treats it as a continuous variable. Now, you can, if you like, create hand-made indicators for the categories of income09 and include a random slope on (all but one of) them. But you can't have income09 acting as discrete in the bottom level and continuous in the top level of your model.

Digging deeper, you probably should reconsider treating non-response to income as a separate category of the income variable. It is well known that doing this biases the results you get for the other categories. You should consider other ways of dealing with this missing data problem.

[John Moran]

#3, Dirk Enzmann query

I think that "scoregrp", SJ 13-4, st0321 may do what you want

john moran

[Emily Mann]

This Q&A is good practice for me. I'll do my best with my reply as I'm new to regression and multilevel regression and following steps suggested by my supervisors (who do use Stata but normally do multilevel modelling in MLWin). There are some things that are not making sense to me, so this is why I am asking the forum for help while they are away.

1. It makes sense what you say about the LR test and -testparm-

2. Why cross level interactions, and therefore why testing this hypothesis? It's not to check on model selection (yes I have read many of your posts on this). I am examining several individual SES measures (education, occupation and income) in separate linear regression models. These results are being presented for each. Then we use a multilevel model to see if the neighbourhood SES contributes to change in the dependent variable (mental well-being), these results are also presented. The cross-level interaction is one step further to see if there is a 'double disadvantage effect' i.e. do people with the lowest level of education & live in the most disadvantaged areas have worse mental well-being compared than those with the highest level of education, who live in the with the most advantaged neighbourhood. I have used -margins- to calculate the predicted means to plot on bar chart to demonstrate these trends (will probably have to revisit this given your comment I address in point 4). However, there still seems to be this need to report in the results section of a paper one line stating if any of the cross level interactions are statistically significant.

3. I don't really understand why vce robust is being used. I was told it's because the data are clustered. If I took it away from the model looking at the cross level interactions, then would I have to remove it from my single level regression models?

4. Picking up the problem with income09 as a categorical is very helpful! So this only affects the cross-level interaction models. In the example I've used there are 5 income categories, not including missing. You are saying I'd have to create 5 new income variables, one for each existing category. Then would I run one model with each new income variable at the bottom and at the top, or do I put each new variables at the bottom and top in the one model?

5. I did have the missing/don't want to answer/don't know for income split into separate categories, but then was told to combine them because if left separate some of the cell sizes would be too small (less than 15) to assess the interaction between neighbourhood and education/occupation/income. I'm not sure what are the more important factors to consider my decision. Maybe I have to split out again.

6. I've had many roundabout conversations about missing data. There seems to be many differing views about how to deal with these (a) to drop from the sample if less than 5%, (b) impute, or (c) leave in and not report the stats for missing. I understand that many statisticians, as you suggest, dislike leaving missing data in due to bias. We know the data are missing at random. The missing for education and income are both less than 5%. The way forward at this point is to leave them all in, but not report the missing data in the tables (making a note of this in the results of course). The main reason is to keep as many people in the analysis as possible.

If you know of any good links about cross level interactions with examples in Stata please feel free to post them I've looked in the Multilevel Stata text book but didn't find it all that helpful.
The other useful resource I'l have to rewatch is the You Tube video (not sure if I should post link here) on "HLM testing and probing within and cross level interactions using Stata" by Mike Crowson.

[Dirk Enzmann]

#7 (John Moran's response): No, unfortunately not - scoregrp does not allow to add variables to the model (it requires that all variables specified belong to the list of predictors).

[Clyde Schechter]

Responding to #6 to the extent I can.

2. I think you misunderstood me. I was not questioning why you are exploring cross-level interactions. That makes perfect sense to me in this context. I was simply observing that simultaneously testing for cross-level and within-bottom-level interactions is unusual. It's an odd conjunction of hypotheses that are typically tested separately. That said, there is nothing inherently wrong with it.

3. Use of -vce(robust)- does not help with clustering. There is a separate -vce(cluster cluster_variable)- option for that. (But it, too, invalidates likelihood ratio tests. I think there may be some confusion here because, starting with, I think, version 13, Stata allows people to specify -vce(robust)- with -xtreg, fe- (and a couple of other fixed-effects estimators) and it is automatically converted to -vce(cluster panelid)-, where panelid is the panel variable from the -xtset- command. So this may cause some people to confuse robust and cluster-robust variance estimators. But they are different, and the automatic substitution of cluster robust for robust applies only to a limited number of commands, -mixed- not being one of them.

4. The problem here is that you need income09 to be treated the same way in the fixed and random effects parts of the model. Unfortunately, Stata does not support factor variable notation in the random effects part. So you can't just write -mixed ... i.income09 ... || habneigh2: i.income09...- That would be ideal, but isn't supported. So this is one of those situations where factor variable notation fails and you have to "roll your own" indicator variables. You could do something like this:

tag income09, gen(income_indicators)
drop income_indicators1
mixed ... income_indicators*... || habneigh2: income_indicators*...

You could also dig out the old -xi- to do this as:

xi: mixed ... i.income09 || habneigh2: i.income09...

5. My concern here is not that you have combined missing/ don't want to answer/ don't know into a single category rather than keeping them separate. My concern is that you have assigned them to any categories at all. These are not known values of income. They are non-response, for various reasons. The actual income distribution among the non-responders is not known, but it is presumably some mixture of the actual income values that are available for the responders. So, let's say somebody has a real income that would, if he or she responded, put that person in category 2 of income09. But that person doesn't respond, so is now in category 9. The point is that category 9 is an unknown mixture of people who would belong in the other categories if we had that information. They are probably not a homogeneous group, and they have removed a biased sample (because there is likely something peculiar about the non-responders that may well be associated with their actual incomes) from the other categories. The result is that none of these indicators actually represents what it purports to represent, so that results on this entire variable are potentially biased.

6. Yes, discussions about missing data tend to be roundabout. Despite the current popularity of multiple imputation, in the real world there really aren't good solutions to the missing data problem because missingness at random is pretty unusual. I know you say that you know your data to be missing at random. But I wonder how you know that. Given the nature of the data you have disclosed, I would be willing to make at least a modest, if not a large, wager that they are not missing at random, and a larger one that you do not in fact "know" that they are, even if they really are. That's because it is difficult to ever know whether your data are missing at random or not--and, in particular, there is no way to know this based on the data themselves: you can only know this from an external study of the missingness mechanism, and a deep knowledge of how the variables with missing values relate to other variables in your data. My own inclination, with less than 5% missing data would be to just analyze complete cases and hope for the best,* and perhaps do some kind of robustness analysis. If I could persuade myself that missingness at random is at least plausible, then I might also do multiple imputation and see if things come out more or less the same (though I would not know what to believe if they came out appreciably different.) And frankly, with a socially sensitive variable like income, I would have a hard time believing that its missingness was actually at random unless I at least had a good hard look at the missingness process, or had several other variables in the analysis that are good predictors of income.

*I know that many people would vehemently disagree with this recommendation, and I respect their contrary opinions on this matter. It is a recommendation I make with a heavy heart--I agree that a complete cases analysis with 5% missing data runs a serious risk of having serious bias. It's not a good approach. It's just that I don't really believe that the alternatives are really any better. In my line of work, at least, missingness at random is uncommon, and the uncritical use of multiple imputation is, in my view, magical thinking. It creates the illusion, but not the reality, of dealing with the problem. I prefer to have problems be obvious, not hidden.

[Emily Mann]

simultaneously testing for cross-level and within-bottom-level interactions is unusual.

Ah okay. I see what you are saying. I didn't realise I was testing interactions both within and across levels. I thought I was only testing the top level. Regardless, I also didn't know that testing both at once is unusual practice. Now I do.

After pulling a few hairs out, I realised the code should be:

tab income09, gen(income_indicators)

not

tag income09, gen(income_indicators)

Though something is not right with what I've done next. I also omitted -var vce(r)- from the code as it sounds like it's not necessary, unless I'm interpreting you wrong:

. mixed WEMWBStotal i.agecat09 income09 i.NHses09 NHses09#income_indicators* if sex2==2 || habneigh2:income_indicators*, cov(unstr)

note: 5.NHses09#1.income_indicators6 omitted because of collinearity

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log likelihood = -14097.752  
Iteration 1:   log likelihood = -14083.779  
Iteration 2:   log likelihood = -14083.314  
Hessian is not negative semidefinite
r(430);

In reply to points 5 &6
I take your point about non-response being a non-response regardless if that person ticks 'don't know' on a survey or ticks nothing at all, and that this can bias results. I was under the impression that it was 'true' missing--i.e. people did not fill in any option on the survey, verses someone that picked 'don't know'--that I need to be concerned about and then 'treat' the 'true missing' non-response data in a particular way.

If I look at non-response in total, rather than 'true missing', then I have about 10-12% non-response for occupation and income. Education is at about 2% missing.

My mistake, we ran a sensitivity analysis using logistic regression to test 'responders to 14 items' vs 'responders to 0-13 items' of the 14 item scale used to calculate the outcome variable. Cases were dropped who did not fill in all 14 questions. The regression showed that people with low SES (using occupation and income) were not were answering all 14 items.

As for the SES predictors, I'm not sure if they are MAR, MNAR or MCAR. Though yes you are most probably right that there will be some social patterning related to the non-response of these indicators and therefore won't be MAR.

So given I can't play the <%5 missing card for occupation and income, I suppose it's back to the drawing board about how to 'treat' non-response items for occupation and income.

[Clyde Schechter]

Sorry about the -tab- vs -tag- thing; that was purely a typo.

For your -mixed- command it should look like this:

mixed WEMWBStotal i.agecat09  NHses09##income_indicators* if sex2==2 || habneigh2:income_indicators*, cov(unstr)

Where you went wrong was including the old income09 variable. First, you included it as continuous, and second, it doesn't belong there anyway because income_indicators* has now taken over its role throughout the code for this command. Note also the use of the ## operator (not #). This makes life much simpler because you don't have to write out all of the individual income_indicator variables, and you don't have to write out NHses09 again.

[Emily Mann]

No worries.

mmm apparently that code doesn't work either...
Does it matter that NHses09 is a factor variable?

. mixed WEMWBStotal i.agecat09 NHses09##income_indicators* if sex2==2 || habneigh2:income_indicators*, cov(unstr)

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -14097.752  
Iteration 1:   log likelihood = -14083.779  
Iteration 2:   log likelihood = -14083.314  
numerical derivatives are approximate
nearby values are missing
Iteration 3:   log likelihood = -14083.284  
numerical derivatives are approximate
nearby values are missing
Hessian is not negative semidefinite
r(430);

[Clyde Schechter]

Well, you can try seeing if adding the -difficult- option to your command helps. But I think the problem is that you are cutting the data too fine. Looking at your earlier outputs, I see you have a sample size of about 4,000 which is divided into about 200 neighborhoods, each of which has, on average, about 20 observations. Now, by trying to estimate a neighborhood # income interaction, you are, in effect, subdividing each neighborhood into 5 (or 6 if you're still using the 9 for non-response category) subsets defined by the income groups. So on average each of those will have only 4 observations. But it's worse than that because it's not as if every neighborhood as exactly 20 observations. Some have more, and some have fewer, as few as 3. It follows that at least some of these combinations of neighborhood and income group will have no observations at all, and many will have only a handful. That is the likely cause of the problem. Your data set is simply not large enough to withstand this kind of fine-grained subdivision.

In general, my advice would be to eliminate the cross-level interaction, but since that is a primary purpose of your study, I won't advise you to go there just yet. But you do need, then, to coarsen your income variable. Why don't you combine some adjacent income categories. Get rid of the 9 (non-response) category altogether if you haven't already. And then combine some others to reduce the number of categories to 2 (probably splitting at or near the median). Then try to estimate your model. If that works OK and you would like to be a bit more fine-grained, try re-organizing the income variable into 3 groups, and see if you can get away with that. (You probably won't, but nothing to lose by trying.)

If that doesn't work, you could also look into consolidating some neighborhoods to reduce the number of those. This gets more complicated because you would have to decide which neighborhoods are sufficiently similar in whatever are the relevant characteristics that it is reasonable to combine them. Let's just hope it doesn't come to that.

[Emily Mann]

Because we are interested in the quartiles of NHses, which already reduces the 200 neighbourhoods into 5 categories, I thought we'd already have the issue of small cell size covered. Is there any reason why NHses can't go in the position of habneigh2? Or does habneigh2 have to remain at the top level in the code because it's the initial sampling frame and most fine grain option available.

No luck with the -difficult- option for the existing code with all the categorical variables, even after I dropped the missing.

I then tried dichotomous vars for education and income, and tri for occupation which I already had. Income worked, but there are issues with education and occupation - the random effect CIs don't look quite right/SE non existent (see output below).

I don't know why income worked though because as you can see the numbers are very small or 0 for example, for income, <$26,000, and even for $26000 and above. I know you suggested split at the median, but this split is based on that we are most interested in what's happening in the most disadvantaged group. At the moment I can't think of any substantive reason as to why splitting income in half is meaningful.

. tab habneigh2 incomedi_indicators if sex2==2

   698:w2: |
   HABITAT |  incomedi09==$26000
neighbourh |       and above
ood number |         0          1 |     Total
-----------+----------------------+----------
         1 |         0          3 |         3 
         2 |         3         13 |        16 
         3 |         3          3 |         6 
         4 |         4         10 |        14 
         5 |         5          5 |        10 
         6 |         1          8 |         9 
         7 |         4          8 |        12 
         8 |         2          5 |         7 
         9 |         4          7 |        11 
        10 |         2          9 |        11 
        11 |         5         19 |        24 
        12 |         4          9 |        13 
        13 |         2         10 |        12 
        14 |         2          8 |        10 
        15 |         8          6 |        14 
        16 |         7          7 |        14 
        17 |         3          9 |        12 
        18 |         5         13 |        18 
        19 |         1          7 |         8 
        20 |         2         12 |        14 
        21 |         3         10 |        13 
        22 |         0         10 |        10 
        23 |         2         10 |        12 
        24 |         2         17 |        19 
        25 |         3          7 |        10 
        26 |         5          9 |        14 
        27 |         4         12 |        16 
        28 |         3          9 |        12 
        29 |         1         13 |        14 
        30 |         2          4 |         6 
        31 |         3          4 |         7 
        32 |         3         10 |        13 
        33 |         4          9 |        13 
        34 |         6          5 |        11 
        35 |         5          9 |        14 
        36 |         5         10 |        15 
        37 |         2          6 |         8 
        38 |         3         11 |        14 
        39 |         4          6 |        10 
        40 |         1         14 |        15 
        41 |        13         27 |        40 
        42 |         7         14 |        21 
        43 |         4         13 |        17 
        44 |         3          8 |        11 
        45 |         4          9 |        13 
        46 |         5         22 |        27 
        47 |         6         24 |        30 
        48 |         6         13 |        19 
        49 |         4          5 |         9 
        50 |         4         13 |        17 
        51 |         1          6 |         7 
        52 |         3          4 |         7 
        53 |         6          4 |        10 
        54 |         5         20 |        25 
        55 |         3         16 |        19 
        56 |         5          7 |        12 
        57 |         4         11 |        15 
        58 |         4         11 |        15 
        59 |         3          4 |         7 
        60 |         4         12 |        16 
        61 |         3         10 |        13 
        62 |         2          9 |        11 
        63 |        11         25 |        36 
        64 |         5         16 |        21 
        65 |         2          8 |        10 
        66 |         6         21 |        27 
        67 |         4         21 |        25 
        68 |         2         13 |        15 
        69 |         2         23 |        25 
        70 |         4         13 |        17 
        71 |         4         14 |        18 
        72 |         9         26 |        35 
        73 |         2          7 |         9 
        74 |         4         10 |        14 
        75 |         3         17 |        20 
        76 |         7         14 |        21 
        77 |         4          7 |        11 
        78 |         7          9 |        16 
        79 |         3          5 |         8 
        80 |         7         10 |        17 
        81 |         7          9 |        16 
        82 |         6         20 |        26 
        83 |         5          9 |        14 
        84 |        14         31 |        45 
        85 |         4          7 |        11 
        86 |         2         12 |        14 
        87 |         6          9 |        15 
        88 |         8         12 |        20 
        89 |         4         10 |        14 
        90 |        14         24 |        38 
        91 |         5         10 |        15 
        92 |         7         11 |        18 
        93 |         4         11 |        15 
        94 |         5         20 |        25 
        95 |        16         32 |        48 
        96 |         4         18 |        22 
        97 |         7         10 |        17 
        98 |         3         13 |        16 
        99 |         7         16 |        23 
       100 |        12         18 |        30 
       101 |         4         11 |        15 
       102 |        12         15 |        27 
       103 |         7         18 |        25 
       104 |        14         42 |        56 
       105 |         9          7 |        16 
       106 |         3          6 |         9 
       107 |         5         10 |        15 
       108 |        14         35 |        49 
       109 |         3         14 |        17 
       110 |         6         12 |        18 
       111 |        18         46 |        64 
       112 |         2         25 |        27 
       113 |         4          9 |        13 
       114 |        11         54 |        65 
       115 |        11         12 |        23 
       116 |         5         11 |        16 
       117 |         7         23 |        30 
       118 |         8          5 |        13 
       119 |        10         41 |        51 
       120 |         5          9 |        14 
       121 |        13         43 |        56 
       122 |         3         22 |        25 
       123 |         1          8 |         9 
       124 |         5         10 |        15 
       125 |         3         10 |        13 
       126 |         8         12 |        20 
       127 |        12         37 |        49 
       128 |         6         23 |        29 
       129 |        10         12 |        22 
       130 |        10         25 |        35 
       131 |         9         26 |        35 
       132 |        15         16 |        31 
       133 |         7         21 |        28 
       134 |        12         20 |        32 
       135 |        13         24 |        37 
       136 |         6         11 |        17 
       137 |         7          0 |         7 
       138 |         6         21 |        27 
       139 |         8         17 |        25 
       140 |        14         12 |        26 
       141 |         4          4 |         8 
       142 |        13         11 |        24 
       143 |        15         40 |        55 
       144 |        19          6 |        25 
       145 |         7         13 |        20 
       146 |        10         20 |        30 
       147 |         6          3 |         9 
       148 |         2         10 |        12 
       149 |         8          9 |        17 
       150 |         5         10 |        15 
       151 |        16         33 |        49 
       152 |        17         48 |        65 
       153 |         6         25 |        31 
       154 |         8         11 |        19 
       155 |         3         12 |        15 
       156 |         8         20 |        28 
       157 |         6         18 |        24 
       158 |         3          5 |         8 
       159 |         5          3 |         8 
       160 |         3          3 |         6 
       161 |        12         20 |        32 
       162 |        12         11 |        23 
       163 |        14         37 |        51 
       164 |        12         19 |        31 
       165 |         7         21 |        28 
       166 |        13         13 |        26 
       167 |        11         15 |        26 
       168 |         4          3 |         7 
       169 |         2          8 |        10 
       170 |         9          8 |        17 
       171 |         2          9 |        11 
       172 |         8          6 |        14 
       173 |         5          4 |         9 
       174 |         4         10 |        14 
       175 |        18         30 |        48 
       176 |        14         56 |        70 
       177 |         6         12 |        18 
       178 |         3         19 |        22 
       179 |         8         22 |        30 
       180 |        13         29 |        42 
       181 |        11          3 |        14 
       182 |         8         16 |        24 
       183 |         8         22 |        30 
       184 |         6         12 |        18 
       185 |         9         14 |        23 
       186 |        11         19 |        30 
       187 |         6         24 |        30 
       188 |         8         20 |        28 
       189 |         7         18 |        25 
       190 |        10          7 |        17 
       191 |        10         14 |        24 
       192 |        21         32 |        53 
       193 |         4         23 |        27 
       194 |         2          6 |         8 
       195 |         2         10 |        12 
       196 |         5         14 |        19 
       197 |         6         15 |        21 
       198 |         1         10 |        11 
       199 |         3          7 |        10 
       200 |         9         17 |        26 
-----------+----------------------+----------
     Total |     1,267      2,913 |     4,180

. mixed WEMWBStotal i.agecat09 NHses09##incomedi_indicators* if sex2==2 || habneigh2:incomedi_indicators*, cov(unstr) diff

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -14127.921  
Iteration 1:   log likelihood = -14124.344  
Iteration 2:   log likelihood = -14124.227  
Iteration 3:   log likelihood = -14124.227  
Iteration 4:   log likelihood = -14124.227  
Iteration 5:   log likelihood = -14124.227  
Iteration 6:   log likelihood = -14124.227  

Computing standard errors:

Mixed-effects ML regression                     Number of obs     =      4,006
Group variable: habneigh2                       Number of groups  =        200

                                                Obs per group:
                                                              min =          3
                                                              avg =       20.0
                                                              max =         66

                                                Wald chi2(13)     =     130.12
Log likelihood = -14124.227                     Prob > chi2       =     0.0000

----------------------------------------------------------------------------------------------
                 WEMWBStotal |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-----------------------------+----------------------------------------------------------------
                    agecat09 |
                      47-51  |   .2087054   .4108501     0.51   0.611     -.596546    1.013957
                      52-56  |   1.445663   .4115125     3.51   0.000     .6391135    2.252213
                      57-61  |   2.056148   .4170483     4.93   0.000     1.238748    2.873548
                      62-70  |   3.328804   .4316708     7.71   0.000     2.482745    4.174863
                             |
                     NHses09 |
                         Q2  |   .3868772   .8945895     0.43   0.665    -1.366486     2.14024
                         Q3  |  -.4671367   .8831104    -0.53   0.597    -2.198001    1.263728
                         Q4  |  -1.237773     .84147    -1.47   0.141    -2.887024    .4114782
               Q5(most dis)  |  -3.833602   .8707006    -4.40   0.000    -5.540144   -2.127061
                             |
      1.incomedi_indicators2 |   1.174813   .6545352     1.79   0.073     -.108052    2.457679
                             |
NHses09#incomedi_indicators2 |
                       Q2#1  |  -.8350251   .9482119    -0.88   0.379    -2.693486    1.023436
                       Q3#1  |  -.7531751   .9525632    -0.79   0.429    -2.620165    1.113814
                       Q4#1  |  -.0083981   .9236235    -0.01   0.993    -1.818667    1.801871
             Q5(most dis)#1  |   2.073114   .9882831     2.10   0.036     .1361143    4.010113
                             |
                       _cons |   49.69618   .6742557    73.71   0.000     48.37467     51.0177
----------------------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
habneigh2: Unstructured      |
               var(income..) |   1.563727   1.164654      .3632381    6.731787
                  var(_cons) |   2.927282   1.346526       1.18828    7.211246
         cov(income..,_cons) |  -2.139502   1.237891     -4.565724    .2867199
-----------------------------+------------------------------------------------
               var(Residual) |   66.74843   1.524085      63.82714    69.80344
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 9.01                  Prob > chi2 = 0.0291

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

. mixed WEMWBStotal i.agecat09 NHses09##edcatdi_indicators* if sex2==2 || habneigh2:edcatdi_indicators*, cov(unstr) diff

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -14133.148  
Iteration 1:   log likelihood =  -14129.81  
Iteration 2:   log likelihood = -14129.759  
Iteration 3:   log likelihood =  -14129.74  
Iteration 4:   log likelihood =  -14129.74  

Computing standard errors:

Mixed-effects ML regression                     Number of obs     =      4,006
Group variable: habneigh2                       Number of groups  =        200

                                                Obs per group:
                                                              min =          3
                                                              avg =       20.0
                                                              max =         66

                                                Wald chi2(13)     =     134.67
Log likelihood =  -14129.74                     Prob > chi2       =     0.0000

---------------------------------------------------------------------------------------------
                WEMWBStotal |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
----------------------------+----------------------------------------------------------------
                   agecat09 |
                     47-51  |   .2282616   .4126152     0.55   0.580    -.5804493    1.036973
                     52-56  |   1.414073   .4125196     3.43   0.001     .6055496    2.222597
                     57-61  |   2.014854   .4170017     4.83   0.000     1.197546    2.832163
                     62-70  |   3.207187    .426782     7.51   0.000      2.37071    4.043664
                            |
                    NHses09 |
                        Q2  |  -.6247495   .6299482    -0.99   0.321    -1.859425    .6099263
                        Q3  |  -.6553735   .6706877    -0.98   0.328    -1.969897    .6591503
                        Q4  |  -.8005429   .7229341    -1.11   0.268    -2.217468    .6163819
              Q5(most dis)  |  -2.179237   .8580179    -2.54   0.011    -3.860922   -.4975532
                            |
      1.edcatdi_indicators2 |  -1.229189   .5225625    -2.35   0.019    -2.253393   -.2049853
                            |
NHses09#edcatdi_indicators2 |
                      Q2#1  |   .7116481   .7863735     0.90   0.365    -.8296156    2.252912
                      Q3#1  |  -.4351626   .8235036    -0.53   0.597      -2.0492    1.178875
                      Q4#1  |   -.486945   .8617599    -0.57   0.572    -2.175963    1.202073
            Q5(most dis)#1  |   -.531809   .9974043    -0.53   0.594    -2.486686    1.423068
                            |
                      _cons |   51.34106   .4711164   108.98   0.000     50.41769    52.26443
---------------------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
habneigh2: Unstructured      |
               var(edcatd~2) |   .0029032   .0653538      2.00e-22    4.21e+16
                  var(_cons) |   .3470457    .604004      .0114536    10.51554
         cov(edcatd~2,_cons) |   .0317419   .3348349     -.6245224    .6880061
-----------------------------+------------------------------------------------
               var(Residual) |   67.42461   1.538731      64.47521    70.50894
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 1.48                  Prob > chi2 = 0.6878

. mixed WEMWBStotal i.agecat09 NHses09##empdi_indicators* if sex2==2 || habneigh2:empdi_indicators*, cov(unstr) diff

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -14126.064  
Iteration 1:   log likelihood =  -14118.78  
Iteration 2:   log likelihood = -14118.583  
Iteration 3:   log likelihood = -14118.565  
Iteration 4:   log likelihood = -14118.565  

Computing standard errors:
standard-error calculation has failed

Mixed-effects ML regression                     Number of obs     =      4,006
Group variable: habneigh2                       Number of groups  =        200

                                                Obs per group:
                                                              min =          3
                                                              avg =       20.0
                                                              max =         66

                                                Wald chi2(18)     =     146.43
Log likelihood = -14118.565                     Prob > chi2       =     0.0000

-------------------------------------------------------------------------------------------
              WEMWBStotal |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------------------+----------------------------------------------------------------
                 agecat09 |
                   47-51  |   .1233039   .4110411     0.30   0.764    -.6823218    .9289296
                   52-56  |   1.267891   .4112608     3.08   0.002      .461835    2.073948
                   57-61  |    1.76422   .4201375     4.20   0.000     .9407656    2.587674
                   62-70  |   2.526932   .4691344     5.39   0.000     1.607445    3.446418
                          |
                  NHses09 |
                      Q2  |  -.3767645   .4628824    -0.81   0.416    -1.283997    .5304683
                      Q3  |  -.7741833   .4700134    -1.65   0.100    -1.695393     .147026
                      Q4  |  -1.331949    .479685    -2.78   0.005    -2.272114   -.3917831
            Q5(most dis)  |  -1.902071   .5420812    -3.51   0.000     -2.96453   -.8396112
                          |
      1.empdi_indicators2 |  -.3078187   .7541393    -0.41   0.683    -1.785905    1.170267
      1.empdi_indicators3 |   1.033363   .8510056     1.21   0.225    -.6345769    2.701304
                          |
NHses09#empdi_indicators2 |
                    Q2#1  |    .855646   1.140097     0.75   0.453    -1.378903    3.090195
                    Q3#1  |  -2.015661   1.188872    -1.70   0.090    -4.345807    .3144843
                    Q4#1  |  -.6390923   1.154131    -0.55   0.580    -2.901148    1.622964
          Q5(most dis)#1  |  -4.360141   1.179882    -3.70   0.000    -6.672668   -2.047614
                          |
NHses09#empdi_indicators3 |
                    Q2#1  |  -.2129341   1.194731    -0.18   0.859    -2.554564    2.128696
                    Q3#1  |  -.7009247   1.204684    -0.58   0.561    -3.062062    1.660213
                    Q4#1  |  -.0068238   1.158589    -0.01   0.995    -2.277617     2.26397
          Q5(most dis)#1  |   -.517347   1.296921    -0.40   0.690    -3.059266    2.024572
                          |
                    _cons |   50.79966   .4085112   124.35   0.000     49.99899    51.60033
-------------------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
habneigh2: Unstructured      |
               var(empdi_~2) |   .9821291          .             .           .
               var(empdi_~3) |    .082077          .             .           .
                  var(_cons) |   .2997875          .             .           .
      cov(empdi_~2,empdi_~3) |   .2839194          .             .           .
         cov(empdi_~2,_cons) |   .5426139          .             .           .
         cov(empdi_~3,_cons) |   .1568619          .             .           .
-----------------------------+------------------------------------------------
               var(Residual) |   66.82862          .             .           .
------------------------------------------------------------------------------
LR test vs. linear model: chi2(6) = 4.35                  Prob > chi2 = 0.6290

Note: LR test is conservative and provided only for reference.
Warning: standard-error calculation failed