cap log close
clear
log using StataLog, replace 

include "http://192.38.117.59/~linearpredictors/datafiles/readFever.do"

/************************************************
Logistic regressions giving log(odds). Here
alcogroup==0 is automatically selected as
reference group.
*************************************************/

xi: logistic death i.alcogroup
matrix logor = e(b)

quietly xi: logit death i.alcogroup, noconstant 
matrix logodds = e(b) 

quietly mean alco, over(alcogroup)
matrix meanMat = e(b) 

/************************************************
Generating points to be plotted
*************************************************/

local logodds_1 = logor[1,5] // the constant term
forvalues i = 1/5 {
  local j = `i'+1
  local logodds_`j' = logodds[1,`i']
  local means_`i' = meanMat[1,`i']
  local points = "`points' `logodds_`i'' `means_`i''"
}

/************************************************
Calculating regression line
*************************************************/

egen alcomeans = mean(alco), by(alcogroup)  
xi: logit death alcomeans

matrix line = e(b) 
local slope = line[1,1]
local intercept = line[1,2]

/************************************************
The plot
*************************************************/
  
twoway scatteri `points' || function y=`slope'*x+`intercept', /*
*/ range(0 5) xtitle("Drinks per week") ytitle("log(hazard ratio)") /*
*/ legend(off)
log close