Calculating an ICC value after multilevel multinomial logistic regression (gsem command)

[DY Kim]

Hi,

I am conducting the multilevel multinomial logistic regression. The multilevel multinomial logistic regression does not have a error variance term in the first level of the model. I just wonder how I can calculate an ICC value following executing the code below. Please let me know if you need additional information. Thank you.

. gsem (i.outcome <- ib(first).offrace offgender i.command c.ecodisadvantage M1[geoid]@1), mlogit

var(M1[geoid]) = > coefficient = .1563367; standard error = .0424555 ; 95% CI = .0918131 and .2662055

[Joseph Coveney]

I just wonder how I can calculate an ICC value following executing the code below.

I think that it might be a little bit tricky (see below) and maybe that's why Stata doesn't allow, say, an estat icc postestimation command with these models.

.ÿ
.ÿversionÿ18.0

.ÿ
.ÿclearÿ*

.ÿ
.ÿquietlyÿwebuseÿestatus

.ÿ
.ÿgsemÿ(ib3.estatusÿ<-ÿi.hhchildÿc.(ageÿhhincome)ÿi.hhsignoÿi.bwinnerÿM[id]@1,ÿ///
>ÿÿÿÿÿÿÿÿÿmlogit),ÿnocnsreportÿnodvheaderÿnolog

GeneralizedÿstructuralÿequationÿmodelÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿ=ÿ4,761
Logÿlikelihoodÿ=ÿ-4458.4888

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿCoefficientÿÿStd.ÿerr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
1.estatusÿÿÿÿ|
ÿÿÿÿÿhhchildÿ|
ÿÿÿÿÿÿÿÿYesÿÿ|ÿÿÿ.4769288ÿÿÿÿ.097508ÿÿÿÿÿ4.89ÿÿÿ0.000ÿÿÿÿÿ.2858167ÿÿÿÿ.6680409
ÿÿÿÿÿÿÿÿÿageÿ|ÿÿ-.0041234ÿÿÿ.0067956ÿÿÿÿ-0.61ÿÿÿ0.544ÿÿÿÿ-.0174424ÿÿÿÿ.0091957
ÿÿÿÿhhincomeÿ|ÿÿ-.0064771ÿÿÿ.0019207ÿÿÿÿ-3.37ÿÿÿ0.001ÿÿÿÿ-.0102416ÿÿÿ-.0027125
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿhhsignoÿ|
ÿÿÿÿÿÿÿÿYesÿÿ|ÿÿÿÿ.489292ÿÿÿ.0952363ÿÿÿÿÿ5.14ÿÿÿ0.000ÿÿÿÿÿ.3026324ÿÿÿÿ.6759517
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿbwinnerÿ|
ÿÿÿÿÿÿÿÿYesÿÿ|ÿÿ-.4839599ÿÿÿ.0733424ÿÿÿÿ-6.60ÿÿÿ0.000ÿÿÿÿ-.6277084ÿÿÿ-.3402115
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿM[id]ÿ|ÿÿÿÿÿÿÿÿÿÿ1ÿÿ(constrained)
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿ_consÿ|ÿÿ-.3122929ÿÿÿ.2895539ÿÿÿÿ-1.08ÿÿÿ0.281ÿÿÿÿÿ-.879808ÿÿÿÿ.2552223
-------------+----------------------------------------------------------------
2.estatusÿÿÿÿ|
ÿÿÿÿÿhhchildÿ|
ÿÿÿÿÿÿÿÿYesÿÿ|ÿÿÿ.0734195ÿÿÿ.1192622ÿÿÿÿÿ0.62ÿÿÿ0.538ÿÿÿÿ-.1603301ÿÿÿÿ.3071691
ÿÿÿÿÿÿÿÿÿageÿ|ÿÿÿ.0042289ÿÿÿ.0081244ÿÿÿÿÿ0.52ÿÿÿ0.603ÿÿÿÿ-.0116947ÿÿÿÿ.0201525
ÿÿÿÿhhincomeÿ|ÿÿ-.0306996ÿÿÿ.0025545ÿÿÿ-12.02ÿÿÿ0.000ÿÿÿÿ-.0357064ÿÿÿ-.0256928
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿhhsignoÿ|
ÿÿÿÿÿÿÿÿYesÿÿ|ÿÿÿ.1745492ÿÿÿ.1187011ÿÿÿÿÿ1.47ÿÿÿ0.141ÿÿÿÿ-.0581007ÿÿÿÿÿ.407199
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿbwinnerÿ|
ÿÿÿÿÿÿÿÿYesÿÿ|ÿÿ-.2648114ÿÿÿ.0946773ÿÿÿÿ-2.80ÿÿÿ0.005ÿÿÿÿ-.4503755ÿÿÿ-.0792472
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿM[id]ÿ|ÿÿÿÿÿÿÿÿÿÿ1ÿÿ(constrained)
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿ_consÿ|ÿÿÿ.0083575ÿÿÿ.3478268ÿÿÿÿÿ0.02ÿÿÿ0.981ÿÿÿÿ-.6733705ÿÿÿÿ.6900856
-------------+----------------------------------------------------------------
3.estatusÿÿÿÿ|ÿÿ(baseÿoutcome)
-------------+----------------------------------------------------------------
ÿÿÿvar(M[id])|ÿÿÿ1.002235ÿÿÿ.1143039ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.8014762ÿÿÿÿÿ1.25328
------------------------------------------------------------------------------

.ÿdisplayÿinÿsmclÿasÿtextÿ"ICCÿ=ÿ"ÿasÿresultÿ%04.2fÿ_b[/:var(M[id])]ÿ/ÿ///
>ÿÿÿÿÿÿÿÿÿ(_b[/:var(M[id])]ÿ+ÿ_pi^2ÿ/ÿ3)
ICCÿ=ÿ0.23

.ÿ
.ÿprogramÿdefineÿbasEm
ÿÿ1.ÿÿÿÿÿÿÿÿÿversionÿ18.0
ÿÿ2.ÿÿÿÿÿÿÿÿÿsyntaxÿanything(name=baseoutcome)
ÿÿ3.ÿ
.ÿÿÿÿÿÿÿÿÿquietlyÿxtmlogitÿestatusÿi.hhchildÿageÿhhincomeÿi.hhsignoÿi.bwinner,ÿi(id)ÿ///
>ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿbaseoutcome(`baseoutcome')ÿcovariance(shared)ÿnolog
ÿÿ4.ÿÿÿÿÿÿÿÿÿdisplayÿinÿsmclÿasÿtextÿ_newline(1)ÿ"BaseÿOutcomeÿ=ÿ`baseoutcome'"
ÿÿ5.ÿÿÿÿÿÿÿÿÿdisplayÿ"ICCÿ=ÿ"ÿasÿresultÿ%04.2fÿ_b[/:var(u)]ÿ/ÿ///
>ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ(_b[/:var(u)]ÿ+ÿ_pi^2ÿ/ÿ3)
ÿÿ6.ÿend

.ÿ
.ÿquietlyÿlevelsofÿestatus,ÿlocal(baseoutcomes)

.ÿ
.ÿforeachÿbaseoutcomeÿofÿlocalÿbaseoutcomesÿ{
ÿÿ2.ÿÿÿÿÿÿÿÿÿbasEmÿ`baseoutcome'
ÿÿ3.ÿ}

BaseÿOutcomeÿ=ÿ1
ICCÿ=ÿ0.20

BaseÿOutcomeÿ=ÿ2
ICCÿ=ÿ0.16

BaseÿOutcomeÿ=ÿ3
ICCÿ=ÿ0.23

.ÿ
.ÿexit

endÿofÿdo-file


.

[Leonardo Guizzetti]

I agree with Joseph's caution about being careful with what and how you interpret the ICC in a multinomial logistic regression. It's also worth noting that there isn't an error variance in the logit/mlogit models because there isn't a modeled error term analogous to the gaussian (normal) models. The variance being used is the one fitted to the between-subject variability from fitting a distribution to the random intercept, and scaling this result based on the log-odds distibution (the factor is pi^2 / 3).

I do think you could look at it as one kind of measure of predictive accuracy, which would extend from the logistic regression case. Then again, this measure for predictive accuracy is also fairly esoteric compared to more conventional choices. For example, if I collapse the multinomial outcome into a collection of binary variables (coded as one level vs the rest), then you would get 3 distinct ICC values that correspond exactly to those shown by Joseph. In this sense, it would relate to predictive accuracy for one outcome only vs all of the rest.

// level 3 vs (1 or 2)
gsem (estatus3 <- i.hhchild c.(age hhincome) i.hhsigno i.bwinner M[id]@1, logit), nocnsreport nodvheader nolog
display in smcl as text "ICC = " as result %04.2f _b[/:var(M[id])] / (_b[/:var(M[id])] + _pi^2 / 3)

// level 2 vs (1 or 3)
gsem (estatus2 <- i.hhchild c.(age hhincome) i.hhsigno i.bwinner M[id]@1, logit), nocnsreport nodvheader nolog
display in smcl as text "ICC = " as result %04.2f _b[/:var(M[id])] / (_b[/:var(M[id])] + _pi^2 / 3)

// level 1 vs (2 or 3)
gsem (estatus1 <- i.hhchild c.(age hhincome) i.hhsigno i.bwinner M[id]@1, logit), nocnsreport nodvheader nolog
display in smcl as text "ICC = " as result %04.2f _b[/:var(M[id])] / (_b[/:var(M[id])] + _pi^2 / 3)

Relevant output:

. // level 3 vs (1 or 2)
ICC = 0.23

. // level 2 vs (1 or 3)
. display in smcl as text "ICC = " as result %04.2f _b[/:var(M[id])] / (_b[/:var(M[id])] + _pi^2 / 3)
ICC = 0.16

. // level 1 vs (2 or 3)
. display in smcl as text "ICC = " as result %04.2f _b[/:var(M[id])] / (_b[/:var(M[id])] + _pi^2 / 3)
ICC = 0.20

[DY Kim]

Joseph and Leonardo,

You suggestions worked very well. Thank you.

Can I get a p-value for var(M[id])?
Can I get information on the within group variance and corresponding p-value?

[Leonardo Guizzetti]

Can I get a p-value for var(M[id])?

You could, but I don't see what good it will be. Even the Stata developers thought to not show this information.

Can I get information on the within group variance and corresponding p-value?

There's isn't a real within-group variance with these models, only between-group.

[Joseph Coveney]

xtmlogit does give a likelihood-ratio test of the multilevel / hierarchical model compared to that for mlogit, which is "a p-value for var(M[id])".

But in line with what Leonardo is implying, the p-value will vary according to which level (category) of the outcome variable is omitted, just as the ICC does.

version 18.0

clear *

quietly webuse estatus

quietly levelsof estatus, local(baseoutcomes)
foreach baseoutcome of local baseoutcomes {
    xtmlogit estatus i.hhchild c.(age hhincome) i.hhsigno i.bwinner, i(id) covariance(shared) ///
        baseoutcome(`baseoutcome') nolog
}

exit

The within-group variance is constrained for identifiability.

[DY Kim]

The random intercept model: gsem (i.outcome <- ib(first).offrace offgender i.command c.ecodisadvantage M1[geoid]@1), mlogit

The effect of offrace is assumed to vary across geoid. I wonder how to add a randomly varying slope to the above random intercept model.

If offrace interacts with the level 2 variable, c.ecodisadvantage, in affecting the outcome, how I can specify the Stata code. Any advice would be appreciated.