9 14 2014 Kaplan Meier Graphs Improved

KMgraphsWorking.R

Jake H

Sun Sep 14 14:41:29 2014

library(survival)

## Loading required package: splines

kmdab=read.csv("KMdataDabob.csv")
names(kmdab)    

## [1] "Animal"     "Population" "Death"      "Status"

with(kmdab, tapply(Death[Status==1],Population[Status==1],mean))

##     H     N     S 
## 4.385 4.522 4.715

with(kmdab, tapply(Death[Status==1],Population[Status==1],var))

##      H      N      S 
## 0.4524 0.9786 1.6039

fit1=with(kmdab,survfit(Surv(Death,Status)~Population))
plot(fit1, col=c(1:3), xlab="Survival Time from Outplant in Months", ylab="% Surviving")
legend("bottomleft", title="Population", c("Dabob","Fidalgo","Oyster Bay"), fill=c(1:3))

kmman=read.csv("KMdataMan.csv")
names(kmman)

## [1] "Animal"     "Population" "Death"      "Status"

with(kmman, tapply(Death[Status==1],Population[Status==1],mean))

##      H      N      S 
## 10.540  9.368  9.197

with(kmman, tapply(Death[Status==1],Population[Status==1],var))

##     H     N     S 
## 1.519 6.468 6.321

fit2=with(kmman, survfit(Surv(Death,Status)~Population))
plot(fit2, col=c(1:3), xlab="Survival Time from Outplant in Months", ylab="% Surviving")
legend("bottomleft", title="Population", c("Dabob","Fidalgo","Oyster Bay"), fill=c(1:3))

kmfid=read.csv("KMdataFid.csv")
with(kmfid, tapply(Death[Status==1],Population[Status==1],mean))

##     H     N     S 
## 5.921 5.754 6.897

with(kmfid, tapply(Death[Status==1],Population[Status==1],var))

##     H     N     S 
## 2.554 2.189 3.989

fit3=with(kmfid, survfit(Surv(Death,Status)~Population))
plot(fit3, col=c(1:3), xlab="Survival Time from Outplant in Months", ylab="% Surviving")
legend("bottomleft", title="Population", c("Dabob","Fidalgo","Oyster Bay"), fill=c(1:3))

kmoys=read.csv("KMdataOys.csv")
with(kmoys, tapply(Death[Status==1],Population[Status==1],mean))

##     H     N     S 
## 7.550 8.069 8.396

with(kmoys, tapply(Death[Status==1],Population[Status==1],var))

##     H     N     S 
## 7.774 5.805 5.814

fit4=with(kmoys, survfit(Surv(Death,Status)~Population))
plot(fit4, col=c(1:3), xlab="Survival Time from Outplant in Months", ylab="% Surviving")
legend("bottomleft", title="Population", c("Dabob","Fidalgo","Oyster Bay"), fill=c(1:3))