* Clustered Outcomes Chapter
drop _all

input ship op1 op2  per1 per2 months incident
1 60 64 60 74   127   0
1 60 64 75 79    63   0
1 65 69 60 74  1095   3
1 65 69 75 79  1095   4
1 70 74 60 74  1512   6
1 70 74 75 79  3353  18
1 75 79 60 74     0   0
1 75 79 75 79  2244  11
2 60 64 60 74 44882  39
2 60 64 75 79 17176  29
2 65 69 60 74 28609  58
2 65 69 75 79 20370  53
2 70 74 60 74  7064  12
2 70 74 75 79 13099  44
2 75 79 60 74     0   0
2 75 79 75 79  7117  18
3 60 64 60 74  1179   1
3 60 64 75 79   552   1
3 65 69 60 74   781   0
3 65 69 75 79   676   1
3 70 74 60 74   783   6
3 70 74 75 79  1948   2
3 75 79 60 74     0   0
3 75 79 75 79   274   1
4 60 64 60 74   251   0
4 60 64 75 79   105   0
4 65 69 60 74   288   0
4 65 69 75 79   192   0
4 70 74 60 74   349   2
4 70 74 75 79  1208  11
4 75 79 60 74     0   0
4 75 79 75 79  2051   4
5 60 64 60 74    45   0
5 60 64 75 79     0   0
5 65 69 60 74   789   7
5 65 69 75 79   437   7
5 70 74 60 74  1157   5
5 70 74 75 79  2161  12
5 75 79 60 74     0   0
5 75 79 75 79   542   1
end
gen co_60_64 = (op1==60 & op2==64)
gen co_65_69 = (op1==65 & op2==69)
gen co_70_74 = (op1==70 & op2==74)
gen co_75_79 = (op1==75 & op2==79)
gen op_60_74 = (per1==60 & per2==74)
gen op_75_79 = (per1==75 & per2==79)
gen double exposure = ln(months)
                                                                                                                             
glm incident op_75_79 co_65_69 co_70_74 co_75_79, family(poiss) link(log) cluster(ship) offset(exposure) eform nolog

xi: glm incident op_75_79 co_65_69 co_70_74 co_75_79 i.ship, family(poiss) link(log) offset(exposure) eform nolog

xtpoisson incident op_75_79 co_65_69 co_70_74 co_75_79, fe i(ship) offset(exposure) eform nolog

xtpoisson incident op_75_79 co_65_69 co_70_74 co_75_79, re normal i(ship) offset(exposure) eform nolog

xtpoisson incident op_75_79 co_65_69 co_70_74 co_75_79, re i(ship) offset(exposure) eform nolog

gllamm incident op_75_79 co_65_69 co_70_74 co_75_79, link(log) family(poisson) i(ship) offset(exposure) eform

xtpoisson incident op_75_79 co_65_69 co_70_74 co_75_79, pa i(ship) offset(exposure) vce(robust) eform nolog

set seed 1
set more off
set obs 500
gen x1 = uniform()*10
gen x2 = uniform()*5
gen x3 = uniform()
gen x4 = uniform()
gen x5 = uniform() 
gen err = invnormal(uniform())
gen y =  1*x1 + 2*x2 + 3*x3 + 4*x4 + 5 + err
gen a1 = x3 + 0.25*invnormal(uniform())
gen a2 = x3 + 0.25*invnormal(uniform())
gen b1 = x4 + 0.25*invnormal(uniform())
gen b2 = x4 + 0.25*invnormal(uniform())

glm y x1 x2 x3 x4, noheader

glm y x1 x2 a1 b1, noheader
glm y x1 x2 a2 b2, noheader

simex (y=x1 x2) (w3: a1 a2) (w4: b1 b2), brep(299) seed(12345)
graph export simexplot.eps, replace