Multi-level logistic regression models and intraclass corellation coefficient

[Rune Lomholt]

Hello everybody. I have used the following syntax to calculate a multi-level multinominal logistic regression model on nested data. students are nested within schools.

Syntax for the model is:

gsem (fjob -> 2.stud, family(multinomial) link(logit)) (fjob -> 3.stud, family(multinomial) link(logit)) (fjob -> 4.stud, family(multinomial) link(logit)) (friends -> 2.stud, family(multinomial) link(logit)) (friends -> 3.stud, family(multinomial) link(logit)) (friends -> 4.stud, family(multinomial) link(logit)) (academic -> 2.stud, family(multinomial) link(logit)) (academic -> 3.stud, family(multinomial) link(logit)) (academic -> 4.stud, family(multinomial) link(logit)) (personal -> 2.stud, family(multinomial) link(logit)) (personal -> 3.stud, family(multinomial) link(logit)) (personal -> 4.stud, family(multinomial) link(logit)) (bildung -> 2.stud, family(multinomial) link(logit)) (bildung -> 3.stud, family(multinomial) link(logit)) (bildung -> 4.stud, family(multinomial) link(logit)) (parents -> 2.stud, family(multinomial) link(logit)) (parents -> 3.stud, family(multinomial) link(logit)) (parents -> 4.stud, family(multinomial) link(logit)) (chal -> 2.stud, family(multinomial) link(logit)) (chal -> 3.stud, family(multinomial) link(logit)) (chal -> 4.stud, family(multinomial) link(logit)) (M28[gymid] -> 2.stud, family(multinomial) link(logit)) (M28[gymid] -> 3.stud, family(multinomial) link(logit)) (M28[gymid] -> 4.stud, family(multinomial) link(logit)) (køn -> 2.stud, family(multinomial) link(logit)) (køn -> 3.stud, family(multinomial) link(logit)) (køn -> 4.stud, family(multinomial) link(logit)) (alder -> 2.stud, family(multinomial) link(logit)) (alder -> 3.stud, family(multinomial) link(logit)) (alder -> 4.stud, family(multinomial) link(logit)) (se -> 2.stud, family(multinomial) link(logit)) (se -> 3.stud, family(multinomial) link(logit)) (se -> 4.stud, family(multinomial) link(logit)) (si -> 2.stud, family(multinomial) link(logit)) (si -> 3.stud, family(multinomial) link(logit)) (si -> 4.stud, family(multinomial) link(logit)), covstruct(_lexogenous, diagonal) nocnsreport cformat(%9.2f) pformat(%5.3f) sformat(%8.2f) latent(M28 ) nocapslatent

The analysis shows that variation at the school level is a significant predictor of the outcome variable. How do i calculate the intraclass correlation coefficient to show how much variation that the school level can account for? or is it hidden somewhere in the output.

Hope you can help

[Clyde Schechter]

Wow, that is the most complicated way I can think of to code this model. And yes, buried in the output somewhere is the variance of M28[gymid], which you could then use to calculate the ICC, using the fact that the lowest level variance in a logistic model is always pi²/3. But it would be much easier to just redo the model in simpler terms:

xtset gymid
xtmlogit stud fjob friends academic personal bildung parents chal køn alder se, re

The output will be much more compact, and near the end of the regression output table you will get a variance component for gymid. Let's say its value is VC. Then the ICC is VC/(VC+(pi²/3).

[Joseph Coveney]

Just to build on Clyde's excellent advice, if you want to use -gsem-, then you'll probably want to constrain the factor loadings on the random effect to be one. That's done with the syntax shown in the user's manual for SEM (see Example 41g) for fitting this model with -gsem-, namely, as follows (in red).

gsem (stud <- fjob friends academic personal bildung parents chal køn alder se M28[gymid]@1), mlogit

If any of your predictors is categorical, keep in mind that you can use Stata's factor variable notation with -gsem-, just like with most other estimation commands (including, of course, -xtmlogit-).

Also, I'm not sure what you had in mind with the covstruct(_lexogenous, diagonal) option, but I think that it doesn't do anything in the context of your particular model.

Last, if you want to use -mlogit- as Clyde recommends, then you'll probably want to specify the -covariance()- option, as the default is to fit a separate random effect to each independent level of the outcome variable, which makes ICC computation problematic.

Using a toy example (only two predictors), I've shown an example of how you might compute the ICC below. Begin at the "Begin here" comment; the top part is just to create the fictitious dataset for illustration.

.ÿ
.ÿversionÿ17.0

.ÿ
.ÿclearÿ*

.ÿ
.ÿ//ÿseedem
.ÿsetÿseedÿ2079322411

.ÿ
.ÿquietlyÿsetÿobsÿ500

.ÿ
.ÿgenerateÿintÿpidÿ=ÿ_n

.ÿgenerateÿdoubleÿpid_uÿ=ÿrnormal(0,ÿ_piÿ/ÿsqrt(3))ÿ//ÿICCÿca.ÿ0.5

.ÿ
.ÿquietlyÿexpandÿ3

.ÿbysortÿpid:ÿgenerateÿbyteÿtimÿ=ÿ_n

.ÿ
.ÿgenerateÿdoubleÿpreÿ=ÿruniform(-0.5,ÿ0.5)

.ÿgenerateÿdoubleÿxbuÿ=ÿpreÿ+ÿ(timÿ-ÿ2)ÿ/ÿ2ÿ+ÿpid_u

.ÿ
.ÿgenerateÿdoubleÿdenÿ=ÿ1ÿ+ÿ3ÿ*ÿexp(xbu)

.ÿgenerateÿdoubleÿp1ÿ=ÿ1ÿ/ÿden

.ÿforvaluesÿiÿ=ÿ2/4ÿ{
ÿÿ2.ÿÿÿÿÿÿÿÿÿgenerateÿdoubleÿp`i'ÿ=ÿexp(xbu)ÿ/ÿden
ÿÿ3.ÿ}

.ÿ
.ÿgenerateÿdoubleÿranduÿ=ÿruniform()

.ÿgenerateÿbyteÿoutÿ=ÿ///
>ÿÿÿÿÿÿÿÿÿcond(randuÿ<ÿp1,ÿ1,ÿ///
>ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿcond(randuÿ<ÿp1ÿ+ÿp2,ÿ2,ÿ///
>ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿcond(randuÿ<ÿp1ÿ+ÿp2ÿ+ÿp3,ÿ3,ÿ4)))

.ÿ
.ÿ*
.ÿ*ÿBeginÿhere
.ÿ*
.ÿgsemÿ(outÿ<-ÿc.preÿi.timÿM[pid]@1),ÿmlogitÿ///
>ÿÿÿÿÿÿÿÿÿnocnsreportÿnodvheaderÿnolog

GeneralizedÿstructuralÿequationÿmodelÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿ=ÿ1,500
Logÿlikelihoodÿ=ÿ-1973.0416

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿCoefficientÿÿStd.ÿerr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
1.outÿÿÿÿÿÿÿÿ|ÿÿ(baseÿoutcome)
-------------+----------------------------------------------------------------
2.outÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿpreÿ|ÿÿÿ.5801009ÿÿÿ.3218026ÿÿÿÿÿ1.80ÿÿÿ0.071ÿÿÿÿ-.0506206ÿÿÿÿ1.210822
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿtimÿ|
ÿÿÿÿÿÿÿÿÿÿ2ÿÿ|ÿÿÿ.5675298ÿÿÿ.2027537ÿÿÿÿÿ2.80ÿÿÿ0.005ÿÿÿÿÿÿÿ.17014ÿÿÿÿ.9649197
ÿÿÿÿÿÿÿÿÿÿ3ÿÿ|ÿÿÿ1.125366ÿÿÿ.2058083ÿÿÿÿÿ5.47ÿÿÿ0.000ÿÿÿÿÿÿ.721989ÿÿÿÿ1.528743
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿM[pid]ÿ|ÿÿÿÿÿÿÿÿÿÿ1ÿÿ(constrained)
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿ_consÿ|ÿÿ-.4225354ÿÿÿ.1673218ÿÿÿÿ-2.53ÿÿÿ0.012ÿÿÿÿ-.7504801ÿÿÿ-.0945908
-------------+----------------------------------------------------------------
3.outÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿpreÿ|ÿÿÿ.6452116ÿÿÿ.3252807ÿÿÿÿÿ1.98ÿÿÿ0.047ÿÿÿÿÿ.0076732ÿÿÿÿÿ1.28275
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿtimÿ|
ÿÿÿÿÿÿÿÿÿÿ2ÿÿ|ÿÿÿ.6547623ÿÿÿ.2005069ÿÿÿÿÿ3.27ÿÿÿ0.001ÿÿÿÿÿÿ.261776ÿÿÿÿ1.047749
ÿÿÿÿÿÿÿÿÿÿ3ÿÿ|ÿÿÿ.8413271ÿÿÿ.2111619ÿÿÿÿÿ3.98ÿÿÿ0.000ÿÿÿÿÿ.4274573ÿÿÿÿ1.255197
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿM[pid]ÿ|ÿÿÿÿÿÿÿÿÿÿ1ÿÿ(constrained)
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿ_consÿ|ÿÿ-.4145052ÿÿÿ.1670863ÿÿÿÿ-2.48ÿÿÿ0.013ÿÿÿÿ-.7419884ÿÿÿÿ-.087022
-------------+----------------------------------------------------------------
4.outÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿpreÿ|ÿÿÿ.6151436ÿÿÿ.3230071ÿÿÿÿÿ1.90ÿÿÿ0.057ÿÿÿÿ-.0179387ÿÿÿÿ1.248226
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿtimÿ|
ÿÿÿÿÿÿÿÿÿÿ2ÿÿ|ÿÿÿ.7129254ÿÿÿ.2039018ÿÿÿÿÿ3.50ÿÿÿ0.000ÿÿÿÿÿ.3132852ÿÿÿÿ1.112566
ÿÿÿÿÿÿÿÿÿÿ3ÿÿ|ÿÿÿ1.195315ÿÿÿ.2084614ÿÿÿÿÿ5.73ÿÿÿ0.000ÿÿÿÿÿ.7867377ÿÿÿÿ1.603891
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿM[pid]ÿ|ÿÿÿÿÿÿÿÿÿÿ1ÿÿ(constrained)
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿ_consÿ|ÿÿ-.5148074ÿÿÿ.1700996ÿÿÿÿ-3.03ÿÿÿ0.002ÿÿÿÿ-.8481964ÿÿÿ-.1814183
-------------+----------------------------------------------------------------
ÿÿvar(M[pid])|ÿÿÿÿ3.43865ÿÿÿÿ.613552ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ2.423872ÿÿÿÿ4.878275
------------------------------------------------------------------------------

.ÿxtmlogitÿoutÿc.preÿi.tim,ÿi(pid)ÿreÿcovariance(shared)ÿnocnsreportÿnolog

Random-effectsÿmultinomialÿlogisticÿregressionÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿ=ÿÿ1,500
Groupÿvariable:ÿpidÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿgroupsÿ=ÿÿÿÿ500

Randomÿeffectsÿu_iÿ~ÿGaussianÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿObsÿperÿgroup:
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿminÿ=ÿÿÿÿÿÿ3
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿavgÿ=ÿÿÿÿ3.0
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿmaxÿ=ÿÿÿÿÿÿ3

Integrationÿmethod:ÿmvaghermiteÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿIntegrationÿpts.ÿ=ÿÿÿÿÿÿ7

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿWaldÿchi2(9)ÿÿÿÿÿ=ÿÿ46.11
Logÿlikelihoodÿ=ÿ-1973.0416ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>ÿchi2ÿÿÿÿÿÿ=ÿ0.0000

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿoutÿ|ÿCoefficientÿÿStd.ÿerr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
1ÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿ(baseÿoutcome)
-------------+----------------------------------------------------------------
2ÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿpreÿ|ÿÿÿ.5801009ÿÿÿ.3218026ÿÿÿÿÿ1.80ÿÿÿ0.071ÿÿÿÿ-.0506206ÿÿÿÿ1.210822
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿtimÿ|
ÿÿÿÿÿÿÿÿÿÿ2ÿÿ|ÿÿÿ.5675298ÿÿÿ.2027537ÿÿÿÿÿ2.80ÿÿÿ0.005ÿÿÿÿÿÿÿ.17014ÿÿÿÿ.9649197
ÿÿÿÿÿÿÿÿÿÿ3ÿÿ|ÿÿÿ1.125366ÿÿÿ.2058083ÿÿÿÿÿ5.47ÿÿÿ0.000ÿÿÿÿÿÿ.721989ÿÿÿÿ1.528743
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿ_consÿ|ÿÿ-.4225354ÿÿÿ.1673218ÿÿÿÿ-2.53ÿÿÿ0.012ÿÿÿÿ-.7504801ÿÿÿ-.0945908
-------------+----------------------------------------------------------------
3ÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿpreÿ|ÿÿÿ.6452116ÿÿÿ.3252807ÿÿÿÿÿ1.98ÿÿÿ0.047ÿÿÿÿÿ.0076732ÿÿÿÿÿ1.28275
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿtimÿ|
ÿÿÿÿÿÿÿÿÿÿ2ÿÿ|ÿÿÿ.6547623ÿÿÿ.2005069ÿÿÿÿÿ3.27ÿÿÿ0.001ÿÿÿÿÿÿ.261776ÿÿÿÿ1.047749
ÿÿÿÿÿÿÿÿÿÿ3ÿÿ|ÿÿÿ.8413271ÿÿÿ.2111619ÿÿÿÿÿ3.98ÿÿÿ0.000ÿÿÿÿÿ.4274573ÿÿÿÿ1.255197
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿ_consÿ|ÿÿ-.4145052ÿÿÿ.1670863ÿÿÿÿ-2.48ÿÿÿ0.013ÿÿÿÿ-.7419884ÿÿÿÿ-.087022
-------------+----------------------------------------------------------------
4ÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿpreÿ|ÿÿÿ.6151436ÿÿÿ.3230071ÿÿÿÿÿ1.90ÿÿÿ0.057ÿÿÿÿ-.0179387ÿÿÿÿ1.248226
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿtimÿ|
ÿÿÿÿÿÿÿÿÿÿ2ÿÿ|ÿÿÿ.7129254ÿÿÿ.2039018ÿÿÿÿÿ3.50ÿÿÿ0.000ÿÿÿÿÿ.3132852ÿÿÿÿ1.112566
ÿÿÿÿÿÿÿÿÿÿ3ÿÿ|ÿÿÿ1.195315ÿÿÿ.2084614ÿÿÿÿÿ5.73ÿÿÿ0.000ÿÿÿÿÿ.7867377ÿÿÿÿ1.603891
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿ_consÿ|ÿÿ-.5148074ÿÿÿ.1700996ÿÿÿÿ-3.03ÿÿÿ0.002ÿÿÿÿ-.8481964ÿÿÿ-.1814183
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿvar(u)|ÿÿÿÿ3.43865ÿÿÿÿ.613552ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ2.423872ÿÿÿÿ4.878275
------------------------------------------------------------------------------
LRÿtestÿvs.ÿmultinomialÿlogit:ÿchibar2(01)ÿ=ÿ142.86ÿÿÿProbÿ>=ÿchibar2ÿ=ÿ0.0000

.ÿ
.ÿ//ÿandÿhere
.ÿtempnameÿVC

.ÿscalarÿdefineÿ`VC'ÿ=ÿr(table)[1,ÿ`=colsof(r(table))']

.ÿdisplayÿinÿsmclÿasÿtextÿ"ICCÿ=ÿ"ÿ%04.2fÿ`VC'ÿ/ÿ(`VC'ÿ+ÿ_pi^2ÿ/ÿ3)
ICCÿ=ÿ0.51

.ÿ
.ÿexit

endÿofÿdo-file


.

[Rune Lomholt]

Hello both

Thank you so much for your reply. You don't know how helpful this is to me. I am not sure I understand Josephs reply though.

To Clyde:
The reason the code is so complicated is that i used the stata GSEM module to build the model and did build it using syntax.

I have a few questions regarding your approach to building the model:
1. Am I to understand that by declaring

gymid

as a panel variable, you then treat it automatically as level 2 using the xtmlogit syntax?
2. I have begun to wonder about the nested structure of my data. What if schools are nested within educational programs? So educational programs is level 3, and schools level 2 and respondents level 1? How would you model that?
3. Can you provide an example of syntax that could be used calculate the intraclass corellation coefficient based on the variance component for gymid?

[Clyde Schechter]

1. That's correct.

2. Well, if your outcome were binomial rather than multinomial, this would be rather simple, using the -melogit- command. Unfortunately, there is no corresponding simple multi-level model for multnomial logistic regression, and -xtmlogit- is restricted to 2 levels. So this would require you to go back to -gsem- It does not have to be complicated; it would look like this:

gsem (stud <- fjob friends academic personal bildung parents chal køn alder se, ///
    M1[program]@1 M2[program>gymid]@1, mlogit)

By the way, I don't know if your explanatory variables are discrete or continuous. But if any are discrete, those should be prefixed with i. so that Stata will treat them as such.

3. Joseph Coveney shows that in his response at #3.

[Rune Lomholt]

Thank you so muc CLyde.

One last question. Above you provided at formula to calculate the ICC. Is the ICC the same as the amount of variance in the outcome variable that can be explained by variation at the second level?

[Rune Lomholt]

I should at that taking your formula with a variance component of 11.46 i get a ICC = .78. That seems to me like a very high number

[Clyde Schechter]

Is the ICC the same as the amount of variance in the outcome variable that can be explained by variation at the second level?

No. Of the variance explained by the second level and first level, it is the proportion that is explained by the second. The amount at the second level would just be the 11.46 itself. Also, if you are doing a three-level model then you need to add the top-level variance component to both the numerator and denominator of that formula.

That seems to me like a very high number

Yes, but that is pretty typical in logistic models.

[Rune Lomholt]

No. Of the variance explained by the second level and first level, it is the proportion that is explained by the second. The amount at the second level would just be the 11.46 itself. Also, if you are doing a three-level model then you need to add the top-level variance component to both the numerator and denominator of that formula.


Yes, but that is pretty typical in logistic models.

Thank you for your reply Clyde. This has really been helpful to me and really made some aspects of the research I am doing much much better.
Do you have a citation for this problem you mention of inflated ICC in logistic models?

[Clyde Schechter]

No, I don't have a citation for that. And I wouldn't refer to it as "inflated." Remember that the outcome variable in a logistic model is just a 0/1 variable, so its variation is limited. Now, the variation being assessed by the variance components in the multilevel logistic model is not specifically that of the outcome variable, it's that of a latent continuous variable that is related to the outcome by a threshold process. But even so, it kind of makes sense.

All I'm relating to you is my experience with these models that the ICCs tend to be high. In particular, my experience with binomial outcomes, is that if you do both a logistic regression and a linear probability model with the same data, the ICC in the logistic regression is typically higher, often much higher. I have much less experience with multinomial outcomes, but the structure is somewhat similar, and I have observed the same phenomenon.

[Rune Lomholt]

No, I don't have a citation for that. And I wouldn't refer to it as "inflated." Remember that the outcome variable in a logistic model is just a 0/1 variable, so its variation is limited. Now, the variation being assessed by the variance components in the multilevel logistic model is not specifically that of the outcome variable, it's that of a latent continuous variable that is related to the outcome by a threshold process. But even so, it kind of makes sense.

All I'm relating to you is my experience with these models that the ICCs tend to be high. In particular, my experience with binomial outcomes, is that if you do both a logistic regression and a linear probability model with the same data, the ICC in the logistic regression is typically higher, often much higher. I have much less experience with multinomial outcomes, but the structure is somewhat similar, and I have observed the same phenomenon.

Okay CLyde. I undertand. Thank you su much.