Sunday, September 14, 2014

9 14 2014 Kaplan Meier Graphs Improved

KMgraphsWorking.R
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))

plot of chunk unnamed-chunk-1

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))

plot of chunk unnamed-chunk-1

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))

plot of chunk unnamed-chunk-1

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))

plot of chunk unnamed-chunk-1

No comments:

Post a Comment