问题描述
我是R的新手.目前,我正在将对数正态分布拟合到我拥有的某些生存数据,但是,在尝试计算诸如中位数和均值之类的统计信息时,我陷入了困境.这是我到目前为止使用的代码,任何人都可以告诉我接下来我应该键入什么来找到均值吗?
I am a novice with R. Currently I am fitting a log-normal distribution to some survival data I have, however I have become stuck when trying to calculate statistics such as the median and the mean. This is the code I have used so far, can anyone tell me what I should type next to find the mean?
# rm(list=ls(all=TRUE))
library(survival)
data<-read.table("M:\\w2k\\Diss\\Hoyle And Henley True IPD with number at risk known.txt",header=T)
attach(data)
data
times_start <-c( rep(start_time_censor, n_censors), rep(start_time_event, n_events) )
times_end <-c( rep(end_time_censor, n_censors), rep(end_time_event, n_events) )
model <- survreg(Surv(times_start, times_end, type="interval2")~1, dist="lognormal")
intercept <- summary(model)$table[1]
log_scale <- summary(model)$table[2]
这是我被卡住的地方,我已经尝试过:
this is where I got stuck, I have tried:
meantime<-exp(intercept+log_scale/2)
但这似乎并不现实.
推荐答案
查找可行示例的地方是?predict.survreg
. (通常,将帮助系统用于predict
方法是任何回归方法的有效策略.)
The place to look for a worked example is ?predict.survreg
. (In general, using the help system for predict
methods is a productive strategy for any regression method.)
运行最后一个示例应该为您提供足够的基础进行下一步.特别要注意的是,回归系数不是生存时间或分位数的估计值.
Running the last example should give you enough basis to proceed. In particular you should see that the regression coefficients are not estimates of survival times or quantiles.
lfit <- survreg(Surv(time, status) ~ ph.ecog, data=lung)
pct <- 1:98/100 # The 100th percentile of predicted survival is at +infinity
ptime <- predict(lfit, newdata=data.frame(ph.ecog=2), type='quantile',
p=pct, se=TRUE)
matplot(cbind(ptime$fit, ptime$fit + 2*ptime$se.fit,
ptime$fit - 2*ptime$se.fit)/30.5, 1-pct,
xlab="Months", ylab="Survival", type='l', lty=c(1,2,2), col=1)
# The plot should be examined since you asked for a median survival time
abline(h= 0.5)
# You can drop a vertical from the intersection to get that graphically
....或...
str(ptime)
List of 2
$ fit : num [1:98] 9.77 16.35 22.13 27.46 32.49 ...
$ se.fit: num [1:98] 2.39 3.53 4.42 5.16 5.82 ...
您可以使用以下方法从该生存时间序列中提取第50个百分点:
You can extract the 50th percentile from that sequence of survival times with:
ptime$fit[which((1-pct)==0.5)]
# [1] 221.6023
以天为单位进行测量,这就是Therneau除以30.5来显示月份的原因
Measured in days which was why Therneau divided by 30.5 to display months
这篇关于在R的生存分析中寻找对数正态分布的均值的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!