Hand-draw response surface

Here is a more elaborate example. Suppose we want to map out the
nonlinear relationship between two covariates and an outcome, say the
probability from a logistic estimation, by showing the probability ranges
in different colors and the covariates on the X- and Y-axes. Let's use the
familiar auto dataset to do just that.

o Create quadratic terms

We need some nonlinearity to make this interesting, so let's create quadratic
terms for weight and mileage (we will also scale the weight variable). (Click
the following 6 commands.)

. drop _all

. set mem 60m

. sysuse auto

. gen wt = weight / 1000

. gen wt2 = wt^2

. gen mpg2 = mpg^2

o Estimate model

Let's estimate a logistic model with foreign as the dependent variable
using our linear and quadratic terms as covariates.

. logistic foreign mpg mpg2 wt wt2

Now, it is time to produce our plot. Start by producing a lattice of values
for our two covariates. We will first find the range of each variable then
use the range and fillin commands to produce an evenly spaced
square lattice of points covering the range of the two variables

o Find range of data

. sum wt mpg

o Create lattice over range of 2 variables

. drop _all

. range wt 1.76 4.84 80

. range mpg 12 41 80

. fillin wt mpg

We need to recompute our quadratic terms for all of the points, then we can
predict the probability of being foreign at each of the points.

o Recreate quadratic terms

. gen wt2 = wt^2

. gen mpg2 = mpg^2

o Predict probability foreign for all lattice values

. predict prob_foreign

Let's cut these predicted probabilities into 4 ranges, each of which we will
plot in a different color

o Cut predicted probabilities into quartiles

. xtile qt = prob_foreign , nq(4)

We will keep only the variables necessary to draw the graph

. keep wt mpg qt

o Create corners of boxes around lattice points to fill the space

Despite requiring a few commands, we are just creating the points for little
squares that are centered at our lattice points.
Click here to run the next 14 commands in order. Or, if
you prefer a challenge, try clicking each individually in the correct
order.

. local dx = (1/2) * (4.84 - 1.76) / 79

. local dy = (1/2) * (41 - 12) / 79

. gen id = _n

. expand 6

. sort id

. by id: replace wt = wt - `dx' if _n==1
. by id: replace mpg = mpg - `dy' if _n<=2

. by id: replace wt = wt + `dx' if _n==2 | _n==3

. by id: replace mpg = mpg + `dy' if _n==3 | _n==4

. by id: replace wt = wt - `dx' if _n==4

. by id: replace wt = wt - `dx' if _n==5
. by id: replace mpg = mpg - `dy' if _n==5

. by id: replace wt = . if _n==6
. by id: replace mpg = . if _n==6

o Graph the boxes using a different color for each quartile

. twoway area mpg wt if qt==1 , nodropbase bc(navy) ||
area mpg wt if qt==2 , nodropbase bc(green) ||
area mpg wt if qt==3 , nodropbase bc(red) ||
area mpg wt if qt==4 , nodropbase bc(yellow)
scheme(intensity)
ylabel(20 30 40) xlabel(2 3 4)
plotregion(margin(zero))
legend(subtitle(Predicted probability foreign)
order(1 "(0,.25]" 2 "(.25,.5]"
3 "(.5,.75]" 4 "(.75,1)" )
)

We can give the two dimensions equal importance by adding the aspect(1)
option.

. twoway area mpg wt if qt==1 , nodropbase bc(navy) ||
area mpg wt if qt==2 , nodropbase bc(green) ||
area mpg wt if qt==3 , nodropbase bc(red) ||
area mpg wt if qt==4 , nodropbase bc(yellow)
scheme(intensity)
ylabel(20 30 40) xlabel(2 3 4)
plotregion(margin(zero))
legend(subtitle(Predicted probability foreign)
order(1 "( 0,.25]" 2 "(.25,.5]"
3 "(.5,.75]" 4 "(.75,1)" )
)
, aspect(1)

>> Our intensity scheme

<< index >>