[Liang Wang Jack]

You should also fit the real model in Stata.

You don't need to concoct a "specific mask" in order to make the comparisons that you want. Both margins and contrast will give you the comparisons that you desire after fitting the real model. (Run the do-file below to see that it is the case.)

Although I have doubts about its exogeneity—and I'm not sure why you sometimes have three categories and sometimes four in your snippets above—for mask, red (apparently high-probability of white matter hyper-intensity) in the corpus callosum means the same as red in the corona radiata, and in all of the brain tracts. Mask is crossed with brain tract. The model that you fit in Stata should include a term for the interaction of brain tract and mask, and not "specific mask".

version 15.1

clear *

set seed `=strreverse("1484913")'

quietly set obs 14
generate byte grp = _n > _N / 2
label variable grp Group

generate int pid = _n
label variable pid "Particpant ID"
generate double pid_u = rnormal(0, 2)
generate byte age = runiformint(18, 88)
generate byte sex = mod(_n, 2)

quietly expand 5
bysort pid: generate byte btc = _n
label define BrainTracts 1 "Cingulum bundle" 2 "Corona radiata" ///
3 "Corpus callosum" 4 "Internal capsule" 5 "Remaining white matter"
label values btc BrainTracts
label variable btc "Brain Tract"
generate double btc_u = rnormal(0, 1) // Assumption about the data-generating mechanism

quietly expand 4
bysort pid btc: generate byte msk = _n
label define Leukoaraiosis 1 "Red High-probability WMH" 2 "Green Low-probability WMH" 3 "Blue NAWM" 4 "God knows"
label values msk Leukoaraiosis
label variable msk Mask

generate double fai = normal(pid_u + btc_u + rnormal(0, 0.5))
label variable fai "Fractional Anisotropy"

*
* Begin here
*

// Recommend the following transformation in lieu of -vce(robust)-; could use arcsine square-root, too
generate double ifa = invnormal(fai)

anova ifa c.age sex grp / pid|grp btc grp#btc / btc|pid|grp msk btc#msk grp#msk grp#btc#msk

quietly mixed ifa c.age i.sex i.grp i.btc i.grp#i.btc i.msk i.btc#i.msk i.grp#i.msk i.grp#i.btc#i.msk ///
|| pid: || btc: , reml dfmethod(satterthwaite) nolrtest

margins btc#msk, dydx(grp) df(`=e(df_max)') mcompare(bonferroni)

contrast r.grp@btc#msk, df(`=e(df_max)') mcompare(bonferroni)

exit

Thank you so much for your suggestion!

Best regards,
Jack Liang