How to conduct Rao/LM score test after gsem with random effects

[Sharad Kumar]

I want to conduct a score test to check whether the inclusion of a random intercept improves the model. I am running a duration model with frailty at the individual level (this is the random intercept). I currently fit the model as follows:

 gsem _t x1 x2 M1[x3], latent(M1) family(exponential, fail(failureVar) lt(_to)) nocapslatent)

where x3 is an identifier for individual.

gsem doesn't provide scoretests as a postestimation option (unlike SEM which allows estat scoretests). Does boottest or scoretest work for my situation? Can someone please let me know what syntax I should use for scoretest if I want to check whether the inclusion of the random intercept M1[x3] leads to better model fit? If scoretest doesn't allow this, do I have other options?

[Sharad Kumar]

David Roodman, I am tagging you in case you are around and able to help with this. I greatly appreciate any help.

[David Roodman]

I would answer this question with a question. To perform a score test for the hypothesis that the random effect exists, we need to tell gsem to run the model with the random effect included yet with its variance constrained to zero. And for scoretest to work with gsem, it must be possible to impose this constraint using gsem's constraints() option. Is there a way to do this? I haven't been able to find it.

Maybe it would just be better to do a likelihood ratio test?

--David

[Sharad Kumar]

Thanks, David. GSEM documentation implies that it is possible to constrain the variance within the model (see the "zero" option for structure in page 2 of this document: https://www.stata.com/manuals13/sems...vstructure.pdf). However, this doesn't work in practice and Stata returns the error "zero valued variances are not allowed" when I include covstructure(_LEx, zero) as an option in the gsem estimation. I also tried imposing this restriction through gsem's "var" option with the following syntax: var(_LEx@0), but I ran into the same error again. So I guess the simple answer to your question is No, I am unable to constrain the random effect to have a variance of 0 within the gsem model. I wonder why the documentation claims otherwise.

I did run likelihood ratio tests as you suggest. I ran into another potential issue with the LR tests (I say potential because I am actually not sure whether what I have is an issue). In short, the LR tests for my duration model yields a chi-sq of 0. This is because the RE model has a log likelihood that is almost identical to the no-RE model. This seems to be suggesting that there is no unobserved heterogeneity at the individual level (the level at which I introduce the RE). However, the fact that the two likelihoods are almost identical gives me pause. Such overlap seems unlikely to me, but I could be mistaken. This is the reason why I started exploring score tests. I was hoping to support the results of the LR tests with the score test. But alas, it looks like that is not an option.

If you are interested in this LR issue, you can view a more detailed post here:https://www.statalist.org/forums/for...se-for-concern. I provide the commands I ran and discuss the results in more detail.

Thanks so much for your help!

p.s. If you are reading the LR post linked above, please note that I use the command "mestreg" instead of "gsem". However, the difference is superficial since mestreg invokes gsem in the background, and I have made sure to specify gsem correctly. I decided to use gsem instead of mestreg because mestreg doesn't give you the option to let the duration parameters in the weibull and exponential baseline models vary by group. gsem allows this, and I might use this parametrization.

[Sharad Kumar]

Quick update. I found a way to constrain the variance using gsem's constraints() option. Unfortunately, I cannot constrain the variance to be exactly 0. I receive this error when I try: "initial values not feasible". I can set the variance to be as close to zero as 1e-200, but that's not exactly zero so I do not know if this is useful.

[David Roodman]

How do you constrain the variance? If you can define it in a constraint using the constraint command, then you can refer to it using the h0() option of scoretest. It could be useful to constrain it near zero but I would want to test for stability when using 1e-4, 1e-5, etc. gsem may well represent the variance in logs for model fitting, which could be why zero is not accepted.