1st-生存分析R语言-Survival Analysis

1260

收藏 2020-10-16

library(survival)
n <- 200
liferate <- 0.01
life <- rexp(n, liferate)  #生成指数分布的随机数200个
round(life[1:10], 1)    #取前10个保留1位小数

censrate <- 0.005
cens <- rexp(n, censrate)
round(cens[1:10], 1) #同上

oldpar <- par(mfrow=c(1, 2))    #用par作两张图
hist(life, main="True lifetimes")
hist(cens, main="Censoring times")

par(oldpar) #重制作图参数

delta <- rep(0, n)
delta[life < cens] <- 1 #把1插入到cens大于life的部分
delta[1:10]          #预览前十个

sum(delta)/n #失败案例的比例

plot(life, cens, main="Censoring times vs. True lifetimes", xlab="lifetimes",
   ylab="censoring times")
abline(0, 1, col="red", lwd=2)
points(life[delta==1], cens[delta==1], pch=16)

times <- pmin(life, cens)  #两者选择更小的数据载入
round(times[1:10], 1)

median.true <- qexp(p=0.5, rate=liferate)
median.true

median.all <- median(times)    #全体数据的中位数
median.all

median.unc <- median(times[delta == 1]) #未删失的观测值的中位数
median.unc

km <- survfit(Surv(times, delta) ~ 1)
plot(km, main="Kaplan-Maier estimate of S(t)", xlab="time")
abline(h=0.5, col="red")

print(km)

median.Surv <- quantile(km, probs=0.5, conf.int=FALSE)
median.Surv

扫码加我拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

全部回复

miracle丶

2020-10-18 05:57:19

后半部分，说实话笔者没看懂这是在干嘛。。。好像是在比较删失数据和不删失数据的准确性？
望能抛砖引玉，才疏学浅，还请多多指教！

n <- 200             # sample size
nsim <- 1000          # simulation sample size
median.all <- median.unc <- median.Surv <- rep(0, nsim)  # vectors to store results
liferate <- 0.01; censrate <- 0.005       # lambda_T and lambda_C
median.true <- qexp(p=0.5, rate=liferate)    # true median
for (j in (1:nsim))
{
  life <- rexp(n, liferate)       # true lifetimes
  cens <- rexp(n, censrate)    # censoring times
  delta <- rep(0, n)          # censoring indicator
  delta[life < cens] <- 1       # 1= uncensored
  times <- pmin(life, cens) # observations (censored or not)
  median.all[j] <- median(times)    # median of all obs (cens + uncens)
  median.unc[j] <- median(times[delta == 1]) # median of uncensored obs.
  km <- survfit(Surv(times, delta) ~ 1)
  median.Surv[j] <- quantile(km, probs=0.5, conf.int=FALSE)
}

hist(median.all, xlim=range(median.all, median.true))
points(median.true, 0, pch=16, col="red", cex=1.5)

hist(median.unc, xlim=range(median.unc, median.true))
points(median.true, 0, pch=16, col="red", cex=1.5)

( MSE.all <- mean((median.all - median.true)^2) ) #MSE均方误差
( MSE.unc <- mean((median.unc - median.true)^2) )
( MSE.Surv <- mean((median.Surv - median.true)^2) )

扫码加我拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

扫码加我 拉你入群

扫码加我 拉你入群

分享

扫码加好友，拉您进群

扫码加我拉你入群

扫码加我拉你入群