* Nonlogit Chapter
drop _all

range eta -3 3 100
gen logit  = exp(eta)/(1+exp(eta))
gen probit = normprob(eta)
gen clog   = 1-exp(-exp(eta))
gen logl   = exp(-exp(-eta))
capture program drop doit1
program define doit1
	version 6.0
	graph logit probit eta, ylab(0,.25,.5,.75,1) xlab(-3,-2,-1,0,1,2,3) /*
		*/ saving(sigmoid, replace) s(o.) c(.l) t1(" ") t2("logit(symbols)    probit(line)")
end
capture program drop doit2
program define doit2
	version 6.0
	graph clog logl eta, ylab(0,.25,.5,.75,1) xlab(-3,-2,-1,0,1,2,3) /*
		*/ saving(sigmoid2, replace) s(o.) c(.l) t1(" ") t2("cloglog(symbols)    loglog(line)")
end
capture program drop doit3
program define doit3
	version 6.0
	graph logit probit eta eta if eta>=-1 & eta<=2, ylab(0,.25,.5,.75,1) /*
		*/ xlab(-1,0,1,2) saving(sigmoid3, replace) s(op.) c(..l) t1(" ") /*
		*/ t2("logit(o)    probit(+)   identity(line)") yline(0,1)
end
doit1
graph export sigmoid.eps, replace
doit2
graph export sigmoid2.eps, replace
doit3
graph export sigmoid3.eps, replace



drop _all
set obs 101
gen eta = -3 + 6*(_n-1)/100
gen logit   = exp(eta)/(1+exp(eta))
gen probit  = normprob(eta)
gen loglog  = exp(-exp(-eta))
gen cloglog = 1-exp(-exp(eta)) 
gen identity= eta 
twoway (scatter logit eta, c(l) s(s) xlabel(-3(1)3)) (scatter probit eta, c(l) s(o) ylabel(0(.25)1))
graph export probitlogit.eps, replace

twoway (scatter loglog eta, c(l) s(s) xlabel(-3(1)3)) (scatter cloglog eta, c(l) s(o) ylabel(0(.25)1))
graph export loglogcloglog.eps, replace

keep if eta >= -1
keep if eta <= 2

twoway  (scatter logit eta , c(l) s(s) xlabel(-1(1)2) yline(1) ylabel(0(.25)1)) (scatter probit eta, c(l) s(o) ylabel(0(.25)1) yline(0)) (scatter identity eta, ylabel(0(.25)1) msymbol(none) c(l)), legend(rows(1))
graph export probitlogitidentity.eps, replace

use heart01, clear
glm death anterior hcabg kk2 kk3 kk4 age2-age4, family(bin) link(probit) nolog

tab death if e(sample)

glm death anterior hcabg kk2 kk3 kk4 age2-age4, family(bin) link(loglog) nolog

glm death anterior hcabg kk2 kk3 kk4 age2-age4, family(bin) link(loglog) notable nolog
glm death anterior hcabg kk2 kk3 kk4 age2-age4, family(bin) link(probit) notable nolog
glm death anterior hcabg kk2 kk3 kk4 age2-age4, family(bin) link(logit) notable nolog
glm death anterior hcabg kk2 kk3 kk4 age2-age4, family(bin) link(cloglog) notable nolog


drop _all
input age total menarche
 9.21  376    0
10.21  200    0
10.58   93    0
10.83  120    2
11.08   90    2
11.33   88    5
11.58  105   10
11.83  111   17
12.08  100   16
12.33   93   29
12.58  100   39
12.83  108   51
13.08   99   47
13.33  106   67
13.58  105   81
13.83  117   88
14.08   98   79
14.33   97   90
14.58  120  113
14.83  102   95
15.08  122  117
15.33  111  107
15.58   94   92
15.83  114  112
17.58 1049 1049
end

gen proportion = menarche/total
label var proportion "Observed prop."

glm menarche age, family(bin total) nolog
predict double logithat
replace logithat = logithat/total
label var logithat "Predicted prop. (logit)"

glm menarche age, family(bin total) link(probit) nolog

predict double probithat
replace probithat = probithat/total
label var probithat "Predicted prop. (probit)"
twoway (scatter probithat age, s(s)) (scatter logithat age, s(o)) (scatter proportion age, s(x)), legend(rows(1))

graph export menarchelogitprobit.eps, replace

glm menarche age, family(bin total) link(logit) nolog notable
glm menarche age, family(bin total) link(probit) nolog notable
glm menarche age, family(bin total) link(loglog) nolog notable
glm menarche age, family(bin total) link(cloglog) nolog notable


drop _all
input trt rep depth damage examine tpaper
1 1 1 120 187 1
1 2 1 145 181 1
1 3 1 123 184 1
1 1 2  60 184 2
1 2 2  85 191 2
1 3 2 113 171 2
1 1 3  35 179 3
1 2 3  23 181 3
1 3 3  33 179 3
1 1 4   5 178 4
1 2 4  40 184 4
1 3 4  18 171 4
1 1 5  66 182 5
1 2 5 104 186 5
1 3 5 106 181 5
2 1 1  97 187 6
2 2 1 128 173 6
2 3 1 112 179 6
2 1 2  85 186 7
2 2 2  95 167 7
2 3 2 137 194 7
2 1 3  72 190 8
2 2 3  89 184 8
2 3 3  92 192 8
2 1 4  49 188 9
2 2 4  50 187 9
2 3 4  65 186 9
2 1 5 132 179 10
2 2 5 122 182 10
2 3 5 120 179 10
3 1 1 138 187 11
3 2 1 169 184 11
3 3 1 163 175 11
3 1 2 156 176 11
3 2 2 169 190 11
3 3 2 159 175 11
end

tab tpaper, gen(T)
tab rep, gen(R)

gen prop = damage/examine
label var tpaper "Treatment"
label var prop   "Observed proportion"
egen mprop = mean(prop), by(tpaper)
twoway scatter prop tpaper, xlabel(1(1)11)
graph export carrot.eps, replace

local sv = 0 

* This is the exact answer quoted in Famoye's paper
glm damage R2 R3 T1 T2 T3 T4 T5 T6 T7 T8 T9 T10, family(bin examine) link(clog) 
predict double glm
replace glm = glm/examine
egen mglm = mean(glm), by(tpaper)

matrix b0 = e(b)
local ll0 = `e(ll)'
matrix z0 = (0) 
matrix b0 = b0, z0

capture program drop doit_ll
program define doit_ll
	version 9.0
	args todo b lnf 
	tempvar  eta pi 
	tempname phi
	mleval `eta' = `b', eq(1)
	mleval `phi' = `b', eq(2) scalar
	local y  "$ML_y1"   /* Dep var name is in global macro by ml  */
	local M  "examine"  /* Specify directly -- problem specific   */
	gen double `pi'    = 1-exp(-exp(`eta'))  /* Cloglog link here */
	mlsum `lnf' = log(`M') + `y'*log(`pi') +              /*
		*/ (`M'+`phi'*`y'-`y')*log(1+`phi'*`pi'-`pi') - /*
		*/ (`M'+`phi'*`y')*log(1+`phi'*`pi') -          /*
		*/ log(`M'+`phi'*`y') +                         /*
		*/ lngamma(`M'+`phi'*`y'+1) - lngamma(`y'+1) -  /*
		*/ lngamma(`M'+`phi'*`y'-`y'+1)
end

ml model d0 doit_ll (gbr: damage = R2 R3 T1 T2 T3 T4 T5 T6 T7 T8 T9 T10) /phi, search(off) init(b0, copy) max iter(50)
ml display, title("Generalized binomial regression")

predict double xb
gen double gb = 1-exp(-exp(xb))
egen mgbr = mean(gb), by(tpaper)

list mprop mglm mgbr