*analysis with censoring: analysis 4 in validation study example
*Risk periods: 13-41 days after each dose

use intuss, clear
sort indiv eventday

*keep only first event
duplicates drop indiv, force
replace indiv = _n
set more off

*remove exposures whose risk period falls before the start of observation
gen before = 1 if agep1 + 41 < cutp1
replace before = 2 if agep2 + 41 < cutp1
replace before = 3 if agep3 + 41 < cutp1
replace agep1 = . if before==3
replace agep1 = agep3 if before==2
replace agep1 = agep2 if before==1
replace agep2 = . if before!=.
replace agep2 = agep3 if before==1
replace agep3 = . if before!=.
drop before

local nindivs = 500

*generate the age group cut points
gen cutp3 = 57
gen cutp4 = 87
gen cutp5 = 117
gen cutp6 = 148
gen cutp7 = 179
gen cutp8 = 210
gen cutp9 = 241
gen cutp10 = 272
gen cutp11 = 303
gen cutp12 = 334

*number of age group cut points
local nage = 10

*generate the exposure group cut points
gen cutp13 = agep1 + 13
gen cutp14 = agep1 + 41

gen cutp15 = agep2 + 13
gen cutp16 = agep2 + 41

gen cutp17 = agep3 + 13
gen cutp18 = agep3 + 41

drop agep1 agep2 agep3

*number of exposure group cut points
local nexgr = 6

foreach i of numlist 13/17{
local j = `i'+1
replace cutp`i' = cutp`j' if cutp`i' > cutp`j' 
}

keep indiv eventday cutp*
compress

local a=cutp3
local b=`nage'+2
while `b'>3{
local c = cutp`b'
local d = ",`c'`d'"
local b = `b' - 1
}
local a = "`a'`d'"

foreach var of varlist cutp*{
replace `var' = cutp1 if `var' < cutp1
replace `var' = cutp2 if `var' > cutp2
}

compress

* expand data and count number of events in interval
sort indiv eventday
reshape long cutp, i(indiv eventday) j(type)
sort indiv eventday cutp type
by indiv: generate int nevents = 1 if eventday > cutp[_n-1]+0.5 & eventday <= cutp[_n]+0.5
collapse (sum) nevents, by(indiv cutp type)

* interval length
by indiv: generate int interval = cutp[_n] - cutp[_n-1]

* age groups
by indiv: generate int agegr = irecode(cutp, `a')

* number exposure intervals - 1, 2, 3, 4, 5... 
local e = 2+`nage'
generate exgr = type-`e'-1 if type>`e'
while `e'>0{
by indiv: replace exgr = exgr[_n+1] if exgr==.
local e = `e'-1
}
replace exgr = `nexgr' if exgr==.

* tidy up
drop cutp* type
drop if interval ==0 | interval==.
gen loginterval = log(interval)
drop interval

gen bootindiv = indiv
gen eventweight = 0
bysort indiv: egen exev = max(nevents*exgr)
compress
save intuss_cens, replace

* data after 3rd exposure

use intuss_cens, clear
drop if exgr < 5
drop if exev < 5
recode exgr (5=3) (6=0)
replace indiv = indiv + (`nindivs'*3)
compress
save 3rdexpo, replace

* data after 2nd exposure

use intuss_cens, clear
drop if exgr < 3
drop if exev < 3
replace eventweight = 3 if nevents == 1 & exgr == 5
recode exgr (3=2) (4=0) (5=0) (6=0)
replace indiv = indiv + (`nindivs'*2)
compress
save 2ndexpo, replace

* data after 1st exposure

use intuss_cens, clear
drop if exgr < 1
drop if exev < 1
replace eventweight = 2 if nevents == 1 & exgr == 3
replace eventweight = 3 if nevents == 1 & exgr == 5
recode exgr (1=1) (2=0) (3=0) (4=0) (5=0) (6=0)
replace indiv = indiv + `nindivs'
compress
save 1stexpo, replace

* all data - to estimate age effects

use intuss_cens, clear
replace eventweight = 1 if nevents == 1 & exgr == 1
replace eventweight = 2 if nevents == 1 & exgr == 3
replace eventweight = 3 if nevents == 1 & exgr == 5
replace exgr = 0

* put data after each exposure together 
append using 1stexpo
append using 2ndexpo
append using 3rdexpo
drop exev
bysort indiv: gen ninterval = _n
save intuss_cens, replace

do intuss_cens_pseudo
estimates replay, eform

* use this command to get bootstrap confidence intervals
* bootstrap "do intuss_cens_pseudo" _b, reps(1000) cluster(bootindiv)