Writing a new program for plotting interactions

[Brian Burgess]

Hello,

I am an avid stata user and researcher who conducts tests of moderation regularly. When testing moderation in my field, it is customary to plot the interaction at +1 and -1 standard deviation above the mean. As far as I am aware, there is not a command that does this easily or automatically, so I generated the below code to make the plots I want.

In this code, the local functions are utilized to specify which variables go into each position (e.g., IV, DV, Moderator).

*-----CODE
reg dependentvar c.independentvar##c.moderatorvar

local mod "moderatorvar"
local iv "independentvar"
local dv "dependentvar"

summarize `mod'
local Modplus1sd = r(mean) + (1 * r(sd) ) // 1 specifies 1 SD above the mean.
local Modminus1sd = r(mean) - (1 * r(sd) )
local modmean=r(mean)
summarize `iv'
local IVplus1sd = r(mean) + r(sd)
local IVminus1sd = r(mean) - r(sd)
margins, dydx(`iv') at(`mod'=(`Modminus1sd' `modmean' `Modplus1sd')) // Test of Simple Slopes
margins, at(`iv'=(`IVminus1sd' `IVplus1sd') `mod'=(`Modminus1sd' `Modplus1sd')) atmeans asbalanced plot (ytitle(`dv') xtitle(`iv') title(Interaction of `iv' and `mod' on `dv')) nose

*---- END CODE


However, this block of code can become repetitive when I have a lot of plots to graph, it takes up a lot of code space, and it makes my do files look messy. Is it possible to write a new program and create a new command that runs this code. For example, it would be great if this program could be written into a command like this that outputs the plot and the test of simple slopes:


modplot depvar indepvar modvar, atsd(1)


Can someone help me with writing/learning to write a new command that runs this code?

Thank you!

[Clyde Schechter]

This is an implementation of what you want:

program define modplot
    syntax, iv(varname) dv(varname) moderator(varname) [atsd(real 1.0)]
    reg `dv' c.`iv'##c.`moderator'
    summarize `moderator'
    local moderatorplussds = r(mean) + (`atsd' * r(sd) ) // `atsd' specifies # sds above the mean.
    local moderatorminussds = r(mean) - (`atsd' * r(sd) )
    local moderatormean=r(mean)
    summarize `iv'
    local IVplussds = r(mean) + r(sd)
    local IVminussds = r(mean) - r(sd)
    margins, dydx(`iv') at(`moderator'=(`moderatorminussds' `moderatormean' `moderatorplussds')) // Test of Simple Slopes
    margins, at(`iv'=(`IVminussds' `IVplussds') `moderator'=(`moderatorminussds' `moderatorplussds')) ///
        atmeans asbalanced plot (ytitle(`dv') xtitle(`iv') ///
        title(Interaction of `iv' and `moderator' on `dv')) nose
    exit
end

I have modified your proposed syntax here. Instead of having a list of variables after -modplot- that have to be in the right order, I chose to make all three of those variables options (though they are required, not optional). A typical invocation of this command would be:

sysuse auto
modplot, iv(price) dv(mpg) moderator(rep78) atsd(2.5)

The -atsd()- option is not required. If you omit it, it will default to 1.0. Do note that, much as I would have liked to abbreviate -moderator()- to -mod()-, that cannot be done because -mod()- is already a built-in Stata function and attempting to use it in this context confuses the parser.

There are two ways you can use this code. You can just copy it into the do-file from which you plan to call it. Or you can save it as an .ado file in your ado/PERSONAL directory so that it will be callable from any do-file running on your system.

If I were really going to make a "production quality" do file out of this, there are other features I would add, such as allowing -if- and -in- options, and also allowing a list of covariate names to include in the regression. Probably I would also make the desired graph title an option to the -modplot- command rather than built in to the -margins- command, as I can envision use of a command like this in a do-file that repeatedly estimates this same model on different subsets of the data--which would lead to the plots all having the same title. Similarly, I would probably give -modplot- an option to specify a file in which to save the graph or store it in memory. Without this last, if the program is invoked more than once in the program, the graphs will all be overwritten, except for the last, unless each -modplot- call is followed by a -graph save- or -graph rename- command.

Anyway, this gives you the basic functionality you were looking for, and if you invest time in learning how to write programs, you can add on these other features later.

[Daniel Schaefer]

I typed something up myself, but it looks very much like what Clyde has already posted.

Do note that, much as I would have liked to abbreviate -moderator()- to -mod()-, that cannot be done because -mod()- is already a built-in Stata function and attempting to use it in this context confuses the parser.

I will just briefly say I chose to abbreviate moderator to md for this very reason. I think md works here because it follows the two letter abbreviation style of iv and dv, but this is purely a matter of personal preference. You could just say m() as shorthand if you think iv and dv are two letters because they abbreviate two words.

OP, you should also read -help program- and -help syntax- to learn more. Other useful resources are here and here.

[daniel klein]

This is an implementation of what you want:
Do note that, much as I would have liked to abbreviate -moderator()- to -mod()-, that cannot be done because -mod()- is already a built-in Stata function and attempting to use it in this context confuses the parser.

How?

. program mod
  1.     
.     version 18
  2.     
.     syntax , mod(varname)
  3.     
.     display "option mod() works and has been specified as `mod'"
  4.     
. end

. 
. sysuse auto , clear
(1978 automobile data)

. mod , mod(rep78)
option mod() works and has been specified as rep78

The same works as

. program mod
  1.     
.     version 18
  2.     
.     syntax , MODerator(varname)
  3.     
.     display "option mod() works and has been specified as `moderator'"
  4.     
. end

. 
. sysuse auto , clear
(1978 automobile data)

. mod , mod(rep78)
option mod() works and has been specified as rep78

Can you provide an example that fails?

[Clyde Schechter]

That is indeed strange. When writing #2, I initially tried to use -mod()- as the option, and Stata gave me "syntax error" originating out of the -syntax- statement itself. I was puzzled by that. But I noticed that in the do-file editor, mod was showing up in blue, reminding me that it is a defined Stata function. When I changed mod to moderator, the syntax error went away. I cannot explain why you are not encountering the same problem.

In any case, even if this is not universally a problem in all Stata installations, it is not good programming practice to use the names of existing functions for other purposes in Stata code.

[daniel klein]

In any case, even if this is not universally a problem in all Stata installations, it is not good programming practice to use the names of existing functions for other purposes in Stata code.

I think it should not be a problem in any Stata installation -- context matters. If syntax chokes on option names even though they are not declared reserved words, then that's a bug.

I tend to agree with your point about programming style, though I can definitely imagine exceptions.

[Brian Burgess]

Clyde Schechter This is incredibly helpful. My peers and I run this type of plot extremely often. Thank you so much for your help.

[Brian Burgess]

Thank you Daniel and Clyde, that code is super helpful and it works exactly as described. I have one more question.

My question is, is it possible to add additional elements to this code that allow for dichotomous Ivs or moderators? The above code specifies c. prior to each variable. However, having this means that the modplot syntax will only work for continuous predictors. Is there a way to state "if option continuous is specified, then run syntax with c., if else, include c. as a default?"

Alternatively, would you be able to point me to where I can find this in the programming manual that Daniel shared, or to examples of code that have nested options like this?

I tried looking at the materials that Daniel linked, but the code I generated didn't work. I am not sharing it here because I don't want to share faulty syntax for others to copy. Thank you!

[Daniel Schaefer]

Hi Brian, I believe something like the following should be simplest (based on Clyde's solution in #2). Please do double check to make sure I've got the two margins commands in the first code block (after the if) correct. Note that since the C is capitalized in the syntax command, you can abbreviate "continuous" as "c".

Clyde Schechter I'm curious why we need -exit- at the end of the program. Doesn't execution return automatically anyway?

capture program drop modplot
program define modplot
    syntax, iv(varname) dv(varname) moderator(varname) [atsd(real 1.0)] [Continuous]
    if(missing("`continuous'")){
        reg `dv' c.`iv'##i.`moderator'
        summarize `iv'
        local IVplussds = r(mean) + r(sd)
        local IVminussds = r(mean) - r(sd)
        qui levelsof `moderator', local(modlevels)
        margins, dydx(`iv') at(`moderator'=(`modlevels')) // Test of Simple Slopes
        margins, at(`iv'=(`IVminussds' `IVplussds') `moderator'=(`modlevels')) ///
            atmeans asbalanced plot (ytitle(`dv') xtitle(`iv') ///
            title(Interaction of `iv' and `moderator' on `dv')) nose
    }
    else{
        reg `dv' c.`iv'##c.`moderator'
        summarize `moderator'
        local moderatorplussds = r(mean) + (`atsd' * r(sd) ) // `atsd' specifies # sds above the mean.
        local moderatorminussds = r(mean) - (`atsd' * r(sd) )
        local moderatormean=r(mean)
        summarize `iv'
        local IVplussds = r(mean) + r(sd)
        local IVminussds = r(mean) - r(sd)
        margins, dydx(`iv') at(`moderator'=(`moderatorminussds' `moderatormean' `moderatorplussds')) // Test of Simple Slopes
        margins, at(`iv'=(`IVminussds' `IVplussds') `moderator'=(`moderatorminussds' `moderatorplussds')) ///
            atmeans asbalanced plot (ytitle(`dv') xtitle(`iv') ///
            title(Interaction of `iv' and `moderator' on `dv')) nose
    }
    exit
end

sysuse auto, clear
modplot, dv(weight) iv(length) moderator(trunk) continuous
modplot, dv(weight) iv(length) moderator(foreign)

[Clyde Schechter]

I'm curious why we need -exit- at the end of the program. Doesn't execution return automatically anyway?

We don't. It's entirely optional. Over my career, I've programmed in a number of different languages, and most of them require an explicit exit or return or something like that, so I have a habit of doing that. Also, my own "taste" in programming language design leans towards preferring that everything be as explicit in the code as possible, because I think it makes the code more transparent. So it's a matter of my style preferences. Feel free not to use -exit-.

Thank you for responding to #8. I'm having a really busy morning at my day job, and couldn't set aside the time needed to do that.

[Daniel Schaefer]

Of course, thank you for taking a moment to explain. I tend to agree that it is good practice to be explicit.

[Brian Burgess]

Daniel Schaefer and Clyde Schechter -- Thank you so much. This is awesome.

[Matthew Williams]

Hi Brian, I believe something like the following should be simplest (based on Clyde's solution in #2). Please do double check to make sure I've got the two margins commands in the first code block (after the if) correct. Note that since the C is capitalized in the syntax command, you can abbreviate "continuous" as "c".

Clyde Schechter I'm curious why we need -exit- at the end of the program. Doesn't execution return automatically anyway?

capture program drop modplot
program define modplot
syntax, iv(varname) dv(varname) moderator(varname) [atsd(real 1.0)] [Continuous]
if(missing("`continuous'")){
reg `dv' c.`iv'##i.`moderator'
summarize `iv'
local IVplussds = r(mean) + r(sd)
local IVminussds = r(mean) - r(sd)
qui levelsof `moderator', local(modlevels)
margins, dydx(`iv') at(`moderator'=(`modlevels')) // Test of Simple Slopes
margins, at(`iv'=(`IVminussds' `IVplussds') `moderator'=(`modlevels')) ///
atmeans asbalanced plot (ytitle(`dv') xtitle(`iv') ///
title(Interaction of `iv' and `moderator' on `dv')) nose
}
else{
reg `dv' c.`iv'##c.`moderator'
summarize `moderator'
local moderatorplussds = r(mean) + (`atsd' * r(sd) ) // `atsd' specifies # sds above the mean.
local moderatorminussds = r(mean) - (`atsd' * r(sd) )
local moderatormean=r(mean)
summarize `iv'
local IVplussds = r(mean) + r(sd)
local IVminussds = r(mean) - r(sd)
margins, dydx(`iv') at(`moderator'=(`moderatorminussds' `moderatormean' `moderatorplussds')) // Test of Simple Slopes
margins, at(`iv'=(`IVminussds' `IVplussds') `moderator'=(`moderatorminussds' `moderatorplussds')) ///
atmeans asbalanced plot (ytitle(`dv') xtitle(`iv') ///
title(Interaction of `iv' and `moderator' on `dv')) nose
}
exit
end

sysuse auto, clear
modplot, dv(weight) iv(length) moderator(trunk) continuous
modplot, dv(weight) iv(length) moderator(foreign)

Thank you for your useful code. How to modify the program if I want to add more independent variables? and cases where the interaction term is between two binary variables? Thanks!