In timonelmer/dena: DENA - dynamic ego-centric network analysis

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

load(file = "../m_cmm_vignette.RData")


# and other packages that we need for visualization
library(ggplot2)
library(survival)
library(coxme)
library(survminer)
require(coxed)

# load the dena package with the data and functions
library(dena)

Basic concepts of Survival analysis

survival Function (with plot)

data(simdat) # the  Rcode to simulate this dataset can be found at the bottom of this script
survivalFunction(dat = simdat,timeVar =  "time", eventVar = "event", returnDF = F, verbose = F)

survivalFunction(dat = simdat,timeVar =  "time", eventVar = "event", returnDF = F, verbose = F, plotType = "Nelson-Aalen")

Cox Proportional Hazard Model

fit a coxph model

fit <-coxph(formula = Surv(time, event) ~ X, data = simdat)
summary(fit)
coxR2(fit) # computes the pseudo r-squared (based on the cox-snell method)

## Schoenfeld test for proportional hazard assumtion ##
# should be non-significant for proportional hazard assumption to hold
cox.zph(fit) 
survminer::ggcoxzph(cox.zph(fit)) # a plot of the residuals

Frailty models

data(simdat2) # the  Rcode to simulate this dataset can be found at the bottom
# of this script


frailty_model <- coxme(Surv(time, event) ~ Covariate1 + Covariate2 + (1|id),dat = simdat2) 
summary(frailty_model)

we can visualize the fixed effects coefficients and the frailty distribution with

plot.coxme(frailty_model)
plot.coxme(frailty_model, type ="frailty")

the pseudo-R^2 (using the Cox and Snell method) can be computed with

coxR2(frailty_model) # computes the pseudo r-squared (based on the cox-snell method)

comparison of estimates with coxme, coxph, and frailtyEM package

# with the survival package
m.coxph <- coxph(Surv(time, event) ~ Covariate1 + Covariate2 + frailty(id), simdat2)

# with the frailtyEM
m.frailtyEM <- frailtyEM::emfrail(Surv(time, event) ~ Covariate1 + Covariate2 + cluster(id), simdat2)

# compare estimates
tmp <- cbind(frailty_model$coefficients[1:2],m.coxph$coefficients[1:2], m.frailtyEM$coefficients)
colnames(tmp) <- c("coxme","coxph","fraityEM")
tmp

Model Diagnostics

## proportional hazard assumption
cox.zph(m.coxph)
ggcoxzph(cox.zph(m.coxph)) # a plot of the residuals

## Testing non-linearity of continuous variables
ggcoxdiagnostics(m.coxph, type = "martingale")

## outlier observations
ggcoxdiagnostics(m.coxph, type = "deviance")

Multistate model

Now, let us take the type of interaction into account. The following histogram shows the waiting times per interaction type:

ggplot(simdat2, aes(x = time, fill = type), alpha = 0.2) + geom_histogram(position = "identity", alpha = 0.7)

with the cmm function we can estimate a multistate model with one starting state and multiple target states (variable type):

m.cmm <- dena::cmm(Surv(time, type) ~ Covariate2 + Covariate3 + (1|id), dat = simdat2)


m.coxph2 <- coxph(Surv(time, event) ~ Covariate2 + Covariate3, simdat2) # without frailty

simdat2[nrow(simdat2)+1,] <- simdat2[nrow(simdat2),]
simdat2[nrow(simdat2),"type"] <- "(censored)"
simdat2[nrow(simdat2),"event"] <- 0
simdat2$type <- factor(simdat2$type, levels =c("(censored)","family","friend","partner"))

m.survfit <- survfit(Surv(time, as.factor(type)) ~ I(Covariate3<median(Covariate3)), simdat2)
m.survfit$transitions

m.survfit # order of plotting
plot(m.survfit, col = c(1,1,2,2,3,3)+1, lty = c(1,2,1,2,1,2),
     xlab ="Minutes since last interaction", ylab = "Probability in state"
     ) # also see p.7 https://cran.r-project.org/web/packages/survival/vignettes/compete.pdf
legend(0, .6, c("> Cov3 median","< Cov3 median","family","friend","partner"),
 col = c(1,1,2,3,4), lty = c(1,2,1,1,1), lwd=1, bty='n')

m.cmm[[1]] # summary of the estimates

we can also visualize the fixed effects coefficients

plot.cmm(m.cmm[[1]])

we see that Covariate3 is highly associated with sooner interactions with partners. From the following plot, that we also see descriptively that for interactions with partners Covariate3 is highly correlated with the waiting time.

# ggplot(simdat2, aes(x = Covariate3, fill = type), alpha = 0.2) + geom_histogram(position = "identity", alpha = 0.7)

psych::pairs.panels(simdat2[simdat2$type %in% "partner",c("time","Covariate3")])

other details of the results can be accessed too

m.cmm[[2]] # list of coxme fit for each type
m.cmm[[3]] # pseudo r^2
m.cmm[[4]] # frailty estimates

m.cmm <- cmm(Surv(time, type) ~ Covariate2 + Covariate3 + frailty(id), dat = simdat2, verbose = T)  # alternative writing of frailty term

Model Diagnostics

# proportional hazard assumption
ggcoxzph(cox.zph(m.cmm[[2]][[1]]), caption ="Submodel: alone to family") # for the alone to family model; Figure in the Supplementary Materials
ggcoxzph(cox.zph(m.cmm[[2]][[2]]), caption ="Submodel: alone to friend") # for the alone to friend model
ggcoxzph(cox.zph(m.cmm[[2]][[3]]), caption ="Submodel: alone to partner") # for the alone to partner model

## Testing non-linearity of continuous variables
ggcoxdiagnostics(m.cmm[[2]][[1]], type = "martingale", title ="Submodel: alone to family") # Figure in the Supplementary Materials
ggcoxdiagnostics(m.cmm[[2]][[2]], type = "martingale", title ="Submodel: alone to frined")
ggcoxdiagnostics(m.cmm[[2]][[3]], type = "martingale", title ="Submodel: alone to partner")

## outlier observations (Xue and Schifano, 2017 use +- 1.96 as outlier indication)
ggcoxdiagnostics(m.cmm[[2]][[1]], type = "deviance", title ="Submodel: alone to family")#  Figure in the Supplementary Materials
ggcoxdiagnostics(m.cmm[[2]][[2]], type = "deviance", title ="Submodel: alone to frined")
ggcoxdiagnostics(m.cmm[[2]][[3]], type = "deviance", title ="Submodel: alone to partner")

Simulation of the simdata dataset

set.seed(44)
dat.sim <- coxed::sim.survdata(N=200, T=500, num.data.frames=1, xvars = 1, 
                        mu = 0.5,
                        sd = 0.1,
                        censor = 0,
                        beta=c(2))

dat.sim$data$X <- round(dat.sim$data$X)
simdat <- dat.sim$data[,c(2,1,3)]
colnames(simdat) <- c("time","X","event")
simdat <- simdat[sample(1:nrow(simdat),nrow(simdat)),]
simdat$ID <- 1:nrow(simdat)
simdat <- simdat[,c(4,1,3,2)]
rownames(simdat) <- simdat$ID
simdat$event <- simdat$event*1
head(simdat)
#save(simdat, file = "data/simdat.RData")

Simulation of the simdata2 dataset

set.seed(133)
require(frailtySurv)
simdat2 <- genfrail(N = 30, K = 20, 
                                 beta = c(log(2), -log(1)),
                                 frailty = "gamma", theta = 2, 
                                 censor.rate = 0,
                                 lambda_0=function(t, tau=4.6, C=0.01) (tau*(C*t)^tau)/t)

colnames(simdat2) <- c("id","rep","time","event","Covariate1","Covariate2")
simdat2$type <- factor(sample(c("partner","friend","family"), nrow(simdat2), replace = T))

# create fake covariate 3 that is associated with eventtype and time 
simdat2$Covariate3 <- ifelse(simdat2$type == "partner" & simdat2$time < 200,1,0) + 
  rnorm(nrow(simdat2),mean = 0.2, sd = 0.2)

head(simdat2)
#save(simdat2, file = "data/simdat2.RData")