install.packages("Bayes")
library(lme4)
?glmer
?title
setwd("~/Dropbox/Analysen/Aalen_analysis/pmr2")
#--------------------------------------------------------------------------
# Run Aalen's additive hazard model
# Revealing the time course of signals influencing the generation of
# secondary saccades using Aalen's additive hazards model
# Ohl, Sven & Kliegl, Reinhold
# Vision Research 2016
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# Comments
# error_hor2 <- starting position of secondary saccade
# error_hor3 <- ending positioin of secondary saccade
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# Clear workspace
rm(list=ls())
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# Define pathes
path_wd 	 <- getwd()
path_data 	 <- paste(path_wd,"/_datfiles/",sep="")
path_figures <- paste(path_wd,"/_Figures/",sep="")
path_script  <- paste(path_wd,"/_RScripts/",sep="")
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# load packages
source(paste(path_script,"loadPackages.R",sep=""))
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# define colors
source(paste(path_script,"loadColors.R",sep=""))
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# load data and provide colnames
source(paste(path_script,"loadData.R",sep=""))
#--------------------------------------------------------------------------
str(drev)
source(paste(path_script,"genVariables.R",sep=""))
str(drev)
#--------------------------------------------------------------------------
# Run Aalen's additive hazard model
# Revealing the time course of signals influencing the generation of
# secondary saccades using Aalen's additive hazards model
# Ohl, Sven & Kliegl, Reinhold
# Vision Research 2016
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# Comments
# error_hor2 <- starting position of secondary saccade
# error_hor3 <- ending positioin of secondary saccade
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# Clear workspace
rm(list=ls())
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# Define pathes
path_wd 	 <- getwd()
path_data 	 <- paste(path_wd,"/_datfiles/",sep="")
path_figures <- paste(path_wd,"/_Figures/",sep="")
path_script  <- paste(path_wd,"/_RScripts/",sep="")
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# load packages
source(paste(path_script,"loadPackages.R",sep=""))
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# define colors
source(paste(path_script,"loadColors.R",sep=""))
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# load data and provide colnames
source(paste(path_script,"loadData.R",sep=""))
#--------------------------------------------------------------------------
#--------------------------------------------------------------------------
# generate new variables
source(paste(path_script,"genVariables.R",sep=""))
#--------------------------------------------------------------------------
str(drev)
source(paste(path_script,"loadTheme.R",sep=""))
drev$event   <- ifelse(drev$latency<1100,1,0)
table(drev$event)
#---------------------------------------------------------------
# Survival Analysis
#
# see also help for hidden functions:
#?Csmooth2B
#?CsmoothB
#---------------------------------------------------------------
#---------------------------------------------------------------
# Compute Aalen's additive hazard model
fit_aalen3 <- aalen(Surv(latency, event) ~ 1 + ecccond*underovercond + error_hor_abs, data = drev, max.time=1100, cluster=drev$vp, residuals=1, n.sim=10000)
# Plot cumulative coefficients with confidence intervals
par(mfrow=c(1,5),bty="l")
plot(fit_aalen3)
#---------------------------------------------------------------
#---------------------------------------------------------------
# Compute corresponding residuals for continuous covariate
resid_aalen3 <- cum.residuals(fit_aalen3, data=drev,cum.resid=1,n.sim=10000)
# Plot residuals along with 50 simulations
plot(resid_aalen3,score=2)
#---------------------------------------------------------------
#---------------------------------------------------------------
# Plot cumulative coefficients along with smoothed derivatives
# get the data from the aalen model
source(paste(path_script,"plot.cums.sven.R",sep=""))
par(mfrow=c(2,3))
gginfo <- plot.cums.sven(fit_aalen3)
B     <- gginfo[[1]]
ul    <- gginfo[[2]]
nl    <- gginfo[[3]]
est   <- gginfo[[4]]
namen <- dimnames(B)[[2]]
namen <- namen[2:length(namen)]
# construct data.frame
df <- data.frame(
Latency = rep(B[,1],time=c(ncol(B)-1)),
varname = rep(namen,each=nrow(B)),
varest  = as.vector(est),
upper   = as.vector(ul),
lower   = as.vector(nl)
)
# assign name to data.frame
levels(df$varname) <- c("Intercept","ecc","eccxoverunder","error","underover")
df$varname <- as.character(df$varname)
df$varname <- factor(df$varname,levels=c("Intercept","ecc","eccxoverunder","error","underover"))
# identify the different components
idxinter         <- which(df$varname=="Intercept")
idxecc           <- which(df$varname=="ecc")
idxoverunderxecc <- which(df$varname=="eccxoverunder")
idxerror         <- which(df$varname=="error")
idxunderover     <- which(df$varname=="underover")
#-----------------------------------------------------
#-----------------------------------------------------
#0.)
spl <- smooth.spline(df$Latency[idxinter], df$varest[idxinter], df=12)
spl0 <- spl
xx <- unique(sort(c(seq(0,1100, by = 1), kn <- unique(df$Latency[idxinter]))))
i.kn <- match(kn, xx)# indices of knots within xx
pp <- predict(spl, xx, deriv = 1)
idx <- which(pp$y[i.kn]==max(pp$y[i.kn]))
idxx<- kn[idx]
df6x <- data.frame(Latency=pp$x,deriv=pp$y,varname="intercept",maxpoint=idxx,maxderiv=pp$y[idx])
#1)
spl <- smooth.spline(df$Latency[idxecc], df$varest[idxecc], df=12)
spl1 <- spl
xx <- unique(sort(c(seq(0,1100, by = 1), kn <- unique(df$Latency[idxecc]))))
i.kn <- match(kn, xx)# indices of knots within xx
pp <- predict(spl, xx, deriv = 1)
idx <- which(pp$y[i.kn]==max(pp$y[i.kn]))
idxx<- kn[idx]
df6a <- data.frame(Latency=pp$x,deriv=pp$y,varname="eccentricity",maxpoint=idxx,maxderiv=pp$y[idx])
#2)
spl <- smooth.spline(df$Latency[idxoverunderxecc], df$varest[idxoverunderxecc], df=12)
spl2 <- spl
xx <- unique(sort(c(seq(0,800, by = 1), kn <- unique(df$Latency[idxoverunderxecc]))))
i.kn <- match(kn, xx)# indices of knots within xx
pp <- predict(spl, xx, deriv = 1)
idx <- which(pp$y[i.kn]==min(pp$y[i.kn]))
idxx<- kn[idx]
df6b <- data.frame(Latency=pp$x,deriv=pp$y,varname="underoverxecc",maxpoint=idxx,maxderiv=pp$y[idx])
#3)
spl <- smooth.spline(df$Latency[idxerror], df$varest[idxerror], df=12)
spl3 <- spl
xx <- unique(sort(c(seq(0,1100, by = 1), kn <- unique(df$Latency[idxerror]))))
i.kn <- match(kn, xx)# indices of knots within xx
pp <- predict(spl, xx, deriv = 1)
idx <- which(pp$y[i.kn]==max(pp$y[i.kn]))
idxx<- kn[idx]
df6c <- data.frame(Latency=pp$x,deriv=pp$y,varname="error",maxpoint=idxx,maxderiv=pp$y[idx])
#4)
spl <- smooth.spline(df$Latency[idxunderover], df$varest[idxunderover], df=12)
spl4 <- spl
xx <- unique(sort(c(seq(0,1100, by = 1), kn <- unique(df$Latency[idxunderover]))))
i.kn <- match(kn, xx)# indices of knots within xx
pp <- predict(spl, xx, deriv = 1)
idx <- which(pp$y[i.kn]==min(pp$y[i.kn]))
idxx<- kn[idx]
df6d <- data.frame(Latency=pp$x,deriv=pp$y,varname="underover",maxpoint=idxx,maxderiv=pp$y[idx])
# put all data.frames together
df6 <- rbind(df6x,df6a,df6b,df6c,df6d)
levels(df6$varname) <- c("Intercept","eccentricity (0=close; 1=distant)","under-/overshoot x ecc","error","under-/overshoot")
spl <- data.frame(Latency = rep(spl0$x,times=5),
spline  = c(spl0$y,spl1$y,spl2$y,spl3$y,spl4$y),
varname = rep(c("Intercept","eccentricity (0=close; 1=distant)","under-/overshoot x ecc","error","under-/overshoot"),each=length(spl0$x)
))
levels(df$varname) <- c("Intercept","eccentricity (0=close; 1=distant)","under-/overshoot x ecc","error","under-/overshoot")
df6$varname <- as.character(df6$varname)
df6$varname <- factor(df6$varname,levels=c("Intercept","error","eccentricity (0=close; 1=distant)","under-/overshoot","under-/overshoot x ecc"))
df$varname <- as.character(df$varname)
df$varname <- factor(df$varname,levels=c("Intercept","error","eccentricity (0=close; 1=distant)","under-/overshoot","under-/overshoot x ecc"))
spl$varname <- as.character(spl$varname)
spl$varname <- factor(spl$varname,levels=c("Intercept","error","eccentricity (0=close; 1=distant)","under-/overshoot","under-/overshoot x ecc"))
#-----------------------------------------------------
#Baseline
par(mfrow=c(2,5),mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(fit_aalen3$cum[,1],fit_aalen3$cum[,2],type="l",lwd=1,xlim=c(0,1200),ylim=c(-1.25,2),ylab="",xlab="Time",bty="n",tck=-0.04)
abline(h=0)
#Error
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(fit_aalen3$cum[,1],fit_aalen3$cum[,5],type="l",lwd=1,xlim=c(0,1200),ylim=c(-1.25,2),ylab="",yaxt="n",xlab="Time",bty="n",tck=-0.04)
abline(h=0)
#Underovershoot
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(fit_aalen3$cum[,1],fit_aalen3$cum[,4],type="l",lwd=1,xlim=c(0,1200),ylim=c(-1.25,2),ylab="",yaxt="n",xlab="Time",bty="n",tck=-0.04)
abline(h=0)
#Eccentricity
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(fit_aalen3$cum[,1],fit_aalen3$cum[,3],type="l",lwd=1,xlim=c(0,1200),ylim=c(-1.25,2),ylab="",yaxt="n",xlab="Time",bty="n",tck=-0.04)
abline(h=0)
#Underovershoot x Eccentricity
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(fit_aalen3$cum[,1],fit_aalen3$cum[,6],type="l",lwd=1,xlim=c(0,1200),ylim=c(-1.25,2),ylab="",yaxt="n",xlab="Time",bty="n",tck=-0.04)
abline(h=0)
#Baseline
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(df6x$Latency,df6x$deriv,xlim=c(0,1200),ylim=c(-0.01,0.015),type="l",lwd=1,ylab="",xlab="Time",bty="n",tck=-0.04)
lines(c(df6x$maxpoint[1],df6x$maxpoint[1]),c(df6x$maxderiv[1],-0.01),lty=2,col=defgrey)
abline(h=0)
#Error
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(df6c$Latency,df6c$deriv,xlim=c(0,1200),ylim=c(-0.01,0.015),type="l",lwd=1,ylab="",xlab="",yaxt="n",bty="n",tck=-0.04)
lines(c(df6c$maxpoint[1],df6c$maxpoint[1]),c(df6c$maxderiv[1],-0.01),lty=2,col=defgrey)
abline(h=0)
#Underovershoot
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(df6d$Latency,df6d$deriv,xlim=c(0,1200),ylim=c(-0.01,0.015),type="l",lwd=1,ylab="",xlab="",yaxt="n",bty="n",tck=-0.04)
lines(c(df6d$maxpoint[1],df6d$maxpoint[1]),c(df6d$maxderiv[1],-0.01),lty=2,col=defgrey)
abline(h=0)
#Eccentricity
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(df6a$Latency,df6a$deriv,xlim=c(0,1200),ylim=c(-0.01,0.015),type="l",lwd=1,ylab="",xlab="",yaxt="n",bty="n",tck=-0.04)
lines(c(df6a$maxpoint[1],df6a$maxpoint[1]),c(df6a$maxderiv[1],-0.01),lty=2,col=defgrey)
abline(h=0)
#Underovershoot x Eccentricity
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(df6b$Latency,df6b$deriv,xlim=c(0,1200),ylim=c(-0.01,0.015),type="l",lwd=1,ylab="",xlab="",yaxt="n",bty="n",tck=-0.04)
lines(c(df6b$maxpoint[1],df6b$maxpoint[1]),c(df6b$maxderiv[1],-0.01),lty=2,col=defgrey)
abline(h=0)
#-----------------------------------------------------
#---------------------------------------------------------------
# Survival Analysis
#
# see also help for hidden functions:
#?Csmooth2B
#?CsmoothB
#---------------------------------------------------------------
#---------------------------------------------------------------
# Compute Aalen's additive hazard model
fit_aalen3 <- aalen(Surv(latency, event) ~ 1 + ecccond*underovercond + error_hor_abs, data = drev, max.time=1100, cluster=drev$vp, residuals=1, n.sim=10000)
# Plot cumulative coefficients with confidence intervals
par(mfrow=c(1,5),bty="l")
plot(fit_aalen3)
#---------------------------------------------------------------
#---------------------------------------------------------------
# Compute corresponding residuals for continuous covariate
resid_aalen3 <- cum.residuals(fit_aalen3, data=drev,cum.resid=1,n.sim=10000)
# Plot residuals along with 50 simulations
plot(resid_aalen3,score=2)
#---------------------------------------------------------------
#---------------------------------------------------------------
# Plot cumulative coefficients along with smoothed derivatives
# get the data from the aalen model
source(paste(path_script,"plot.cums.sven.R",sep=""))
par(mfrow=c(2,3))
gginfo <- plot.cums.sven(fit_aalen3)
B     <- gginfo[[1]]
ul    <- gginfo[[2]]
nl    <- gginfo[[3]]
est   <- gginfo[[4]]
namen <- dimnames(B)[[2]]
namen <- namen[2:length(namen)]
# construct data.frame
df <- data.frame(
Latency = rep(B[,1],time=c(ncol(B)-1)),
varname = rep(namen,each=nrow(B)),
varest  = as.vector(est),
upper   = as.vector(ul),
lower   = as.vector(nl)
)
# assign name to data.frame
levels(df$varname) <- c("Intercept","ecc","eccxoverunder","error","underover")
df$varname <- as.character(df$varname)
df$varname <- factor(df$varname,levels=c("Intercept","ecc","eccxoverunder","error","underover"))
# identify the different components
idxinter         <- which(df$varname=="Intercept")
idxecc           <- which(df$varname=="ecc")
idxoverunderxecc <- which(df$varname=="eccxoverunder")
idxerror         <- which(df$varname=="error")
idxunderover     <- which(df$varname=="underover")
#-----------------------------------------------------
#-----------------------------------------------------
#0.)
spl <- smooth.spline(df$Latency[idxinter], df$varest[idxinter], df=12)
spl0 <- spl
xx <- unique(sort(c(seq(0,1100, by = 1), kn <- unique(df$Latency[idxinter]))))
i.kn <- match(kn, xx)# indices of knots within xx
pp <- predict(spl, xx, deriv = 1)
idx <- which(pp$y[i.kn]==max(pp$y[i.kn]))
idxx<- kn[idx]
df6x <- data.frame(Latency=pp$x,deriv=pp$y,varname="intercept",maxpoint=idxx,maxderiv=pp$y[idx])
#1)
spl <- smooth.spline(df$Latency[idxecc], df$varest[idxecc], df=12)
spl1 <- spl
xx <- unique(sort(c(seq(0,1100, by = 1), kn <- unique(df$Latency[idxecc]))))
i.kn <- match(kn, xx)# indices of knots within xx
pp <- predict(spl, xx, deriv = 1)
idx <- which(pp$y[i.kn]==max(pp$y[i.kn]))
idxx<- kn[idx]
df6a <- data.frame(Latency=pp$x,deriv=pp$y,varname="eccentricity",maxpoint=idxx,maxderiv=pp$y[idx])
#2)
spl <- smooth.spline(df$Latency[idxoverunderxecc], df$varest[idxoverunderxecc], df=12)
spl2 <- spl
xx <- unique(sort(c(seq(0,1100, by = 1), kn <- unique(df$Latency[idxoverunderxecc]))))
i.kn <- match(kn, xx)# indices of knots within xx
pp <- predict(spl, xx, deriv = 1)
idx <- which(pp$y[i.kn]==min(pp$y[i.kn]))
idxx<- kn[idx]
df6b <- data.frame(Latency=pp$x,deriv=pp$y,varname="underoverxecc",maxpoint=idxx,maxderiv=pp$y[idx])
#3)
spl <- smooth.spline(df$Latency[idxerror], df$varest[idxerror], df=12)
spl3 <- spl
xx <- unique(sort(c(seq(0,1100, by = 1), kn <- unique(df$Latency[idxerror]))))
i.kn <- match(kn, xx)# indices of knots within xx
pp <- predict(spl, xx, deriv = 1)
idx <- which(pp$y[i.kn]==max(pp$y[i.kn]))
idxx<- kn[idx]
df6c <- data.frame(Latency=pp$x,deriv=pp$y,varname="error",maxpoint=idxx,maxderiv=pp$y[idx])
#4)
spl <- smooth.spline(df$Latency[idxunderover], df$varest[idxunderover], df=12)
spl4 <- spl
xx <- unique(sort(c(seq(0,1100, by = 1), kn <- unique(df$Latency[idxunderover]))))
i.kn <- match(kn, xx)# indices of knots within xx
pp <- predict(spl, xx, deriv = 1)
idx <- which(pp$y[i.kn]==min(pp$y[i.kn]))
idxx<- kn[idx]
df6d <- data.frame(Latency=pp$x,deriv=pp$y,varname="underover",maxpoint=idxx,maxderiv=pp$y[idx])
# put all data.frames together
df6 <- rbind(df6x,df6a,df6b,df6c,df6d)
levels(df6$varname) <- c("Intercept","eccentricity (0=close; 1=distant)","under-/overshoot x ecc","error","under-/overshoot")
spl <- data.frame(Latency = rep(spl0$x,times=5),
spline  = c(spl0$y,spl1$y,spl2$y,spl3$y,spl4$y),
varname = rep(c("Intercept","eccentricity (0=close; 1=distant)","under-/overshoot x ecc","error","under-/overshoot"),each=length(spl0$x)
))
levels(df$varname) <- c("Intercept","eccentricity (0=close; 1=distant)","under-/overshoot x ecc","error","under-/overshoot")
df6$varname <- as.character(df6$varname)
df6$varname <- factor(df6$varname,levels=c("Intercept","error","eccentricity (0=close; 1=distant)","under-/overshoot","under-/overshoot x ecc"))
df$varname <- as.character(df$varname)
df$varname <- factor(df$varname,levels=c("Intercept","error","eccentricity (0=close; 1=distant)","under-/overshoot","under-/overshoot x ecc"))
spl$varname <- as.character(spl$varname)
spl$varname <- factor(spl$varname,levels=c("Intercept","error","eccentricity (0=close; 1=distant)","under-/overshoot","under-/overshoot x ecc"))
#-----------------------------------------------------
#Baseline
par(mfrow=c(2,5),mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(fit_aalen3$cum[,1],fit_aalen3$cum[,2],type="l",lwd=1,xlim=c(0,1200),ylim=c(-1.25,2),ylab="",xlab="Time",bty="n",tck=-0.04)
abline(h=0)
#Error
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(fit_aalen3$cum[,1],fit_aalen3$cum[,5],type="l",lwd=1,xlim=c(0,1200),ylim=c(-1.25,2),ylab="",yaxt="n",xlab="Time",bty="n",tck=-0.04)
abline(h=0)
#Underovershoot
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(fit_aalen3$cum[,1],fit_aalen3$cum[,4],type="l",lwd=1,xlim=c(0,1200),ylim=c(-1.25,2),ylab="",yaxt="n",xlab="Time",bty="n",tck=-0.04)
abline(h=0)
#Eccentricity
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(fit_aalen3$cum[,1],fit_aalen3$cum[,3],type="l",lwd=1,xlim=c(0,1200),ylim=c(-1.25,2),ylab="",yaxt="n",xlab="Time",bty="n",tck=-0.04)
abline(h=0)
#Underovershoot x Eccentricity
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(fit_aalen3$cum[,1],fit_aalen3$cum[,6],type="l",lwd=1,xlim=c(0,1200),ylim=c(-1.25,2),ylab="",yaxt="n",xlab="Time",bty="n",tck=-0.04)
abline(h=0)
#Baseline
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(df6x$Latency,df6x$deriv,xlim=c(0,1200),ylim=c(-0.01,0.015),type="l",lwd=1,ylab="",xlab="Time",bty="n",tck=-0.04)
lines(c(df6x$maxpoint[1],df6x$maxpoint[1]),c(df6x$maxderiv[1],-0.01),lty=2,col=defgrey)
abline(h=0)
#Error
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(df6c$Latency,df6c$deriv,xlim=c(0,1200),ylim=c(-0.01,0.015),type="l",lwd=1,ylab="",xlab="",yaxt="n",bty="n",tck=-0.04)
lines(c(df6c$maxpoint[1],df6c$maxpoint[1]),c(df6c$maxderiv[1],-0.01),lty=2,col=defgrey)
abline(h=0)
#Underovershoot
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(df6d$Latency,df6d$deriv,xlim=c(0,1200),ylim=c(-0.01,0.015),type="l",lwd=1,ylab="",xlab="",yaxt="n",bty="n",tck=-0.04)
lines(c(df6d$maxpoint[1],df6d$maxpoint[1]),c(df6d$maxderiv[1],-0.01),lty=2,col=defgrey)
abline(h=0)
#Eccentricity
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(df6a$Latency,df6a$deriv,xlim=c(0,1200),ylim=c(-0.01,0.015),type="l",lwd=1,ylab="",xlab="",yaxt="n",bty="n",tck=-0.04)
lines(c(df6a$maxpoint[1],df6a$maxpoint[1]),c(df6a$maxderiv[1],-0.01),lty=2,col=defgrey)
abline(h=0)
#Underovershoot x Eccentricity
par(mar=c(3,3,1,1)-0.5,ps=7,cex=1)
plot(df6b$Latency,df6b$deriv,xlim=c(0,1200),ylim=c(-0.01,0.015),type="l",lwd=1,ylab="",xlab="",yaxt="n",bty="n",tck=-0.04)
lines(c(df6b$maxpoint[1],df6b$maxpoint[1]),c(df6b$maxderiv[1],-0.01),lty=2,col=defgrey)
abline(h=0)
#-----------------------------------------------------
#---------------------------------------------------------------
# Survival Analysis - Make predictions
#---------------------------------------------------------------
#---------------------------------------------------------------
# Fit simple model
fit_simple  <- aalen(Surv(latency, event) ~ 1 + error_hor_abs, data = drev, max.time=1100, residuals=0, n.sim=0,resample.iid=1)
fit_complex <- aalen(Surv(latency, event) ~ 1 + error_hor_abs + underovercond + ecccond + ecccond*underovercond, data = drev, max.time=1100, residuals=0, n.sim=0,resample.iid=1)
par(mfrow=c(1,2))
plot(fit_simple)
#---------------------------------------------------------------
#---------------------------------------------------------------
# Prepare data.frame for prediction
newdata <- data.frame(error_hor_abs = rep(seq(0.5,2.5,length.out=5),times=4),underovercond = rep(c(-1,1),each=2*5),ecccond=rep(rep(c(-1,1),each=5),times=2))
out.prediction <- predict(fit_complex,newdata=newdata,resample.iid=1)
str(out.prediction)
#---------------------------------------------------------------
# Survival Analysis - Make predictions
#---------------------------------------------------------------
#---------------------------------------------------------------
# Fit simple model
fit_simple  <- aalen(Surv(latency, event) ~ 1 + error_hor_abs, data = drev, max.time=1100, residuals=0, n.sim=0,resample.iid=1)
fit_complex <- aalen(Surv(latency, event) ~ 1 + error_hor_abs + underovercond + ecccond + ecccond*underovercond, data = drev, max.time=1100, residuals=0, n.sim=0,resample.iid=1)
par(mfrow=c(1,2))
plot(fit_simple)
#---------------------------------------------------------------
#---------------------------------------------------------------
# Prepare data.frame for prediction
newdata <- data.frame(error_hor_abs = rep(seq(0.5,2.5,length.out=5),times=4),underovercond = rep(c(-1,1),each=2*5),ecccond=rep(rep(c(-1,1),each=5),times=2))
out.prediction <- predict(fit_complex,newdata=newdata,resample.iid=1)
# Grouping of error_hor
dx <- drev[which(drev$event==1),]
dx$error_hor_group <- cut(dx$error_hor_abs,breaks=c(quantile(dx$error_hor_abs,probs=c(0,0.5,1))),labels=c(1,2))
# Compute means and confidence intervals
m0 <- melt(dx,id=c("vp","error_hor_group","underovershoot","eccentricity"),measure=c("latency"))
c0 <- cast(m0,vp + error_hor_group + underovershoot + eccentricity ~ variable,median)
m1 <- melt(c0,id=c("error_hor_group","underovershoot","eccentricity"),measure=c("latency"))
c1 <- cast(m1,error_hor_group + underovershoot + eccentricity ~ variable,c(mean,sd,length))
c1$se <- c1$latency_sd/sqrt(c1$latency_length)
c1$ci_upper <- c1$latency_mean+1.96*c1$se
c1$ci_lower <- c1$latency_mean-1.96*c1$se
# Compute median collapsed over all data
m0 <- melt(dx,id=c("error_hor_group","underovershoot","eccentricity"),measure=c("latency"))
c0 <- cast(m0,error_hor_group + underovershoot + eccentricity ~ variable,length)
# Plot survival curves, density and means
par(mfrow=c(2,4),mar=c(3,1,1,1)-0.5,ps=7,cex=1)
#SURVIVAL
plot(out.prediction$time,out.prediction$S0[1,],ylim=c(0,1),type="n",lwd=1,ylab=c("Survival"),xlab="Time",bty="n",tck=-0.04)
for(i in 1:5){
lines(out.prediction$time,out.prediction$S0[i,],lty=2,col=defgrey)
}
par(mar=c(3,1,1,1)-0.5,ps=7,cex=1)
plot(out.prediction$time,out.prediction$S0[1,],ylim=c(0,1),type="n",lwd=1,ylab=c("Survival"),xlab="Time",bty="n",tck=-0.04)
for(i in 6:10){
lines(out.prediction$time,out.prediction$S0[i,],lty=2,col=defgrey)
}
par(mar=c(3,1,1,1)-0.5,ps=7,cex=1)
plot(out.prediction$time,out.prediction$S0[1,],ylim=c(0,1),type="n",lwd=1,ylab=c("Survival"),xlab="Time",bty="n",tck=-0.04)
for(i in 11:15){
lines(out.prediction$time,out.prediction$S0[i,],lty=2,col=defgrey)
}
par(mar=c(3,1,1,1)-0.5,ps=7,cex=1)
plot(out.prediction$time,out.prediction$S0[1,],ylim=c(0,1),type="n",lwd=1,ylab=c("Survival"),xlab="Time",bty="n",tck=-0.04)
for(i in 16:20){
lines(out.prediction$time,out.prediction$S0[i,],lty=2,col=defgrey)
}
#DENSITY
plot(density(dx$latency[which(dx$event==1 & dx$underovershoot=="undershoot" & dx$eccentricity=="close target" & dx$error_hor_group==2)],from=0,to=1100),col="black",lwd=1,ylim=c(0,0.018),ylab=c("Density"),xlab="Time",bty="n",tck=-0.04,main="")
lines(density(dx$latency[which(dx$event==1 & dx$underovershoot=="undershoot" & dx$eccentricity=="close target" & dx$error_hor_group==1)],from=0,to=1100),col="grey",lwd=1,ylab=c("Density"),xlab="Time",bty="n",tck=-0.04,main="")
abline(v=c(c0$latency[1],c0$latency[5]),col=c("grey","black"),lty=2)
plot(density(dx$latency[which(dx$event==1 & dx$underovershoot=="undershoot" & dx$eccentricity=="distant target" & dx$error_hor_group==2)],from=0,to=1100),col="black",lwd=1,ylim=c(0,0.018),ylab=c("Density"),xlab="Time",bty="n",tck=-0.04,main="")
lines(density(dx$latency[which(dx$event==1 & dx$underovershoot=="undershoot" & dx$eccentricity=="distant target" & dx$error_hor_group==1)],from=0,to=1100),col="grey",lwd=1,ylab=c("Density"),xlab="Time",bty="n",tck=-0.04,main="")
abline(v=c(c0$latency[2],c0$latency[6]),col=c("grey","black"),lty=2)
plot(density(dx$latency[which(dx$event==1 & dx$underovershoot=="overshoot" & dx$eccentricity=="close target" & dx$error_hor_group==2)],from=0,to=1100),col="black",lwd=1,ylim=c(0,0.018),ylab=c("Density"),xlab="Time",bty="n",tck=-0.04,main="")
lines(density(dx$latency[which(dx$event==1 & dx$underovershoot=="overshoot" & dx$eccentricity=="close target" & dx$error_hor_group==1)],from=0,to=1100),col="grey",lwd=1,ylab=c("Density"),xlab="Time",bty="n",tck=-0.04,main="")
abline(v=c(c0$latency[3],c0$latency[7]),col=c("grey","black"),lty=2)
plot(density(dx$latency[which(dx$event==1 & dx$underovershoot=="overshoot" & dx$eccentricity=="distant target" & dx$error_hor_group==1)],from=0,to=1100),col="grey",lwd=1,ylim=c(0,0.018),ylab=c("Density"),xlab="Time",bty="n",tck=-0.04,main="")
lines(density(dx$latency[which(dx$event==1 & dx$underovershoot=="overshoot" & dx$eccentricity=="distant target" & dx$error_hor_group==2)],from=0,to=1100),col="black",lwd=1,ylab=c("Density"),xlab="Time",bty="n",tck=-0.04,main="")
abline(v=c(c0$latency[4],c0$latency[8]),col=c("grey","black"),lty=2)
#---------------------------------------------------------------
