rm(list=ls())

library(survival)
library(KMsurv)
#library(OIsurv)

#####set directory ####
setwd<-"C:/Users/jsweka/Documents/BrookTrout/SurvivalEstimation"

data.in<-read.csv("troutsurvival_daily.csv")


#######Brook Trout Survival - overall, no covariates##############
bkt<-subset(data.in, species == "brook")
head(bkt)

s<-Surv(bkt$time1, bkt$time2, bkt$death)#creates a survival object

fit<-survfit(s~1, conf.int = 0.95, conf.type="log")#fits kaplan-meier survival curve

summary(fit)

plot(fit, xlab = "Day", ylab = "Probability of surviving to day", main = "Brook Trout")

####Brook Trout Hazard Function #####
H.hat<- -log(fit$surv)
plot(fit$time,H.hat, xlab = "Day", ylab = "Cumulative hazard",
	main = "Brook Trout", type = "s")


##############Brook Trout Survival - stream order covariate#############
order_effect<-coxph(s~bkt$order)
summary(order_effect)

##############Brook Trout Survival - cum drainage covariate#############
cumdrain_effect<-coxph(s~bkt$cumdrain)
summary(cumdrain_effect)

##############Brook Trout Survival - stream covariate#############
stream_effect<-coxph(s~bkt$stream)
summary(stream_effect)

##############Brook Trout Survival - fish size covariate#############
size_effect<-coxph(s~bkt$Wgt)
summary(size_effect)

##############Brook Trout Survival - movement rate covariate#############
move_effect<-coxph(s~bkt$move_rate)
summary(move_effect)

##################################################################################################################

#######Brown Trout Survival - overall, no covariates##############
bnt<-subset(data.in, species == "brown")
head(bnt)

s<-Surv(bnt$time1, bnt$time2, bnt$death)#creates a survival object

fit<-survfit(s~1, conf.int = 0.95, conf.type="log")#fits kaplan-meier survival curve

summary(fit)

plot(fit, xlab = "Day", ylab = "Probability of surviving to day", main = "Brown Trout")

####Brown Trout Hazard Function #####
H.hat<- -log(fit$surv)
plot(fit$time,H.hat, xlab = "Day", ylab = "Cumulative hazard",
	main = "Brook Trout", type = "s")


##############Brown Trout Survival - stream order covariate#############
order_effect<-coxph(s~bnt$order)
summary(order_effect)

##############Brown Trout Survival - cum drainage covariate#############
cumdrain_effect<-coxph(s~bnt$cumdrain)
summary(cumdrain_effect)

##############Brown Trout Survival - stream order covariate#############
stream_effect<-coxph(s~bnt$stream)
summary(stream_effect)

##############Brown Trout Survival - fish size covariate#############
size_effect<-coxph(s~bnt$Wgt)
summary(size_effect)

##############Brown Trout Survival - movement rate covariate#############
move_effect<-coxph(s~bnt$move_rate)
summary(move_effect)



###########################################################################################################

##################species combined - overall, no covariates######################################

s<-Surv(data.in$time1, data.in$time2, data.in$death)#creates a survival object

fit<-survfit(s~1, conf.int = 0.95, conf.type="log")#fits kaplan-meier survival curve

summary(fit)

plot(fit, xlab = "Day", ylab = "Probability of surviving to day", main = "Species Combined")


############## species covariate################################
species_effect<-coxph(s~data.in$species)
summary(species_effect)

##############  fish size covariate ############################
size_effect<-coxph(s~data.in$Wgt)
summary(size_effect)

############## stream covariate  ###############################
stream_effect<-coxph(s~data.in$stream)
summary(stream_effect)


############# stream order covariate  #########################
order_effect<-coxph(s~data.in$order)
summary(order_effect)

############## cum drainage covariate#############
cumdrain_effect<-coxph(s~data.in$cumdrain)
summary(cumdrain_effect)

############# movement rate covariate #########################
move_effect<-coxph(s~data.in$move_rate)
summary(move_effect)


#############  Wilcoxon tests for species differences in movement rate and size  ##############
wilcox.test(move_rate~species, data = data.in)
wilcox.test(Wgt~species, data = data.in)