st: graph the fitted values in a 3-way interaction with additional covariates

["Gentile, Elisabetta" <eg7@njit.edu>]

Dear StataListers,

I have a dataset of firms grouped by country and sector of operation.
I run a linear probability model with a three-way interaction, where
my dependent variable Y is obviously a binary indicator, my
independent variable X is continuos on [0,1], the modifier W is a
binary indicator, and Z is a categorical variable that splits the data
into three groups.

I want to graph the slopes for each of the interactions, so that if I
have the simple regression:

regress Y c.X##W##Z

then for example I can do:

twoway lfit Y X, by(W Z)

which generates six charts with the slopes for each group.

The problem is that my model has additional covariates, such as sector
and year fixed effects, and other continuous variables both at the
country and firm level.

How do I graph the slopes for each of the interactions for the model
including all the covariates?  The first thing that came to mind was
avplot, but it does not accept the by() option, so I am at a loss
here.  There must be a simple solution that I am not seeing!  Your
input is much appreciated.

Regards,
Elisa
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/