2013-09-15 4 views
2

Я новичок с R. В настоящее время я подгоняю лог-нормальное распределение некоторых данных о выживании, которые у меня есть, однако я застрял при попытке рассчитать статистику таких как медиана и среднее. Это код, который я использовал до сих пор, может ли кто-нибудь сказать мне, что я должен напечатать рядом, чтобы найти среднее значение?Поиск среднего значения логарифмически нормального распределения в анализе выживаемости в R

# 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] 

это где я застрял, я попытался:

meantime<-exp(intercept+log_scale/2) 

, но это, кажется, не дают реалистичное среднее.

+3

Было бы очень полезно, если бы вы включили копию своих данных с 'dput (data)'. Если ваши данные очень большие, возможно, вы можете включить первые 50 строк с 'dput (head (data))'. Кроме того, в этом вопросе укажите фактическое число, которое вы получили за это время, и то, что вы ожидали. – nograpes

ответ

1

Место для поиска обработанного примера - ?predict.survreg. (В общем случае использование справочной системы для методов predict является продуктивной стратегией для любого метода регрессии.)

Запуск последнего примера должен дать вам достаточно оснований для продолжения. В частности, вы должны видеть, что коэффициенты регрессии не являются оценками времени выживания или квантилей.

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-процентиль из этой последовательности времени выживания с:

ptime$fit[which((1-pct)==0.5)] 
# [1] 221.6023 

Измеряется в днях, который был почему Therneau деленная на 30,5 для отображения месяцев