How to plot multilevel coefficients from mixed ? Can you use coefplot? If not a loop? macro?

[Brennan RBratton]

Hello Statalist & happy holidays!

I need to plot the curves from my mixed analysis. But I'm not quite sure what is the best method. Should I do this manually? Is there a command I can use? A loop? A macro?

The exposure is categorical so I have 3 curves to plot. Finding a way to do this automatically would be great.

Thank you for your insights.

[Joseph Coveney]

I need to plot the curves from my mixed analysis. But I'm not quite sure what is the best method. Should I do this manually? Is there a command I can use?

It's not clear whether you want to plot the coefficients, themselves (as your post's title indicates), or the model's marginal predictions (you've got interaction terms in your model, and so plotting the main effects coefficients might not be very informative).

If the latter, you'd use marginsplot in conjunction with margins as a postestimation command. You can find more information by tracing forward through the hyperlinks, starting with the help file for mixed, scrolling down until you come to the hyperlink for mixed postestimation and clicking on it.

[Brennan RBratton]

Joseph, thank you for your suggestion. I am hoping to create a figure like the one below that is from the above results. I used xtgraph to create this figure just to demonstrate what I am hoping to create with the coefficients for each type of exposure. I am guessing that I have to create 3 new variables with the following equations (one for each type of exposure) then plot the variables over time. But I was hoping that there a simple commands or code to do so

Not sure if this is necessary or needed but here is the model:
mixed outcome exposure_cat##c.Time exposure_cat##c.T_5 exposure_cat##c.T_10 exposure_cat##c.T_19 covariate1 covariate2 covariate3, || group: Time, residuals(ar1, t(T)) reml

Exposure C (reference group)
line for Exposure C = intercept + Exposure=C + time=1990-2014 + (Exposure=C*time=1990-2014) + time=1990-1995 + (Exposure=C*time=1990-1995) + time=1996-2000 + (Exposure=C*time=1996-2000) + time=2001-2008 + (Exposure=C*time=2001-2008) + time=2009-2014 + (Exposure=C*time=2009-2014) +Covariate1+ Covariate2 + Covariate3

Exposure A
line for Exposure A = intercept + Exposure=A + time=1990-2014 + (Exposure=A*time=1990-2014) + time=1990-1995 + (Exposure=A*time=1990-1995) + time=1996-2000 + (Exposure=A*time=1996-2000) + time=2001-2008 + (Exposure=A*time=2001-2008) + time=2009-2014 + (Exposure=A*time=2009-2014) +Covariate1+ Covariate2 + Covariate3

Exposure B
line for Exposure B = intercept + Exposure=B + time=1990-2014 + (Exposure=B *time=1990-2014) + time=1990-1995 + (Exposure=B *time=1990-1995) + time=1996-2000 + (Exposure=B *time=1996-2000) + time=2001-2008 + (Exposure=B *time=2001-2008) + time=2009-2014 + (Exposure=B *time=2009-2014) +Covariate1+ Covariate2 + Covariate3


Not sure if this is necessary or needed but here is the Stata code for the model:
mixed outcome exposure_cat##c.Time exposure_cat##c.T_5 exposure_cat##c.T_10 exposure_cat##c.T_19 covariate1 covariate2 covariate3, || group: Time, residuals(ar1, t(T)) reml

Stata code for the below graph:
xtgraph outcome, group(Exposure)

[Brennan RBratton]

My apologies for any confusion with my models or my request for assistance.

The model
includes time continuously (T) and three other time variables that represent the time periods of using 4 time1996-2014, 2001-2014, 2010-2014. Below is the code used to generate these variables. Using these time periods within this model help differentiate how the exposure impacts the outcome variable during each of these periods instead of only looking at time continuously.

T= Year-1990

These are the four time periods:
gen T = year-1990
label variable year " Year 0-24"

gen T_5 = (T-5)*(T>=5)
label variable T_5 "Time period 1: 1996-2014"

gen T_10 = (T-10)*(T>=10)
label variable T_10 "Time period 2: 2001-2014"

gen T_19 = (T-19)*(T>=19)
label variable T_19 "Time period 3: 2010-2014"

[Clyde Schechter]

OK, so you have generated a homebrew marginal spline for time. These are difficult to work with when you need -margins-, but there is no genuinely easier substitute. So we'll go with what you've done. The problem here is that -margins- does not know about the relationships among T, T_5, T_10, and T_19, so if you just naively use the -margins- command the way you ordinarily would, you will get seriously incorrect results. So you have to build those relationships into your -at()- options directly. To get started:

mixed outcome exposure_cat##c.(Time T_5 T_10 T_19) covariate1 covariate2 covariate3 ///
    || group: Time, residuals(ar1, t(T)) reml // THIS IS THE SAME AS ORIGINAL, JUST MORE COMPACTLY WRITTEN
local ats
forvalues y = 0/24 {
    foreach x in 5 10 19 {
        summ T_`x' if T == `y', meanonly
        local T_`x' = r(mean)
    }
    local ats `ats' at(T = `y' T_5 = `T_5' T_10 = `T-10' T_19 = `T-19')
}
margins exposure_cat, `ats' saving(margins_output, replace)
use margins_output, clear
gen year = 1990 + T

Now, as I do not have example data to work with, I have not tested this code and it may contain some errors. At this point you have in memory the data you need to generate the plot you want. The variable names are some special names created by -margins-, and without looking at an example, I don't know exactly what they will be. If you browse that data set it will be clear which variables are which, and if you like, you can rename them to whatever you find more convenient. From that point you want to use -graph two way line- and overlay -twoway rcap- with the appropriate variables, and a -by()- option specifying whatever the exposure category variable is now called, to get the graph you want.

If you need more specific advice, provide example data when you post back. Be sure to use -dataex- for that. If you are running version 16 or a fully updated version 15.1 or 14.2, -dataex- is already part of your official Stata installation. If not, run -ssc install dataex- to get it. Either way, run -help dataex- to read the simple instructions for using it. -dataex- will save you time; it is easier and quicker than typing out tables. It includes complete information about aspects of the data that are often critical to answering your question but cannot be seen from tabular displays or screenshots. It also makes it possible for those who want to help you to create a faithful representation of your example to try out their code, which in turn makes it more likely that their answer will actually work in your data.

[Brennan RBratton]

Dear Clyde,
Thank you so much for your suggestion. Below please find the output from dataex.

Also, instead of the "homebrew marginal spline" should I have created this using mkspline instead? Or another technique?

Thank you again for your insights.

* Example generated by -dataex-. To install: ssc install dataex
clear
input int year float(Time Tsq Tcube T_5 T_10 T_19 T_2008 outcome exposure_cat group)
1990  0   0     0  0  0 0 0  2.721215 0 101
1991  1   1     1  0  0 0 0  2.786076 0 101
1992  2   4     8  0  0 0 0  3.033087 0 101
1993  3   9    27  0  0 0 0 3.2524495 0 101
1994  4  16    64  0  0 0 0 3.2674875 0 101
1995  5  25   125  0  0 0 0  3.544667 0 101
1996  6  36   216  1  0 0 0  3.510141 0 101
1997  7  49   343  2  0 0 0  3.434771 0 101
1998  8  64   512  3  0 0 0 3.4271665 0 101
1999  9  81   729  4  0 0 0  3.419504 0 101
2000 10 100  1000  5  0 0 0  3.353254 0 101
2001 11 121  1331  6  1 0 0 3.4522934 0 101
2002 12 144  1728  7  2 0 0  3.726234 0 101
2003 13 169  2197  8  3 0 0  3.811199 0 101
2004 14 196  2744  9  4 0 0  3.812931 0 101
2005 15 225  3375 10  5 0 0  3.838551 0 101
2006 16 256  4096 11  6 0 0  3.963615 0 101
2007 17 289  4913 12  7 0 0  4.132956 0 101
2008 18 324  5832 13  8 0 1  4.160219 0 101
2009 19 361  6859 14  9 0 1 4.1819854 0 101
2010 20 400  8000 15 10 1 1 4.0172124 0 101
2011 21 441  9261 16 11 2 1 4.2183957 0 101
2012 22 484 10648 17 12 3 1 4.2067947 0 101
2013 23 529 12167 18 13 4 1  4.160219 0 101
2014 24 576 13824 19 14 5 1   4.17237 0 101
1990  0   0     0  0  0 0 0 2.8405066 0 102
1991  1   1     1  0  0 0 0  2.850923 0 102
1992  2   4     8  0  0 0 0 3.0349705 0 102
1993  3   9    27  0  0 0 0  3.190352 0 102
1994  4  16    64  0  0 0 0  3.194035 0 102
1995  5  25   125  0  0 0 0  3.440436 0 102
1996  6  36   216  1  0 0 0  3.408316 0 102
1997  7  49   343  2  0 0 0  3.402364 0 102
1998  8  64   512  3  0 0 0  3.402364 0 102
1999  9  81   729  4  0 0 0 3.4259636 0 102
2000 10 100  1000  5  0 0 0  3.393368 0 102
2001 11 121  1331  6  1 0 0   3.45753 0 102
2002 12 144  1728  7  2 0 0 3.6403224 0 102
2003 13 169  2197  8  3 0 0  3.751548 0 102
2004 14 196  2744  9  4 0 0  3.757844 0 102
2005 15 225  3375 10  5 0 0  3.820681 0 102
2006 16 256  4096 11  6 0 0 3.9495356 0 102
2007 17 289  4913 12  7 0 0  4.085008 0 102
2008 18 324  5832 13  8 0 1 4.1478486 0 102
2009 19 361  6859 14  9 0 1  4.137884 0 102
2010 20 400  8000 15 10 1 1  3.951261 0 102
2011 21 441  9261 16 11 2 1 4.1478486 0 102
2012 22 484 10648 17 12 3 1 4.1520886 0 102
2013 23 529 12167 18 13 4 1  4.122022 0 102
2014 24 576 13824 19 14 5 1  4.116191 0 102
1990  0   0     0  0  0 0 0 3.2943344 1 103
1991  1   1     1  0  0 0 0 3.3155365 1 103
1992  2   4     8  0  0 0 0 3.5344756 1 103
1993  3   9    27  0  0 0 0 3.7463195 1 103
1994  4  16    64  0  0 0 0 3.8037305 1 103
1995  5  25   125  0  0 0 0 4.1457057 1 103
1996  6  36   216  1  0 0 0  4.073004 1 103
1997  7  49   343  2  0 0 0  4.046424 1 103
1998  8  64   512  3  0 0 0 3.9767094 1 103
1999  9  81   729  4  0 0 0  3.991047 1 103
2000 10 100  1000  5  0 0 0 3.9056516 1 103
2001 11 121  1331  6  1 0 0 4.0328646 1 103
2002 12 144  1728  7  2 0 0  4.271586 1 103
2003 13 169  2197  8  3 0 0 4.3809705 1 103
2004 14 196  2744  9  4 0 0 4.3761396 1 103
2005 15 225  3375 10  5 0 0  4.383377 1 103
2006 16 256  4096 11  6 0 0 4.5117874 1 103
2007 17 289  4913 12  7 0 0  4.636826 1 103
2008 18 324  5832 13  8 0 1 4.6824126 1 103
2009 19 361  6859 14  9 0 1 4.6859717 1 103
2010 20 400  8000 15 10 1 1 4.4817424 1 103
2011 21 441  9261 16 11 2 1  4.658967 1 103
2012 22 484 10648 17 12 3 1 4.6405506 1 103
2013 23 529 12167 18 13 4 1  4.596859 1 103
2014 24 576 13824 19 14 5 1 4.5692983 1 103
1990  0   0     0  0  0 0 0 3.5519204 1 104
1991  1   1     1  0  0 0 0  3.587767 1 104
1992  2   4     8  0  0 0 0 3.7967255 1 104
1993  3   9    27  0  0 0 0 3.9979155 1 104
1994  4  16    64  0  0 0 0 3.9500594 1 104
1995  5  25   125  0  0 0 0 4.2450676 1 104
1996  6  36   216  1  0 0 0  4.165329 1 104
1997  7  49   343  2  0 0 0  4.070019 1 104
1998  8  64   512  3  0 0 0 4.0524745 1 104
1999  9  81   729  4  0 0 0 4.0524745 1 104
2000 10 100  1000  5  0 0 0 3.9979155 1 104
2001 11 121  1331  6  1 0 0 4.0999994 1 104
2002 12 144  1728  7  2 0 0 4.3555245 1 104
2003 13 169  2197  8  3 0 0 4.4898725 1 104
2004 14 196  2744  9  4 0 0 4.5152617 1 104
2005 15 225  3375 10  5 0 0 4.5400224 1 104
2006 16 256  4096 11  6 0 0 4.6887236 1 104
2007 17 289  4913 12  7 0 0  4.809884 1 104
2008 18 324  5832 13  8 0 1   4.86635 1 104
2009 19 361  6859 14  9 0 1   4.88769 1 104
2010 20 400  8000 15 10 1 1 4.6793118 1 104
2011 21 441  9261 16 11 2 1  4.864388 1 104
2012 22 484 10648 17 12 3 1   4.86045 1 104
2013 23 529 12167 18 13 4 1  4.811956 1 104
2014 24 576 13824 19 14 5 1  4.778269 1 104
end
label values T_2008 T_2008lab
label def T_2008lab 0 "pre-economic crisis", modify
label def T_2008lab 1 "post-economic crisis", modify
label values exposure_cat exp_lab
label def exp_lab 0 "ref", modify
label def exp_lab 1 "A", modify

[Clyde Schechter]

There were some typos in the -local ats …- command in #10. Those are fixed. And I follow with the massaging of the -margins- output and a start on the graph. The graph command gives you the plots you want. You can tinker with the labeling and other options affecting the appearance of the graph.

* Example generated by -dataex-. To install: ssc install dataex
clear
input int year float(Time Tsq Tcube T_5 T_10 T_19 T_2008 outcome exposure_cat group)
1990  0   0     0  0  0 0 0  2.721215 0 101
1991  1   1     1  0  0 0 0  2.786076 0 101
1992  2   4     8  0  0 0 0  3.033087 0 101
1993  3   9    27  0  0 0 0 3.2524495 0 101
1994  4  16    64  0  0 0 0 3.2674875 0 101
1995  5  25   125  0  0 0 0  3.544667 0 101
1996  6  36   216  1  0 0 0  3.510141 0 101
1997  7  49   343  2  0 0 0  3.434771 0 101
1998  8  64   512  3  0 0 0 3.4271665 0 101
1999  9  81   729  4  0 0 0  3.419504 0 101
2000 10 100  1000  5  0 0 0  3.353254 0 101
2001 11 121  1331  6  1 0 0 3.4522934 0 101
2002 12 144  1728  7  2 0 0  3.726234 0 101
2003 13 169  2197  8  3 0 0  3.811199 0 101
2004 14 196  2744  9  4 0 0  3.812931 0 101
2005 15 225  3375 10  5 0 0  3.838551 0 101
2006 16 256  4096 11  6 0 0  3.963615 0 101
2007 17 289  4913 12  7 0 0  4.132956 0 101
2008 18 324  5832 13  8 0 1  4.160219 0 101
2009 19 361  6859 14  9 0 1 4.1819854 0 101
2010 20 400  8000 15 10 1 1 4.0172124 0 101
2011 21 441  9261 16 11 2 1 4.2183957 0 101
2012 22 484 10648 17 12 3 1 4.2067947 0 101
2013 23 529 12167 18 13 4 1  4.160219 0 101
2014 24 576 13824 19 14 5 1   4.17237 0 101
1990  0   0     0  0  0 0 0 2.8405066 0 102
1991  1   1     1  0  0 0 0  2.850923 0 102
1992  2   4     8  0  0 0 0 3.0349705 0 102
1993  3   9    27  0  0 0 0  3.190352 0 102
1994  4  16    64  0  0 0 0  3.194035 0 102
1995  5  25   125  0  0 0 0  3.440436 0 102
1996  6  36   216  1  0 0 0  3.408316 0 102
1997  7  49   343  2  0 0 0  3.402364 0 102
1998  8  64   512  3  0 0 0  3.402364 0 102
1999  9  81   729  4  0 0 0 3.4259636 0 102
2000 10 100  1000  5  0 0 0  3.393368 0 102
2001 11 121  1331  6  1 0 0   3.45753 0 102
2002 12 144  1728  7  2 0 0 3.6403224 0 102
2003 13 169  2197  8  3 0 0  3.751548 0 102
2004 14 196  2744  9  4 0 0  3.757844 0 102
2005 15 225  3375 10  5 0 0  3.820681 0 102
2006 16 256  4096 11  6 0 0 3.9495356 0 102
2007 17 289  4913 12  7 0 0  4.085008 0 102
2008 18 324  5832 13  8 0 1 4.1478486 0 102
2009 19 361  6859 14  9 0 1  4.137884 0 102
2010 20 400  8000 15 10 1 1  3.951261 0 102
2011 21 441  9261 16 11 2 1 4.1478486 0 102
2012 22 484 10648 17 12 3 1 4.1520886 0 102
2013 23 529 12167 18 13 4 1  4.122022 0 102
2014 24 576 13824 19 14 5 1  4.116191 0 102
1990  0   0     0  0  0 0 0 3.2943344 1 103
1991  1   1     1  0  0 0 0 3.3155365 1 103
1992  2   4     8  0  0 0 0 3.5344756 1 103
1993  3   9    27  0  0 0 0 3.7463195 1 103
1994  4  16    64  0  0 0 0 3.8037305 1 103
1995  5  25   125  0  0 0 0 4.1457057 1 103
1996  6  36   216  1  0 0 0  4.073004 1 103
1997  7  49   343  2  0 0 0  4.046424 1 103
1998  8  64   512  3  0 0 0 3.9767094 1 103
1999  9  81   729  4  0 0 0  3.991047 1 103
2000 10 100  1000  5  0 0 0 3.9056516 1 103
2001 11 121  1331  6  1 0 0 4.0328646 1 103
2002 12 144  1728  7  2 0 0  4.271586 1 103
2003 13 169  2197  8  3 0 0 4.3809705 1 103
2004 14 196  2744  9  4 0 0 4.3761396 1 103
2005 15 225  3375 10  5 0 0  4.383377 1 103
2006 16 256  4096 11  6 0 0 4.5117874 1 103
2007 17 289  4913 12  7 0 0  4.636826 1 103
2008 18 324  5832 13  8 0 1 4.6824126 1 103
2009 19 361  6859 14  9 0 1 4.6859717 1 103
2010 20 400  8000 15 10 1 1 4.4817424 1 103
2011 21 441  9261 16 11 2 1  4.658967 1 103
2012 22 484 10648 17 12 3 1 4.6405506 1 103
2013 23 529 12167 18 13 4 1  4.596859 1 103
2014 24 576 13824 19 14 5 1 4.5692983 1 103
1990  0   0     0  0  0 0 0 3.5519204 1 104
1991  1   1     1  0  0 0 0  3.587767 1 104
1992  2   4     8  0  0 0 0 3.7967255 1 104
1993  3   9    27  0  0 0 0 3.9979155 1 104
1994  4  16    64  0  0 0 0 3.9500594 1 104
1995  5  25   125  0  0 0 0 4.2450676 1 104
1996  6  36   216  1  0 0 0  4.165329 1 104
1997  7  49   343  2  0 0 0  4.070019 1 104
1998  8  64   512  3  0 0 0 4.0524745 1 104
1999  9  81   729  4  0 0 0 4.0524745 1 104
2000 10 100  1000  5  0 0 0 3.9979155 1 104
2001 11 121  1331  6  1 0 0 4.0999994 1 104
2002 12 144  1728  7  2 0 0 4.3555245 1 104
2003 13 169  2197  8  3 0 0 4.4898725 1 104
2004 14 196  2744  9  4 0 0 4.5152617 1 104
2005 15 225  3375 10  5 0 0 4.5400224 1 104
2006 16 256  4096 11  6 0 0 4.6887236 1 104
2007 17 289  4913 12  7 0 0  4.809884 1 104
2008 18 324  5832 13  8 0 1   4.86635 1 104
2009 19 361  6859 14  9 0 1   4.88769 1 104
2010 20 400  8000 15 10 1 1 4.6793118 1 104
2011 21 441  9261 16 11 2 1  4.864388 1 104
2012 22 484 10648 17 12 3 1   4.86045 1 104
2013 23 529 12167 18 13 4 1  4.811956 1 104
2014 24 576 13824 19 14 5 1  4.778269 1 104
end
label values T_2008 T_2008lab
label def T_2008lab 0 "pre-economic crisis", modify
label def T_2008lab 1 "post-economic crisis", modify
label values exposure_cat exp_lab
label def exp_lab 0 "ref", modify
label def exp_lab 1 "A", modify

xtset group year
mixed outcome exposure_cat##c.(Time T_5 T_10 T_19) ///
    || group: Time, residuals(ar1, t(Time)) reml
   
//  CALCULATE MARGINS
local ats
forvalues y = 0/24 {
    foreach x in 5 10 19 {
        summ T_`x' if Time == `y', meanonly
        local T_`x' = r(mean)
    }
    local ats `ats' at(Time = `y' T_5 = `T_5' T_10 = `T_10' T_19 = `T_19')
}
margins exposure_cat, `ats' saving(margins_output, replace)


//  CLEAN UP MARGINS OUTPUT FOR GRAPHING
use margins_output, clear
rename _at Time
gen year = 1989 + Time
label var year Year
rename _m1 exposure_cat
label var exposure_cat "Exposure Category"
rename _margin outcome
label var outcome "Outcome"

keep outcome _ci* year exposure_cat

reshape wide outcome _ci_lb _ci_ub, i(year) j(exposure_cat)

graph twoway (line outcome* year, sort) (rcap _ci_lb0 _ci_ub0 year) ///
    (rcap _ci_lb1 _ci_ub1 year)

I have some concerns, however. In your example data I see variables Tsq and Tcube that have not previously been mentioned, and they appear to be quadratic and cubic terms in Time. You don't seem to use these in your regression command, so perhaps I shouldn't worry about them. But if you are using them in your regression command, it probably makes little sense to also use the T_5, T_10, and T_19 variables at the same time. If it does make sense to do that, then your model is going to need to be extremely complicated, and the code shown here would also have to be modified to accommodate it. But breaking up a time span of 24 years into 4 cubic polynomials does not strike me as a good idea.

There is also this unexplained variable T_2008, apparently an indicator for before and after the financial crisis: is this supposed to a covariate in the model? (Your original -mixed- command calls for three covariates, but this seems to be the only candidate for that role in your example data.)

[Clyde Schechter]

Also, instead of the "homebrew marginal spline" should I have created this using mkspline instead? Or another technique?

No, -mkspline- would have given you the same thing as your homebrew. I was just moaning about the difficulties of using spline with -margins-, not about how they were computed.

[Brennan RBratton]

Dear Clyde,
Please ignore those other time-related variables. I am not using those in this analysis and should have removed them from the sample, my apologies.

Thank you for explaining your comment about the spline and helping me with my quest to graph the results.
Be well