d <- data.frame(time=c(3,2,1,4),event = c(1,0,1,1),
z1=c(1,1,1,0),
z2=c(0,0,NA,0),
z3=c(0,NA,NA,1),
z4=c(NA,NA,NA,0))
d<-d[order(d$time),]
d
## time event z1 z2 z3 z4
## 3 1 1 1 NA NA NA
## 2 2 0 1 0 NA NA
## 1 3 1 1 0 0 NA
## 4 4 1 0 0 1 0
z <- d[3:6]
loc <- which(d$event == 1)
logPL <- function(b){
nomi <- exp(sapply(c(1,3,4),FUN = function(x) z[[x]][x])*b)
zr <- sapply(c(1,3,4),FUN = function(x) z[[x]][x:4])
denomi<-sapply(zr, function(x) sum(exp(x*b)))
log(prod(nomi/denomi))
}
plot(-4:4,sapply(-4:4,logPL),type = 'b')
optimize(function(x) -logPL(x),
c(-4,4))
## $minimum
## [1] -0.54932620658373565
##
## $objective
## [1] 2.010105077578161
U_prime <- function(b){
zr <- sapply(c(1,3,4),FUN = function(x) z[[x]][x:4])
nomi <- sapply(zr, function(x) sum(x^2*exp(x*b))*sum(exp(x*b)) -
sum(x*exp(x*b))^2 )
denomi<-sapply(zr, function(x) sum(exp(x*b))^2 )
sum(-nomi/denomi)
}
sqrt(-U_prime(-0.5493048727361014))#standard error
## [1] 0.68125003863310529
NO INTERACTION
library(survival)
## Warning: package 'survival' was built under R version 3.2.5
suppressMessages(library(dplyr))
## Warning: package 'dplyr' was built under R version 3.2.5
heart_age<-mutate(heart , age_group = ifelse(age > 0 ,1,0))
cox.strata <- coxph(Surv(start,stop,event) ~ transplant +strata(age_group),data = heart_age,method = 'breslow')
cox.strata
## Call:
## coxph(formula = Surv(start, stop, event) ~ transplant + strata(age_group),
## data = heart_age, method = "breslow")
##
## coef exp(coef) se(coef) z p
## transplant1 0.048904430381 1.050119986377 0.310580638897 0.15746 0.87488
##
## Likelihood ratio test=0.02 on 1 df, p=0.8747293573058
## n= 172, number of events= 75
the coefficient beta for transplant is 0.0489 and its coresponding wald statistic is 0.15746 which is not significant. So we can suggest that the beta of transplant is close to zero.
INTERACTION
cox.strata.interactive <- coxph(Surv(start,stop,event) ~ transplant * strata(age_group),data = heart_age,method = 'breslow')
cox.strata.interactive
## Call:
## coxph(formula = Surv(start, stop, event) ~ transplant * strata(age_group),
## data = heart_age, method = "breslow")
##
## coef exp(coef)
## transplant1 -0.16763657302 0.84566111491
## transplant1:strata(age_group)age_group=1 0.52575435061 1.69173452814
## se(coef) z p
## transplant1 0.40277477551 -0.41620 0.67726
## transplant1:strata(age_group)age_group=1 0.63523616548 0.82765 0.40787
##
## Likelihood ratio test=0.72 on 2 df, p=0.6993745561392
## n= 172, number of events= 75