# 図1と図2のような図示をするためのコード。 # developed by Minato Nakazawa scale <- function(X) { # 2年未満を0-0.5に,2年から20年までの18年間を0.5-1にスケーリング ifelse(X<2*365.24,0.5*X/(2*365.24),0.5+0.5*(X-2*365.24)/(18*365.24)) } dat <- read.delim("surv.txt") # データの読み込み datop <- subset(dat,!is.na(ageop1)) # ageop1(手術日)が欠損でないデータのみのサブセット datnop <- subset(dat,is.na(ageop1)) # ageop1(手術日)が欠損なデータのみのサブセット datdead <- subset(dat,dead==1) # 死亡例だけのサブセット # # ** variables ** Patient,agepres,agelast,ageop1,dead,sex,paanat,adfol, # followup,opfpres,unopage,unopfpre,preopded,hadop,dedlyrpp,agepresx # # 年齢そのまま。ただし2年目までを長くスケーリング png("fig1fig2.png",width=640,height=1280) par(mfrow=c(2,1),las=1) plot(scale(dat$agepres),31-dat$Patient,axes=F,pch=ifelse(dat$paanat==0,5,18), xlab="Age in years",ylab="Patient",xlim=c(0,1), main="Fig.1. Histories of 30 selected patients from birth") points(scale(datop$ageop1),31-datop$Patient,pch=1) points(scale(datdead$agelast),31-datdead$Patient,pch=3) axis(1,scale(0:20*365.24),0:20) axis(2,1:30,30:1) segments(scale(datnop$agepres),31-datnop$Patient,scale(datnop$agelast),31-datnop$Patient,lty=1) segments(scale(dat$agepres),31-dat$Patient,scale(dat$ageop1),31-dat$Patient,lty=1) segments(scale(dat$ageop1),31-dat$Patient,scale(dat$agelast),31-dat$Patient,lty=4) xgrids <- scale(c(1,2,10)*365.24) segments(xgrids,0,xgrids,30,lty=3) legend("right",c("first seen, PA anat=0","first seen, PA anat=1","operation","death"), pch=c(5,18,1,3)) legend(scale(11.5*365.24),12.6,c("pre-op follow-up","post-op follow-up"),lty=c(1,4)) # 出現日をゼロに揃える timepres <- dat$agepres - dat$agepres timeoppres <- datop$agepres - datop$agepres timenoppres <- datnop$agepres - datnop$agepres timeop1 <- dat$ageop1 - dat$agepres timeopop1 <- datop$ageop1 - datop$agepres timelast <- datdead$agelast - datdead$agepres timeoplast <- datop$agelast - datop$agepres timenoplast <- datnop$agelast - datnop$agepres plot(scale(timepres),31-dat$Patient,axes=F,pch=ifelse(dat$paanat==0,5,18), xlab="Years after presentation",ylab="Patient",xlim=c(0,1), main="Fig.2. Histories of 30 selected patients from presentation") points(scale(timeopop1),31-datop$Patient,pch=1) points(scale(timelast),31-datdead$Patient,pch=3) axis(1,scale(0:20*365.24),0:20) axis(2,1:30,30:1) segments(scale(timenoppres),31-datnop$Patient,scale(timenoplast),31-datnop$Patient,lty=1) segments(scale(timeoppres),31-datop$Patient,scale(timeopop1),31-datop$Patient,lty=1) segments(scale(timeopop1),31-datop$Patient,scale(timeoplast),31-datop$Patient,lty=4) xgrids <- scale(c(1,2,10)*365.24) segments(xgrids,0,xgrids,30,lty=3) legend("right",c("first seen, PA anat=0","first seen, PA anat=1","operation","death"), pch=c(5,18,1,3)) legend(scale(11.5*365.24),12.6,c("pre-op follow-up","post-op follow-up"),lty=c(1,4)) dev.off()