capture log close
log using prothrobin, replace text
set more off

********************************
** Marker data

insheet using prothro.dat, clear

* Divide prothrobin by 10
replace pro = pro/10

* save data in Stata format for later use
save junk, replace

********************************
** Read and expand survival data

clear
set mem 100m
insheet using prothros.dat, clear
* Define data as survival data for st commands
stset time, fail(death) id(id)

* Fit exponential model before expanding data (to compare with results after expanding data)
streg  treat, dist(exp) nohr

* Expand data
stsplit , at(failures)  riskset(set)

* Fit exponential model after expanding data
streg  treat, dist(exp) nohr

replace death = 0 if death==.
sort id time

* Create offset
sort id time
qui by id: gen w=cond(_n>1,time - time[_n-1],time)
gen lnw = ln(w)

* Fit exponential model using Poisson regression
poisson death treat, offset(lnw) 

* Mean-center time and express in years
summ time
local mean = r(mean)
gen t=(time-`mean')/365

* fracpoly poisson death t treat, offset(lnw) degree(3)
* third order polynomial got chosen

gen t2=t^2
gen t3=t^3
poisson death treat t t2 t3, offset(lnw)

********************************
** Merge survival data with marker data

gen var=2

append using junk
sort id var t
qui by id: gen occ=_n

* Mean-center time (same mean as for survival data)
replace t=(time-`mean')/365  if t==.
replace t2=0 if t2==.
replace t3=0 if t3==.

* resp = death if var==2 and resp = pro if var==1
gen resp = death
replace resp = pro if var==1

* Dummy variables for response types
tab var, gen(d)

* Interactions of dummy d2 with covariates for survival model
gen d2_t = d2*t
gen d2_t2 = d2*t2
gen d2_t3 = d2*t3
gen d2_treat = d2*treat


* No offset in marker model:
gen d2_lnw = cond(d2==1,lnw,0)


** Set up gllamm model: joint model for marker and survival

* equation for level-2 latent variable (latent marker)
eq lev2: d1 d2

* equation for level-3 random intercept
gen zero = 0
eq lev3: zero

* equation for fixed part of structural model
eq f1:  t treat 

* B matrix for adding level-3 random intercept to structural model
matrix B = (0,1\0,0)

* constraints for structural model - residual variance=0 and B-matrix contains fixed 1
cons def 1 [b1_2]_cons = 1
cons def 2 [occ1_1]d1 = 0

* starting values are estimates using nip(1 10)
matrix a=(8.1640282,-8.0389266,-.16640546,.09631052,-.00768988,-.15605308,.61928368, /*
 */ -.3532174,0,1.8409467,1,.16676205,-.79390748)
matrix coleq a=resp resp resp resp resp resp lns1 occ1_1l occ1_1 id2_1 b1_2  f1 f1
matrix colnames a=d1 d2 d2_t d2_t2 d2_t3 d2_treat _cons d2 d1 zero _cons  t treat


gllamm resp d1 d2 d2_t d2_t2 d2_t3 d2_treat, i(occ id) link(ident log) fam(gauss poiss) /*
  */ lv(var) fv(var) nocons offset(d2_lnw) eqs(lev2 lev3) nip(1 20) constr(1 2) geqs(f1) /*
  */ bmat(B) from(a) long trace

log close
exit