# Analysis of MMR vaccination and hospital admission for ITP

############################################################
# Data input and variable definition
############################################################

# Columns are: caseno, age at event, age at start of observation - 1, age at end of observation, age at MMR
# All ages are in days. 
# Note that time intervals are all of the form (a,b]. Since time is discrete here, (a, b] = [a+1, b].

dat<-matrix(c(
1,	691,	453,	730,	670,
2,	722,	365,	730,	868,
3,	442,	365,	730,	540,
4,	429,	365,	730,	378,
5,	414,	365,	730,	710,
5,	418,	365,	730,	710,
6,	708,	438,	730,	487,
7,	615,	365,	730,	461,
8,	463,	365,	730,	526,
9,	440,	365,	730,	529,
9,	473,	365,	730,	529,
10,	477,	365,	730,	458,
11,	396,	365,	730,	374,
12,	676,	365,	730,	428,
13,	480,	365,	730,	446,
14,	633,	365,	730,	423,
15,	403,	365,	730,	365,
16,	419,	365,	730,	369,
16,	443,	365,	730,	369,
17,	553,	365,	730,	889,
17,	666,	365,	730,	889,
18,	705,	365,	730,	389,
19,	419,	365,	730,	389,
20,	402,	365,	730,	385,
21,	406,	365,	730,	458,
22,	494,	365,	730,	468,
23,	374,	365,	730,	819,
23,	389,	365,	730,	819,
23,	452,	365,	730,	819,
23,	522,	365,	730,	819,
23,	564,	365,	730,	819,
24,	598,	365,	730,	430,
25,	409,	365,	730,	384,
26,	612,	365,	730,	398,
27,	381,	365,	723,	427,
28,	438,	365,	730,	427,
29,	425,	365,	677,	647,
30,	543,	365,	677,	422,
31,	609,	365,	674,	860,
32,	412,	365,	730,	387,
33,	407,	365,	730,	396,
34,	484,	365,	730,	408,
34,	623,	365,	730,	408,
35,	411,	365,	730,	383), ncol=5, byrow=T)

case <-dat[,1]
evday<-dat[,2]
start<-dat[,3]
end  <-dat[,4]
mmr  <-dat[,5]

############################################################
# Baseline parametric model
############################################################

# Define risk periods and age boundaries

ex1 <- mmr - 1
ex2 <- mmr + 14
ex3 <- mmr + 28
ex4 <- mmr + 42
expo <- cbind(ex1, ex2, ex3, ex4)
expolev <- c(1, 2, 3, 0)

age <- c(426, 487, 548, 609, 670)

# Expand data

ncuts <- ncol(expo) + length(age) + 2
nevents <- length(evday)
ind <- rep(1:nevents, times = ncuts)
cutp <- c(as.matrix(start), as.matrix(end), expo, rep(age, each = nevents))
o <- order(ind, cutp)
ind = as.factor(ind[o])
cutp = cutp[o]

interval <- c(0, cutp[2:length(ind)]-cutp[1:length(ind)-1])
interval <- ifelse(cutp<=start[ind], 0, interval)
interval <- ifelse(cutp>end[ind], 0, interval)

event <- ifelse(evday[ind]>cutp-interval, 1, 0)
event <- ifelse(evday[ind]<=cutp, event, 0)

agegr <- cut(cutp, breaks = c(min(start), age, max(end)), labels = FALSE)

exgr <- rep(0, nevents*ncuts)
for(i in 1:ncol(expo)){
exgr <- ifelse(cutp > expo[,i][ind], expolev[i], exgr)
}

# Fit model

chopdat <- data.frame(indiv = ind[interval!=0], event = event[interval!=0], interval = interval[interval!=0], agegr = agegr[interval!=0], exgr = exgr[interval!=0], loginterval = log(interval[interval!=0]))

require(survival)
mod1 <- clogit(event ~ factor(exgr) + factor(agegr) + strata(indiv) + offset(loginterval), data = chopdat)
summary(mod1)

# Fit baseline model without age (likelihood ratio test for age effect is based on summaries of mod1 and mod1a)

mod1a <- clogit(event ~ factor(exgr) + strata(indiv) + offset(loginterval), data = chopdat)
summary(mod1a)

# Combine risk periods into a single 0 - 42 day risk period and refit baseline model

exgr1 <- ifelse(exgr>0,1,0)
chopdat1 <- data.frame(indiv = ind[interval!=0], event = event[interval!=0], interval = interval[interval!=0], agegr = agegr[interval!=0], exgr1 = exgr1[interval!=0], loginterval = log(interval[interval!=0]))

mod1b <- clogit(event ~ factor(exgr1) + factor(agegr) + strata(indiv) + offset(loginterval), data = chopdat1)
summary(mod1b)


############################################################
# Model for first events only
############################################################

# Select data for first events

first<-c(1,diff(case))
evday1<-evday[first==1]
start1<-start[first==1]
end1  <-end[first==1]
mmr1  <-mmr[first==1]

# Apply model to these data

ex1 <- mmr1 - 1
ex2 <- mmr1 + 14
ex3 <- mmr1 + 28
ex4 <- mmr1 + 42
expo <- cbind(ex1, ex2, ex3, ex4)
expolev <- c(1, 2, 3, 0)

age <- c(426, 487, 548, 609, 670)

ncuts <- ncol(expo) + length(age) + 2
nevents <- length(evday1)
ind <- rep(1:nevents, times = ncuts)
cutp <- c(as.matrix(start1), as.matrix(end1), expo, rep(age, each = nevents))
o <- order(ind, cutp)
ind = as.factor(ind[o])
cutp = cutp[o]

interval <- c(0, cutp[2:length(ind)]-cutp[1:length(ind)-1])
interval <- ifelse(cutp<=start1[ind], 0, interval)
interval <- ifelse(cutp>end1[ind], 0, interval)

event <- ifelse(evday1[ind]>cutp-interval, 1, 0)
event <- ifelse(evday1[ind]<=cutp, event, 0)

agegr <- cut(cutp, breaks = c(min(start1), age, max(end1)), labels = FALSE)

exgr <- rep(0, nevents*ncuts)
for(i in 1:ncol(expo)){
exgr <- ifelse(cutp > expo[,i][ind], expolev[i], exgr)
}

# Fit model for first events

chopdat2 <- data.frame(indiv = ind[interval!=0], event = event[interval!=0], interval = interval[interval!=0], agegr = agegr[interval!=0], exgr = exgr[interval!=0], loginterval = log(interval[interval!=0]))

require(survival)
mod2 <- clogit(event ~ factor(exgr) + factor(agegr) + strata(indiv) + offset(loginterval), data = chopdat2)
summary(mod2)

############################################################
# Model with pre-MMR 'risk' period
############################################################

# Define new risk period for period 1 to 61 days pre-MMR inclusive

ex0 <- mmr - 62
ex1 <- mmr - 1
ex2 <- mmr + 14
ex3 <- mmr + 28
ex4 <- mmr + 42
expo <- cbind(ex0, ex1, ex2, ex3, ex4)
expolev <- c(1, 2, 3, 4, 0)

# Expand data and run model

age <- c(426, 487, 548, 609, 670)

ncuts <- ncol(expo) + length(age) + 2
nevents <- length(evday)
ind <- rep(1:nevents, times = ncuts)
cutp <- c(as.matrix(start), as.matrix(end), expo, rep(age, each = nevents))
o <- order(ind, cutp)
ind = as.factor(ind[o])
cutp = cutp[o]

interval <- c(0, cutp[2:length(ind)]-cutp[1:length(ind)-1])
interval <- ifelse(cutp<=start[ind], 0, interval)
interval <- ifelse(cutp>end[ind], 0, interval)

event <- ifelse(evday[ind]>cutp-interval, 1, 0)
event <- ifelse(evday[ind]<=cutp, event, 0)

agegr <- cut(cutp, breaks = c(min(start), age, max(end)), labels = FALSE)

exgr <- rep(0, nevents*ncuts)
for(i in 1:ncol(expo)){
exgr <- ifelse(cutp > expo[,i][ind], expolev[i], exgr)
}

# Fit model

chopdat3 <- data.frame(indiv = ind[interval!=0], event = event[interval!=0], interval = interval[interval!=0], agegr = agegr[interval!=0], exgr = exgr[interval!=0], loginterval = log(interval[interval!=0]))

require(survival)
mod3 <- clogit(event ~ factor(exgr) + factor(agegr) + strata(indiv) + offset(loginterval), data = chopdat3)
summary(mod3)

############################################################
# Semiparametric model
############################################################

# Apply semiparametricmodel to these data

ex1 <- mmr - 1
ex2 <- mmr + 14
ex3 <- mmr + 28
ex4 <- mmr + 42
expo <- cbind(ex1[!duplicated(case)], ex2[!duplicated(case)], ex3[!duplicated(case)], ex4[!duplicated(case)])
expolev <- c(1, 2, 3, 0)

nindivs <- length(unique(case)) 
ncases <-nrow(dat) 

eventdays <- sort(unique(evday)) 
n <- length(eventdays)  
eventdays <- rep(eventdays, times=nindivs) 
indiv <- rep(unique(case), each=n)
events <- as.vector(table(evday, case))

keep <- start[!duplicated(case)][indiv]>=eventdays|end[!duplicated(case)][indiv]<eventdays
indiv <- indiv[keep==0]
eventdays <- eventdays[keep==0]
events <- events[keep==0]

exgr <- rep(0, times = length(indiv))
for(i in 1:ncol(expo)){
exgr <- ifelse(eventdays > expo[,i][indiv], expolev[i], exgr)
}

# Fit model

semipdata <- data.frame(indiv = indiv, day = eventdays, exgr = exgr, events = events)

require(survival)
spmodel <- clogit(events ~ factor(exgr) + factor(day) + strata(indiv), data = semipdata)
summary(spmodel)

############################################################
# End
############################################################