Using suest with Multi-level mixed effects models

[Alison Casey]

I am running some stratified multiple regression models where I am using the suest and the posestimation "test" command to see if the relationship between my DV and IV is significantly different across my two subsamples (first-gen students and continuing-gen students). Prior to multiple imputation, I had been using the suest command followed by the test command to accomplish this with the following code:

logit everenrollcollege bioguardian singleparent guardianonly twobioparent if firstgen==0
est sto cgen
logit everenrollcollege bioguardian singleparent guardianonly twobioparent if firstgen==1
est sto fgen
suest cgen fgen, vce(cluster schoolid) cformat(%9.3f)
test [cgen_everenrollcollege]singleparent=[fgen_everenrollcollege]singleparent

This gives me what I'm looking for with no issues. However, if I try to run multi-level mixed effects regression, store the results, and use suest to estimate both models simultaneously, I get an error message stating that meglm is not supported by suest r(322);. The below code produces the error message:

melogit everenrollcollege bioguardian singleparent guardianonly twobioparent if firstgen==0 || id:
est sto cgen
melogit everenrollcollege bioguardian singleparent guardianonly twobioparent if firstgen==1 || id:
est sto fgen
suest cgen fgen, vce(cluster schoolid) cformat(%9.3f)

I essentially just need to be able to test to see whether the differences in the coefficients for the "singleparent" variable between continuing-gen and first-gen students are significant. Any help on this is much appreciated!

 * Example generated by -dataex-. To install: ssc install dataex
clear
input float(firstgen everenrollcollege bioguardian singleparent guardianonly twobioparent)
1 1 0 0 0 1
1 1 0 0 0 1
1 1 0 1 0 0
1 1 0 0 0 1
1 0 0 1 0 0
0 1 0 0 0 1
1 1 0 0 0 1
1 1 0 0 0 1
0 1 1 0 0 0
1 1 1 0 0 0
0 1 0 0 0 1
0 1 0 1 0 0
1 0 0 0 0 1
0 1 0 0 0 1
1 1 0 0 0 1
1 0 0 1 0 0
1 1 0 0 0 1
1 1 0 0 0 1
0 1 0 0 1 0
1 1 0 0 0 1
1 1 0 0 0 1
1 1 0 0 0 1
1 0 1 0 0 0
1 0 0 1 0 0
1 1 0 0 0 1
0 0 1 0 0 0
1 1 0 0 0 1
1 1 0 1 0 0
1 1 0 0 0 1
1 0 1 0 0 0
1 1 0 1 0 0
1 1 0 1 0 0
1 0 0 1 0 0
1 1 0 0 0 1
1 0 1 0 0 0
1 1 1 0 0 0
0 1 0 0 0 1
1 1 0 1 0 0
1 1 0 1 0 0
end  
label values firstgen firstgencategory
label def firstgencategory 0 "Non-first-gen", modify  
label def firstgencategory 1 "First-gen", modify  
label values everenrollcollege everattendcategory  
label def everattendcategory 0 "0.No", modify  
label def everattendcategory 1 "1.Yes", modify

[Clyde Schechter]

You can't run -suest- after -me- commands. You can determine whether coefficients differ between the firstgen groups by doing an interaction model:

melogit everenrollcollege i.firstgen##(i.bioguardian i.singleparent i.guardianonly i.twobioparent) || id:

The difference in the coefficient of variable X between firstgen == 1 and firstgen == 0 will be found in the 1.firstgen#1.X row in the -melogit- output, along with standard error, confidence interval, and test statistics.

[Lanny Martin]

Given the -suest- limitation after -me-, here's a workaround using -gsem- (drawing on the fact that it allows you to stack equations).

First, you'll have to create two separate dependent variables corresponding to your groups, e.g.:

gen everenrollcollege_firstno = everenrollcollege if firstgen == 0
gen everenrollcollege_firstyes = everenrollcollege if firstgen == 1

Then estimate with -gsem-:

gsem (everenrollcollege_firstno <- bioguardian singleparent guardianonly twobioparent M1[id]) ///
(everenrollcollege_firstyes <- bioguardian singleparent guardianonly twobioparent M2[id]), ///
family(bernoulli) link(logit) vce(robust) cov(M1[id]*M2[id]@0)

And for the coefficient equality test:

test [everenrollcollege_firstno]singleparent = [everenrollcollege_firstyes]singleparent

[Jeff Pitblado (StataCorp)]

Also see help gsem group options.

[Joseph Coveney]

. . . here's a workaround using -gsem- . . .

I notice that you didn't constrain the two random effects to have the same variance. Doing so will give you the same results as the interaction model that Clyde shows above in #2 (you can run the code below to see what I mean).

In contrast, you use the vce(robust) option. That follows the same tack that suest uses for non-multilevel models, but it's not clear to me what that buys you in this context. Is it to accommodate possible heteroscedasticity in the random effects variances?

And if so does it? More generally, is there a consensus as to how to fulfill OP's request to conduct what seems essentially to be a Chow test when the models are hierarchical / multilevel?

version 18.0

clear *

// seedem
set seed 687922220

quietly set obs 250
generate int pid = _n
generate double pid_u = rnormal()

generate byte trt = mod(_n, 2)

quietly expand 10
generate double pre = runiform(-0.5, 0.5)

generate double xbu = trt - 0.5 + pid_u
generate byte out = rbinomial(1, invlogit(xbu))

// Interaction model
melogit out i.trt##c.pre || pid: , nolog
lincom _b[1.trt#c.pre]

// Equivalent -gsem- model
quietly generate byte out0 = out if !trt
quietly generate byte out1 = out if  trt

constraint define 1 _b[/:var(M1[pid])] = _b[/:var(M0[pid])]
gsem ///
    (out0 <- c.pre M0[pid]) ///
    (out1 <- c.pre M1[pid]), logit ///
        constraints(1) covariance(M0[pid]*M1[pid]@0) ///
            nocnsreport nodvheader nolog
lincom _b[out1:pre] - _b[out0:pre]

// Same, using -gsem-'s -group()- options
gsem (out <- c.pre M[pid]), logit ///
    group(trt) ginvariant(means covex) ///
        nocnsreport nodvheader nolog
lincom _b[out:1.trt#c.pre] - _b[out:0.trt#c.pre]

exit

[Lanny Martin]

On constraining the variance, I suppose one could let the data speak to that. Running your code, there's no difference in goodness of fit between the constrained and unconstrained model, but that may not be the case with the OP's data. On the VCE option, I was just mimicking -suest-, but right, not sure if it buys much. Here's a link to a pretty old paper that suggests that -vce(robust)- might be a good idea depending on sample size at level 2 and the size of the ICC. https://www.joophox.net/publist/sn04.pdf I don't know if it's a settled issue.

[Joseph Coveney]

Yes, gsem makes fitting separate variances easy. It can be done with melogit, too, but I think that it's not conventionally done, perhaps because of the extra effort and that the focus is nearly always on the fixed effects.

Reading that paper—thank you, by the way, I was not aware of it—kind of deflates my estimation of the benefit of robust standard errors. The only noticeable benefit there (linear models) was in the standard errors of the variance components when the residual error distribution is grossly mispecified, and again interest typically lies elsewhere, in fixed effects.