/* Liver data from Nelson et al with gamma random effects */

data liver;
input INST SURVT treat heart;
log_s = log(SURVT);
cens=1;
datalines;
1 23.286 1 1
1 6.429 0 0
1 26.857 1 0
1 11.143 0 1
2 3.857 0 1
2 9.000 1 0
2 8.714 0 0
3 1.143 0 1
3 23.143 1 0
3 2.571 1 0
3 2.857 0 1
4 76.429 1 0
4 35.857 1 1
4 25.857 0 0
4 52.286 0 0
5 25.143 1 1
5 25.143 0 0
6 29.286 1 1
6 28.857 0 0
7 1.857 0 1
7 14.286 1 0
run;

*proc nlmixed data=liver method=GAUSS NOAD ;
proc nlmixed data=liver;
parms theta_1 1 b0 -2.8 b_trt -.54 b_hrt .18;
/*starting values under independence for b0,b_trt,b_hrt */
p_i = CDF('NORMAL',a_i) ; /* These next three lines create */
g_i2 = GAMINV(p_i,1/theta_1) ; /* gamma random effect g_i */
g_i = theta_1*g_i2 ;
xb = log_s + b0 + b_trt*treat + b_hrt*heart + log(g_i) ;
lambda = exp(xb);
model cens ~ poisson(lambda);
random a_i ~ normal(0,1) subject=inst;
run;