Merlin in R

data(heart.valve, package = “joineRML”) heart.valve <- heart.valve[!is.na(heart.valve\(grad), ] print(heart.valve[heart.valve\)num %in% c(1:2, 13),c(1:3, 25, 5:6, 4, 8, 10)],row.names = FALSE)

heart.valve\(id <- heart.valve\)num heart.valve\(stime <- heart.valve\)fuyrs heart.valve\(stime[duplicated(heart.valve\)id)] <- NA heart.valve\(died <- heart.valve\)status heart.valve\(died[duplicated(heart.valve\)id)] <- NA heart.valve\(type <- as.numeric(heart.valve\)hs) - 1 print(heart.valve[heart.valve$num %in% c(1:2, 13),c(26, 2:3, 29, 27:28, 8, 10, 4)],row.names = FALSE)

library(merlin) m1<-merlin(model = log.grad ~ sex + age + time, family = “gaussian”, data = heart.valve) summary(m1)

heart.valve\(m1res <- heart.valve\)log.grad - predict(m1, stat = “mu”) library(ggplot2) ggplot(heart.valve, aes(x = time, y = m1res)) + geom_point() + geom_hline(yintercept = 0, color = “blue”) + xlab(“Time (years)”) + ylab(“Residual”) + theme_classic() m2<-merlin(model = log.grad ~ sex + age + rcs(time, df=3, orthog = TRUE), timevar = “time”, family = “gaussian”, data = heart.valve) summary(m2)

heart.valve\(m2res <- heart.valve\)log.grad - predict(m2, stat = “mu”) ggplot(heart.valve, aes(x = time, y = m2res)) + geom_point() + geom_hline(yintercept = 0, color = “blue”) + xlab(“Time (years)”) + ylab(“Residual”) + theme_classic() m3<-merlin(model = log.grad ~ sex + age + rcs(time, df=3, orthog = TRUE) + M1[id]1 + time:M2[id]1, timevar = “time”, levels = “id”, covariance = “unstructured”, family = “gaussian”, data = heart.valve) summary(m3)

ldata <- heart.valve[!is.na(heart.valve\(log.grad), ] ldata\)pred1 <- predict(m3, stat = “mu”, predmodel = 1, type = “marginal”) print(ldata[ldata$id %in% c(1:2, 13), c(“id”, “time”, “log.grad”, “pred1”)],row.names = FALSE)

logl_gaussian <- function(gml) {
y <- merlin_util_depvar(gml) xzb <- merlin_util_xzb(gml) se <- exp(merlin_util_ap(gml, 1))

mu <- (sweep(xzb, 1, y, “-”))^2 logl <- ((-0.5 * log(2 * pi) - log(se)) - (mu / (2 * se^2))) return(logl) }

m4 <- merlin(log.grad ~ sex + age + time + ap(1),family=“user”,userf=“logl_gaussian”,data=heart.valve) summary(m4) m5 <-merlin(model= Surv(stime,died)~age+type,family=“weibull”,data=heart.valve) summary(m5)

sdata<-heart.valve[!is.na(heart.valve\(stime),] sdata\)pred2<-predict(m5,stat=“survival”,premodel=1,type=“fixedonly”) print(sdata[sdata$id %in% c(1:2,13),c(“id”,“stime”,“died”,“age”,“type”,“pred2”)],row.names=FALSE)

p_50 <- predict(m5, stat = “survival”, type = “fixedonly”, at = c(age = 50, type = 1)) p_60 <- predict(m5, stat = “survival”, type = “fixedonly”, at = c(age = 60, type = 1)) p_70 <- predict(m5, stat = “survival”, type = “fixedonly”, at = c(age = 70, type = 1)) surv_pred<-data.frame(p_50,p_60,p_70, stime=heart.valve[!is.na(heart.valve\(stime),"stime"]) surv_pred <- surv_pred[!duplicated(surv_pred\)stime),]

library(reshape2) surv_pred <- melt(surv_pred, id.var = “stime”) ggplot(surv_pred, aes(x = stime, y = value, linetype = variable)) + geom_line(size = 0.6) + coord_cartesian(ylim = c(0, 1)) + xlab(“Time (years)”) + ylab(“Survival probability”) + theme_classic() + theme(legend.position = c(0.2, 0.2)) + scale_linetype_discrete(name = “Age”, labels = c(“50 years”, “60 years”, “70 years”))

m6 <- merlin(model = Surv(stime, died) ~ age + type + rcs(stime, df = 3, log = TRUE, event = TRUE), timevar = “stime”, family = “rp”, data = heart.valve) summary(m6)

base_m5 <- predict(m5, stat = “hazard”, type =“fixedonly”,at = c(age = 0, type = 0)) base_m6 <- predict(m6, stat = “hazard”, type =“fixedonly”,at = c(age = 0, type = 0)) base_pred <- data.frame(stime = heart.valve[!is.na(heart.valve\(stime), "stime"],base_m5,base_m6) base_pred <- base_pred[!duplicated(base_pred\)stime),] base_pred <- melt(base_pred, id.var = “stime”)

ggplot(base_pred, aes(x = stime, y = value, linetype = variable)) + geom_line(size = 0.6) + xlab(“Time (years)”) + ylab(“Baseline hazard”) + theme_classic() + theme(legend.position = c(0.2, 0.8)) + scale_linetype_discrete(name = “Model”, labels = c(“Weibull”, “RP - 3df”)) m7 <- merlin(model = Surv(stime, died) ~ age + type + type:fp(stime,powers=c(0))+ rcs(stime, df=3, log=TRUE,event=TRUE),timevar = “stime”,family = “rp”, data = heart.valve) summary(m7)

m8 <- merlin(model = Surv(stime, died) ~ type + fp(age,powers=c(1,1)) + rcs(stime, df=3, log=TRUE,event=TRUE),timevar = “stime”,family = “rp”, data = heart.valve) summary(m8)

mlsurv(formula=Surv(stime,died)~age + type, distribution=“weibull”, data=heart.valve)

set.seed(6342) heart.valve\(cardio[!is.na(heart.valve\)died)] <-rbinom(length(heart.valve\(died[!is.na(heart.valve\)died)]), 1, 0.6) heart.valve\(other <- 1 - heart.valve\)cardio heart.valve\(cardio[heart.valve\)died == 0 & !is.na(heart.valve\(died)] <- 0 heart.valve\)other[heart.valve\(died == 0 & !is.na(heart.valve\)died)] <- 0

m9 <- merlin(model = list(Surv(stime, cardio) ~ type + rcs(stime, df = 3, log = TRUE, event = TRUE), Surv(stime, other) ~ type + rcs(stime, df = 3, log = TRUE, event = TRUE)),timevar = c(“stime”, “stime”),family = c(“rp”, “rp”),data = heart.valve) summary(m9)

card_homo <- predict(m9, stat = “cif”, type = “fixedonly”,predmodel = 1, at = c(age = 50, type = 0)) card_stent <- predict(m9, stat = “cif”, type = “fixedonly”,predmodel = 1, at = c(age = 50, type = 1)) other_homo <- predict(m9, stat = “cif”, type = “fixedonly”,predmodel = 2, at = c(age = 50, type = 0)) other_stent <- predict(m9, stat = “cif”, type = “fixedonly”,predmodel = 2, at = c(age = 50, type = 1)) stime = rep(heart.valve\(stime[!is.na(heart.valve\)stime)], 4) pred = c(card_homo, card_stent, other_homo, other_stent) valve = rep(c(rep(“Homograft valve”, length(card_homo)),rep(“Stentless valve”,length(card_stent))),2) cod = c(rep(“cardio”, length(card_homo)2),rep(“other”,length(card_stent)2)) pred_comp <- data.frame(stime, pred, valve, cod) pred_comp <- pred_comp[!duplicated(pred_comp[, c(1, 3, 4)]),]

ggplot(pred_comp, aes(x = stime, y = pred, fill = cod)) + geom_area() + xlab(“Time (years)”) + ylab(“Cumulative Incidence”) + theme_classic(base_size = 10) + theme(legend.position = c(0.19, 0.83)) + scale_fill_discrete(name = “Cause of death”, labels = c(“Other”, “Cardiovascular disease”), guide = guide_legend(reverse = TRUE)) + facet_grid(. ~ valve)

m10 <- merlin(model = list(Surv(stime, died) ~ type + M1[id], log.grad ~ sex + age + time + M1[id] * 1), timevar = c(“stime”, “time”), levels = c(“id”), family = c(“weibull”, “gaussian”), data = heart.valve) summary(m10)

m11 <- merlin(model = list(Surv(stime, died) ~ type + EV[log.grad],log.grad ~ sex + age + time + M1[id] * 1),timevar = c(“stime”, “time”),levels = c(“id”),family = c(“weibull”, “gaussian”),data = heart.valve) summary(m11)

m12 <- merlin(model = list(Surv(stime, died) ~ type + dEV[log.grad],log.grad ~ sex + age + time + time:M1[id] * 1),timevar = c(“stime”, “time”),levels = c(“id”),family = c(“weibull”, “gaussian”),data = heart.valve) summary(m12)

m13 <- merlin(model=list(Surv(stime,died)~type + EV[log.grad] + EV[log.grad]:fp(stime, powers = c(0)),log.grad ~ time + M1[id]*1),timevar =c(“stime”, “time”),levels = c(“id”),family = c(“weibull”,“gaussian”),data = heart.valve) summary(m13)

heart.valve\(catef <- 0 heart.valve\)catef[heart.valve$ef > 70] <- 1 m14 <- merlin(model = list(Surv(stime, cardio) ~ type + EV[log.grad] + M2[id],Surv(stime, other) ~ type + type:fp(stime, powers = c(0)) + M1[id], log.grad ~ age + type + rcs(time, df = 3, orthog = TRUE) + M1[id] * 1, catef ~ fp(time, powers = c(1)) + M2[id] * 1 ), timevar = c(“stime”, “stime”, “time”, “time”), levels = c(“id”), covariance = “unstructured”, family = c(“weibull”, “weibull”, “gaussian”, “bernoulli”), data = heart.valve, control = list(ip = 9)) summary(m14)