2013-03-05 2 views
2

Мне нужно подгонять параметрическую модель PH (так, а не модель Кокса) с изменяющимися во времени ковариатами. Можем ли мы это сделать в R? Я слышал, что функция супрега не может обрабатывать изменяющиеся во времени ковариаты. Я тщетно искал пакеты, которые могли бы с этим справиться.Установка полностью параметрической модели пропорциональной опасности с изменяющимися во времени ковариатами в R

+0

Может ли этот пакет помочь: http://cran.r-project.org/web/packages/joineR/index.html –

+0

вы могли бы преобразовать свои данные и вставить glm или gam с помощью 'family = poisson' (см. Http : //data.princeton.edu/wws509/notes/c7s4.html) – adibender

+0

Terry Therneau, автор (??) пакета выживания не поставил его в качестве приоритета (см. https://stat.ethz.ch /pipermail/r-help/2008-May/161027.html). Это отлично работает в пакетах анализа статистической выживаемости Stata, если у вас есть доступ к этому. Насколько я знаю, SAS пока не использует данные стиля процесса подсчета в параметрических моделях выживаемости. EDIT: Терри Т. в прошлом реагировал на электронные письма, поэтому, возможно, вы могли бы подать запрос с разумным примером, в котором эта модель нужна - это подход, который будет полезен многим, я подозреваю. –

ответ

0

Как @adibender пишет, вы можете легко оценить модель с постоянной базой с семейством poisson с журнальным временем offset с glm. Вот пример

> # Input parameters 
> n <- 100   # Number of individuals 
> t_max <- 5  # max number of period per individual 
> beta <- c(-1, 1) # true coefficient 
> 
> # Simulate data 
> set.seed(47261114) 
> sim_dat <- replicate(
+ n, 
+ { 
+  out <- data.frame(
+  tstart = rep(NA_integer_, t_max), 
+  tstop = rep(NA_integer_, t_max), 
+  event = rep(NA, t_max), 
+  x  = rnorm(t_max)) 
+  
+  for(i in 1:t_max){ 
+  rate  <- exp(beta %*% c(1, out$x[i])) 
+  tstop <- min(rexp(1, rate), 1) 
+  out[i, ] <- list(i - 1, i - (1 - tstop), tstop < 1, out$x[i]) 
+  if(out$event[i]) 
+   break 
+  } 
+  out[!is.na(out$tstart), ] 
+ }, simplify = FALSE) 
> 
> sim_dat <- do.call(rbind, sim_dat) 
> head(sim_dat) # show final data 
    tstart  tstop event   x 
1  0 0.3018182 TRUE 0.7095841 
2  0 0.6724803 TRUE 1.5152877 
3  0 1.0000000 FALSE 0.1036868 
4  1 2.0000000 FALSE -0.5214508 
5  2 2.4831577 TRUE 1.0101403 
6  0 1.0000000 FALSE 0.1437594 
> 
> # Fit with glm 
> glm(event ~ x + offset(log(tstop - tstart)), sim_dat, family = poisson()) 

Call: glm(formula = event ~ x + offset(log(tstop - tstart)), family = poisson(), 
    data = sim_dat) 

Coefficients: 
(Intercept)   x 
    -0.9053  0.9714 

Degrees of Freedom: 248 Total (i.e. Null); 247 Residual 
Null Deviance:  382.5 
Residual Deviance: 306.4 AIC: 498.4 

Для другого распределения, кажется, что flexsurv пакет может сделать. См. this post.

 Смежные вопросы

  • Нет связанных вопросов^_^