clear
set obs 200
generate age = runiformint(18, 99)

* Cutpoints and Multipliers
* We can define a logit curve in terms of its
*    cutpoint and multiplier.
* The cutpoint is just the point where the
*    probability is 0.5, i.e. where the odds 
*    are 1, i.e. where the log-odds, the linear
*    predictor in GLM terms, is 0.
* The multiplier is difficult to explain in
*    terms of probability, because it changes.
*    However, the multiplier is constant in terms
*    of odds, and it becomes a constant addition
*    in terms of log-odds.

* Let's illustrate with a simulation.
* In this example we will model Wisdom.  Either
*    you are wise, or you are not.  Based on my
*    experience, I would say it is like someone
*    flips on a light switch.
* By age 55, about half of us have become wise
*    (sorry, kids), so this is our cutpoint.
* The probability that you are wise increases as
*    you get older - there is an age multiplier
*    which is greater than 1 in terms of odds, or
*    greater than 0 in terms of log-odds.  (In 
*    terms of probability, the amount of increase
*    varies, depending on the age.)
* For this example, we will take the increase to
*    be about 0.2 per year in log-odds, or about
*    a 1.22 per year multiplier in odds.  (This implies
*    the odds of being wise doubles each 3.5 years.)

* One approach to simulating data begins with our
*    theoretical log-odds equation.
* If our log-odds, the linear predictor, is
*    B0 + B1*age, then
*    B1 = 0.2, given above.
*    B0 is given by B0 + 0.2*55 = 0, so
*    B0 = -0.2*55.
*    (This is the log-odds when age = 0.)

* In Stata, we can take these log-odds and convert
*    them into probabilities.  (We could write this
*    a couple of ways.)
generate wisel = invlogit(-0.2*55 + 0.2*age)
* Draw the theoretical curve we wish to simulate.
twoway (line wisel age, sort), name(g1)
* Draw numbers from a random uniform distribution,
*    compare to our theoretical distribution, and
*    generate random 0's and 1's from our model.
generate wise = wisel > runiform()
* Graph both the data and the theoretical distribution.
twoway (line wisel age, sort) (scatter wise age), name(g2)
* Now turn around and estimate the model based on
*    this particular simulation.
logit wise age

* What does the constant mean?
* The constant gives us the log-odds of being wise
*    at the age of zero, a very small number.
display 1/(1+exp(-_b[_cons]))
* The margins command, with all independent variable
*    values specified, gives us the conditional margin,
*    plus an estimate of the standard error.
margins, at(age=(0))

* We can also investigage the probability of being wise
*    at age 55.
display 1/(1+exp(-(_b[_cons]+55*_b[age])))
margins, at(age=(55))

* Margins makes it easy to graph the estimated model
*    (albeit with delta-method standard errors).
quietly margins, at(age=(18(5)99))
marginsplot, name(g3)

generate female = runiformint(0, 1)
generate wisef = invlogit(0.2*(age-55+female*(5)))

preserve
separate wisef, by(female)
line wisef? age, sort
restore

generate wiseg = wisef > runiform()
logit wiseg age i.female

margins female, at(age=(50,55))
marginsplot

* Cutpoints
local fcut = -(_b[_cons]+_b[1.female])/_b[age]
display `fcut'
local mcut = -(_b[_cons]+_b[0.female])/_b[age]
display `mcut'

margins female, at(age=(`fcut' `mcut'))

quietly margins female, at(age=(18(5)99))
marginsplot

generate wisefi = invlogit(0.1*(1+female)*(age-70+female*(20)))
preserve
separate wisefi, by(female)
line wisefi? age, sort
restore

generate wisegi = wisefi > runiform()
logit wisegi c.age##female

* Cutpoints
local fcut = -(_b[_cons]+_b[1.female])/(_b[age]+_b[1.female#c.age])
display `fcut'
local mcut = -(_b[_cons]+_b[0.female])/(_b[age]+_b[0.female#c.age])
display `mcut'
* Inflection point
local equal = -_b[1.female]/_b[1.female#c.age]
display `equal'

margins female, at(age=(`fcut' `mcut' `equal'))
quietly margins female, at(age=(18(5)99))
marginsplot