# TODO: Add comment
# 
# Author: E.Korsching Sep 6, 2012
###############################################################################



#### survival objects

library(survival)

## samples
# fit a Kaplan-Meier and plot it 

#> aml
#	time status             x
#1     9      1    Maintained
#2    13      1    Maintained
#...
#23   45      1 Nonmaintained

afit1 <- survfit(Surv(time, status) ~ x, data=aml) 
plot(afit1, lty = 2:3) 
legend(100, .8, c("Maintained", "Nonmaintained"), lty = 2:3) 


# fit a Cox proportional hazards model and plot the predicted survival for a 60 year old

#> ovarian
#	futime fustat     age resid.ds rx ecog.ps
#1      59      1 72.3315        2  1       1
#2     115      1 74.4932        2  1       1
#...
#26    377      0 58.3096        1  2       1

ab <- names(ovarian)
attributes(ovarian)

afit2 <- coxph(Surv(futime, fustat) ~ age, data=ovarian) 
afit2 <- coxph(as.formula(paste("Surv(futime, fustat) ~", ab[3])), data=ovarian) 
afit2 <- coxph(Surv(ovarian[,1], ovarian[,2]) ~ ovarian[,3])		#rechnet zwar, aber geht nicht, weil der Name verlohren geht
plot(survfit(afit2, newdata=data.frame(age=60)), xscale=365.25, xlab="Years", ylab="Survival") 


#> attributes(object)		#survfit()
#$names
#[1] "n"         "time"      "n.risk"    "n.event"   "n.censor"  "surv"     
#[7] "type"      "strata"    "std.err"   "upper"     "lower"     "conf.type"
#[13] "conf.int"  "call"     
#
#$class
#[1] "survfit"
#
#> attributes(object)		#coxph()
#$names
#[1] "coefficients"      "var"               "loglik"           
#[4] "score"             "iter"              "linear.predictors"
#[7] "residuals"         "means"             "concordance"      
#[10] "method"            "n"                 "nevent"           
#[13] "terms"             "assign"            "wald.test"        
#[16] "y"                 "formula"           "call"             
#
#$class
#[1] "coxph"


# survfit - difference significant ?   log-rank test

adiff1 <- survdiff(Surv(time, status) ~ x, data=aml)

#> attributes(object)		#survdiff
#$names
#[1] "n"     "obs"   "exp"   "var"   "chisq" "call" 
#
#$class
#[1] "survdiff"

# p value by exec of adiff1  or  explicitely by :
p.val <- 1 - pchisq(adiff1$chisq, length(adiff1$n) - 1)


#### coxph

data(burn); attach(burn)
my.surv <- Surv(T1, D1)
coxph(my.surv ~ Z1 + as.factor(Z11), method='breslow')	# output omitted


####
#local.coxph.test(coxph.fit, pos, C=NA, d=NA, sign.digits=3)
#
#This function is used to compute a local p-value. Specifically, it is probably most useful
#to use it to find the p-value of an 'as.factor' variable but it can also be used to find
#the p-value of a many other more complex local tests.
#Arguments
#
#coxph.fit :: This argument should be a coxph object. That is, it should be the output
#from the function coxph()
#
#pos :: Given the output of coxph.fit, this is a vector of the position numbers of the
#explanatory variables (covariates) that will be in the local test.
#
#C and d :: These two variables define what kind of hypothesis is being checked.
#Specifically, the null hypothesis is of the form: C*beta.hat = d. The default for
#C is the identity matrix, and the default of d is a vector of zeros. That is, the
#default test is a global test on the variables in pos.
#
#sign.digits :: The number of significant digits/figures to output for the p-value.

#===> extract a p-value from coxph <===#
local.coxph.test <- function(coxph.fit, pos, C=NA, d=NA, sign.digits=3){
	# checking H_0: C*beta.hat = c
	# coxph.fit is a coxph model
	# pos is the position numbers of the parameters of interest
	# C is a q-by-p matrix
	# d is a q-by-1 matrix
	n <- length(pos)
	if(is.na(C)){
		C <- matrix(0, n, n)
		for(i in 1:n)
			C[i,i] <- 1
	} else {
		if(dim(C)[1] != n)
			stop("C has improper dimensions\n")
	}
	if(is.na(d))
		d <- matrix(0, n, 1)
	if(dim(d)[1] != dim(C)[1])
		stop("C and d do not have appropriate dimensions\n")
	I. <- coxph.fit$var[pos,pos]
	est <- matrix(as.vector(coxph.fit$coeff[pos]), dim(C)[2])
	X <- as.numeric( t(C%*%est - d) %*%
					solve( t(C) %*% I. %*% C )  %*%
					(C%*%est - d) )
	signif(1-pchisq(X, dim(C)[1]), sign.digits)
}
local.p.test <- local.coxph.test


#### John Fox 2002 coxph .pdf
Rossi <- read.table("/home/korschi/Desktop/Survival_R/Rossi.txt", header=TRUE, sep="\t", dec=".", stringsAsFactors=FALSE)

md <- coxph(Surv(week, arrest) ~ fin + age + race + wexp + mar + paro + prio, data=Rossi)
md
summary(md)
md$loglik		#Likelihood ratio test ??zwei Zahlen
md$wald.test	#Wald test
md$score		#Score (logrank) test

plot(survfit(md), ylim=c(.7, 1), xlab="Weeks", ylab="Proportion Not Rearrested")	# only mean plot of all covariates

# Figure 2: Estimated survival functions for those receiving (fin = 1) and not receiving (fin = 0) financial aid.
# Other covariates are fixed at their average values. Each estimate is accompanied by a point-wise 95-percent
# confidence envelope.

Rossi.fin <- data.frame( fin=c(0,1), age=rep(mean(Rossi$age),2), race=rep(mean(Rossi$race),2), wexp=rep(mean(Rossi$wexp),2),
		mar=rep(mean(Rossi$mar),2), paro=rep(mean(Rossi$paro),2), prio=rep(mean(Rossi$prio),2) )

plot( survfit(md, newdata=Rossi.fin), conf.int=T, lty=c(1,2), ylim=c(.6, 1) )
legend(x="bottomleft", legend=c("fin = 0", "fin = 1"), lty=c(1,2))


###### multiple covariates

# subset the data to include the time variable, indicator, and all other variables you want in the model, and then Surv(time, event) ~ .