Foresplot

[laure rouch]

Hi all,

I would like to reproduce the attached foresplot (done in R) using Stata 15. I have the following linear mixed models (simplified version here) :

mixed score c.age##c.var1 i.sexep || id: age
mixed score c.age##c.var2 i.sexep || id: age
mixed score c.age##c.var3 i.sexep || id: age
etc. until var 6

For each of the 6 'model 1' models (the one with var 1 to the one with var 6), I would like to extract coefficients for var1 (simple effect) and coefficient for interaction var1 with cage as well as 95% CI and p values. I would like to reproduce the attached forest plot with for each of the 6 models (2 coefficients the one for the simple effect and the one for the interaction) with 95% CI and p value in each case. I would like to have the two coefficients for each model closer than the coefficients for the other models. I am a bit struggling with reproducing exactly this. I have used coefplot in the past but not exactly like this.

Any help would be really highly appreciated !
Thank you so much, L.

[Andrew Musau]

Can you re-upload the image? It's not visible to me. Please make sure it's in PNG format, as recommended in the FAQ advice.

[laure rouch]

Hi Andrews, here it is in .png. Thank you!

Best, Laure

[Andrew Musau]

Thanks. You can install coefplot from SSC and do the following:

webuse nlswork, clear
mixed ln_w c.grade##c.age || id:
est sto m1

mixed ln_w c.grade##c.age ttl_exp tenure c.tenure#c.tenure || id:
est sto m2

coefplot (m?, keep(grade)) (m?, keep(c.grade#c.age)), aseq swapname ///
ylab(1 "Model 1" 2 "Model 2", noticks) ///
leg(order(4 "grade" 2 "grade `=ustrunescape("\u00d7")' age"))

Res.:

[laure rouch]

Thank you, Andrew. In the following coefplot code: do you have any idea of how vertically displaying coefficients related to model A (coefficient var 1 and coefficient interaction var 1 with age) closer together than similar coefficients coming from model B (var 2 and interaction var 2 with age)? I would like a bigger space on y axis between coefficients coming from model A compared to coefficients coming from model B.

coefplot (modelA, keep(var1 c.age#var1)) (modelB, keep(var2 c.age#var2)), xline(0)

Thanks a lot!
Best, L.

[Andrew Musau]

If I understand correctly: Assuming the coefficients are named differently across models (but represent the same measurements), you can specify different offsets.

webuse nlswork, clear
mixed ln_w c.grade##c.age || id:
est sto m1

rename (grade age) (grade2 age2)
mixed ln_w c.grade##c.age2 ttl_exp tenure c.tenure#c.tenure || id:
est sto m2

local opts1 mcolor(red) ciopts(lc(red))
local opts2 mcolor(blue) ciopts(lc(blue))

coefplot (m?, keep(grade) offset(-.1) `opts1') (m?, keep(c.grade#c.age) offset(.1) `opts2') ///
(m?, keep(grade2) offset(-.4) `opts1') (m?, keep(c.grade2#c.age2) offset(.4) `opts2'), aseq swapname ///
ylab(1 "Model 1" 2 "Model 2", noticks) ///
leg(order(4 "grade" 2 "grade `=ustrunescape("\u00d7")' age"))

Res.:

[Andrew Musau]

Usually, showing the coefficients, confidence intervals (CIs), and p-values is somewhat redundant, assuming your readers can interpret coefficient plots. However, if you need to include them, here is a program written by Ben Jann that aligns them to the right. You may need to adjust the x-axis range and label gap to get the placement right.

capt program drop coefplot_mlbl2
*! version 1.0.0  10jun2021  Ben Jann
program coefplot_mlbl2, sclass
    _parse comma plots 0 : 0
    syntax [, MLabel(passthru) * ]
    if `"`mlabel'"'=="" local mlabel mlabel(string(@b))
    preserve
    qui coefplot `plots', `options' `mlabel' generate replace nodraw
    sreturn clear
    tempvar touse
    qui gen byte `touse' = 0
    local nplots = r(n_plots)
    forv i = 1/`nplots' {
        qui replace `touse' = __plot==`i' & __at<.
        mata: st_global("s(mlbl`i')", ///
            invtokens((strofreal(st_data(.,"__at","`touse'")) :+ " " :+ ///
            "`" :+ `"""' :+ st_sdata(.,"__mlbl","`touse'") :+ `"""' :+ "'")'))
    }
    sreturn local plots `"`plots'"'
    sreturn local options `"`options'"'
end


webuse nlswork, clear
mixed ln_w c.grade##c.age || id:
est sto m1

mixed ln_w c.grade##c.age ttl_exp tenure c.tenure#c.tenure || id:
est sto m2

coefplot_mlbl2 m?, keep(grade c.grade#c.age) xline(0) mlabel(string(@b,"%9.3f") + " ["+ string(@ll,"%9.3f") + ", " + string(@ul,"%9.3f") + "]; " + "p = "+string(@pval,"%9.3f")) 

coefplot (m?, keep(grade)) (m?, keep(c.grade#c.age)), aseq swapname ///
ylab(1 "Model 1" 2 "Model 2", noticks)  xsc(r(. .085)) ///
leg(order(4 "grade" 2 "grade `=ustrunescape("\u00d7")' age") pos(6) row(1) ) ///
ymlabel(`s(mlbl1)', angle(0) notick labgap(1) labsize(small) axis(2) add custom labcolor(navy)) ///
    ymlabel(`s(mlbl2)', angle(0) notick labgap(1)