version 13.1
* Set up data
webuse nhanes2f, clear
keep if !missing(diabetes, black, female, age)
// Note: In later versions of Stata, some of these labels
// are already defined. When that is so, the capture 
// prefix keeps the program from aborting with an errorlabel define black 0 "nonBlack" 1 "black"
capture label define female 0 "male" 1 "female"
capture label values female female
capture label values black black
quietly logit diabetes i.black i.female age c.age#c.age, nolog
margins, dydx(*)

* Basics
logit diabetes i.black i.female age c.age#c.age, nolog
est store m1
mcp age, show
mcp age, ci
quietly margins, at(age=(20(1)74))
marginsplot, noci

* Select ranges of values

mcp age, var1(20) show
mcp age, at1(20(1)100) show
mcp age, at1(% 10 (1) 90) show

* Transformed X
generate logage = log(age)
logit diabetes i.black i.female logage, nolog
est store m2
lrtest m1 m2, stats
mcp logage
summarize age
range w1 r(min) r(max) 20
generate logw1 = log(w1)
fre w1 logw1, tab(5)
mcp age (logage), var1(w1 (logw1)) show

* Plot multiple groups simultaneously
quietly logit diabetes i.black i.female age c.age#c.age, nolog
mcp age black, plotopts(scheme(sj))
mcp age black, ci
mcp age black, ci plotopts(ycommon)
* Use a continuous variable for X2, but specify plotting values
quietly logit diabetes i.black i.female age c.age#c.age c.weight, nolog
mcp weight age, at1(% 1 (1) 99) at2( 20 45 70) plotopts(scheme(sj)) show
mcp weight age, var1(30) at2( 20 45 70) plotopts(scheme(sj)) show

* mcp with OLS regression & nonlinear terms
sysuse auto, clear
quietly reg price mpg c.mpg#c.mpg
mcp mpg

* Appendix: Analysis with Spline Models
* Thanks to Patrick Royston for providing this example!
set scheme sj
/*
	Spline example, 4 knots at 5 35 65 95 centiles of age.
*/
global stub ex3b
webuse nhanes2f, clear
keep if sex==1
gen map = (bpsystol + 2 * bpdiast)/3
gen bmi = weight / (height/100)^2

// Generate splines for age (3 basis variables, zage1-zage3)
mkspline zage = age, cubic nknots(4)
matrix knots = r(knots)
local knots
forvalues j = 1 / 4 {
	local k = knots[1, `j']
	local knots `knots' `k'
}
regress map zage1 zage2 zage3 i.race bmi hgb race#c.(zage1 zage2 zage3) c.bmi#c.(zage1 zage2 zage3) c.hgb#c.(zage1 zage2 zage3)

// Create spline-transformed plotting ages
local ages 20 22 24 26 28 30 35 40 45 50 55 60 70
gen ages = .
tokenize `ages'
local i 1
while "``i''" != "" {
	qui replace ages = ``i'' in `i'
	local ++i
}
mkspline wage = ages, cubic knots(`knots')

local stuff lp(l - _ _- -.) lwidth(medthick ..)

marginscontplot age (zage1 zage2 zage3), var1(ages(wage1 wage2 wage3)) ci ///
 plotopts(title("(a) Age only", placement(west)) leg(off) name(g1, replace))

marginscontplot age (zage1 zage2 zage3) race, var1(ages(wage1 wage2 wage3)) ///
plotopts(`stuff' title("(b) Age x race", placement(west)) legend(label(1 "White") label(2 "Black") label(3 "Other") row(1)) name(g2, replace))

marginscontplot age (zage1 zage2 zage3) bmi, at2(20(5)40) var1(ages(wage1 wage2 wage3)) ///
 plotopts(`stuff' title("(c) Age x BMI", placement(west)) legend(label(1 "20") label(2 "25") label(3 "30") label(4 "35") label(5 "40") row(1)) name(g3, replace))

marginscontplot age (zage1 zage2 zage3) hgb, at2(13(1)16) var1(ages(wage1 wage2 wage3)) ///
plotopts(`stuff' title("(d) Age x haemoglobin", placement(west)) legend(label(1 "13") label(2 "14") label(3 "15") label(4 "16") label(5 "17") row(1)) name(g4, replace))

graph combine g1 g2 g3 g4, imargin(small) saving($stub, replace)
graph export $stub.eps, replace
graph drop g1 g2 g3 g4