* Author: Nicola Orsini
* Workshop: Going beyond odds ratios: understanding what else logistic regression can offer
* Date: May 24, 2016
* Motivating example by Sara Fritzell.
/* Data: cross-sectional study
y = (0,1) 0 "Good health" 1 "Poor Health"
x = (0,1) 0 "Two parents" 1 "Single-mother"
z = (1,2,..., 8) Number of children
OR xy = 2
Compared to two-parents family, single mothers had a 2-fold
increase odds of poor health regardless of the number of children.
OR zy = 0.8
Every 1 additional children was associated with a 20% lower odds
of poor health, no matter what is the family type.
Baseline risk of being in poor health with 1 child (br) = 20%
To make sure Pr(y=1|x=1) = 0.2, just subtract 1 to x
log(p/(1-p)) = b0 + b1*(x-1) + b2*z
odds = exp(b0 + b1*(x-1) + b2*z)
risk = odds/(1+odds)
b0 = is the log odds of being in poor health for a two-parent
family type and a mother with 1 children
*/
scalar baseline_odds = .2/(1-.2)
scalar b0 = log(baseline_odds)
scalar b1 = log(2)
scalar b2 = log(0.8)
* Graph of the risk function in the two groups
twoway (function risk0 = invlogit(b0 + b1*0 + b2*(x-1)), range(1 8) lc(red)) ///
(function risk1 = invlogit(b0 + b1*1 + b2*(x-1)), range(1 8) lc(blue)) , ///
legend(label(1 "Two Parents") label(2 "Single mother") ring(0) pos(1) col(1)) plotregion(style(none)) ///
ytitle("Risk of poor self-reported health") xtitle("Number of children") ///
ylabel(#10, angle(horiz) format(%3.2f)) xlabel(0(1)8) name(risk, replace)
* Graph of the odds function (on the log-scale) in the two groups
twoway (function odds0 = exp(b0 + b1*0 + b2*(x-1)), range(1 8) lc(red)) ///
(function odds1 = exp(b0 + b1*1 + b2*(x-1)), range(1 8) lc(blue)) , ///
legend(label(1 "Two Parents") label(2 "Single mother") ring(0) pos(1) col(1)) plotregion(style(none)) ///
ytitle("Odds of poor self-reported health") xtitle("Number of children") ///
ylabel(#10, angle(horiz) format(%3.2f)) xlabel(0(1)8) yscale(log) name(odds, replace)
* What is the marginal effect among single mother with 1 child?
scalar p = invlogit(b0 + b1*1 + b2*(1-1))
display p*(1-p)*b1
* What is the marginal effect among two parents mothers with 1 child?
scalar p = invlogit(b0 + b1*0 + b2*(1-1))
display p*(1-p)*b1
* What is the marginal effect among single mother with 2 children?
scalar p = invlogit(b0 + b1*1 + b2*(2-1))
display p*(1-p)*b1
* What is the marginal effect among two parents mothers with 2 children?
scalar p = invlogit(b0 + b1*0 + b2*(2-1))
display p*(1-p)*b1
* To estimate an average marginal effect we would need the actual exposure distribution in the sample.
* Real world data from NEJM, 2005
use http://www.imm.ki.se/biostatistics/data/hyponatremia, clear
* average marginal change for a continuous predictor
* Univariable model
logit nas135 wtdiff
* Estimate the conditional probability of the outcome at any value of the exposure distribution
gen p = invlogit(_b[_cons]+_b[wtdiff]*wtdiff)
* Estimate the marginal effect at any value of the exposure distribution
gen mc = p*(1-p)*_b[wtdiff]
* Average the marginal effects
tabstat mc
* Very convenient to use the post-estimation command (after logistic)
margins , dydx(wtdiff)
* Graph the distribution of the marginal effects
tw (histogram mc, percent width(.02) fcolor(white)) ///
, ylabel(#10, angle(horiz) format(%2.0f)) ///
plotregion(style(none)) xtitle("Marginal effect of weight change") ///
xlabel(0(.02).2, format(%3.2f) ) xline(0.068, lc(black) lp(-)) text(20 .068 "Average Marginal Effect", place(3))
* Multivariable model
* MEM (Marginal effects at Means)
logistic nas135 wtdiff female
su wtdiff female
capture drop p
capture drop mc
gen p = invlogit(_b[_cons]+_b[wtdiff]*-.689473 +_b[female]*.3472527 )
gen mc = p*(1-p)*_b[wtdiff]
su mc
margins , dydx(wtdiff) atmeans
* MEM (Marginal effects at Representative values), example men who did not change weight
capture drop p
capture drop mc
gen p = invlogit(_b[_cons]+_b[wtdiff]*0 +_b[female]*0)
gen mc = p*(1-p)*_b[wtdiff]
su mc
margins , dydx(wtdiff) at(wtdiff=0 female =0)
* AME (Average marginal effects)
logistic nas135 wtdiff female
* Average marginal effect leaving all predictors at observed
capture drop p mc
gen p = invlogit(_b[_cons]+_b[wtdiff]*wtdiff +_b[female]*female) if e(sample)
gen mc = p*(1-p)*_b[wtdiff]
su mc
margins , dydx(wtdiff)
* Average marginal effect fixing one predictor and leaving others at observed
capture drop p mc
gen p = invlogit(_b[_cons]+_b[wtdiff]*wtdiff +_b[female]*0) if e(sample) & female ==0
gen mc = p*(1-p)*_b[wtdiff]
su mc
margins if female == 0 , dydx(wtdiff) at(female=0)
capture drop p mc
gen p = invlogit(_b[_cons]+_b[wtdiff]*wtdiff +_b[female]*1) if e(sample) & female ==1
gen mc = p*(1-p)*_b[wtdiff]
su mc
margins if female == 1, dydx(wtdiff) at(female=1)
* Average Discrete Change for a continuous predictor and 1 unit increase.
logistic nas135 wtdiff i.female
capture drop adc
gen adc = invlogit(_b[_cons]+_b[wtdiff]*(wtdiff+1) +_b[1.female]*female) - ///
invlogit(_b[_cons]+_b[wtdiff]*wtdiff +_b[1.female]*female) if e(sample)
su adc
* Average Discrete Change for a binary predictor
logistic nas135 wtdiff i.female
capture drop adc
gen adc = invlogit(_b[_cons]+_b[wtdiff]*wtdiff +_b[1.female]*1) - ///
invlogit(_b[_cons]+_b[wtdiff]*wtdiff +_b[1.female]*0) if e(sample)
su adc
margins, dydx(female)