Определение структуры дисперсии в модели Кокса со смешанными эффектами в R

Я устанавливаю модель Кокса со смешанными эффектами в R, используя функцию coxme() в пакете coxme. В моей модели у меня есть цензурированное время выживания $ X $, единственная ковариата $ Z $ и группирующая переменная $ Group $. Есть случайный перехват и случайный наклон, т.е. я подбираю модель $ \ lambda (t | Z, b_0, b_1) = \ lambda_0 (t) e ^ {\ beta Z + b_0 + b_1 Z} $.

Хочу уточнить структуру ковариационной матрицы случайных эффектов; в частности, я хотел бы указать, что $ b_0 $ и $ b_1 $ некоррелированы, так что ковариационная матрица диагональна. Кажется, что coxme() довольно гибко определяет структуру ковариации. Я думаю, что мне следует передать список матриц параметру varlist этой функции, но пока мои попытки не увенчались успехом, и я не думаю, что полностью понимаю, как это должно работать.

Также можно передать в эту опцию настраиваемую функцию отклонения, и на самом деле одна из виньеток пакета дает несколько примеров того, как это можно сделать. Однако этот процесс кажется утомительным и (надеюсь) ненужным в данном случае. Итак, мой вопрос: как я могу легко указать структуру диагональной ковариационной матрицы в функции coxme()?

Ниже приведены некоторые смоделированные примеры данных и первая попытка определения ковариационной структуры. Я надеялся, что я сказал coxme() использовать линейную комбинацию $ V = \ sigma ^ 2_1 A + \ sigma ^ 2_2 B $, с $ A $ и $ B $, как определено ниже, и что это будет эффективно соответствовать диагональной ковариации матрица с произвольными диагональными элементами.

> n = 25  # Size of each cluster
> K = 25  # Number of clusters
> N = n*K  # Total number of observations
> 
> Z = rnorm(n=N, mean=0.5, sd=0.5)  # Covariate
> b0 = rep(rnorm(n=K, mean=0, sd=0.5), each=n)  # Random intercept
> b1 = rep(rnorm(n=K, mean=0, sd=0.5), each=n)  # Random slope
> Group = factor(x=rep(1:K, each=n))
> 
> beta = 2
> eta = beta*Z + b0 + b1*Z
> T = rexp(n=N, rate=exp(eta))  # Exponential failure time, conditional on Z, b0, and b1
> C = runif(n=N, min=0, max=2.5)  # Uniform censoring time to get about 20% censoring
> 
> time = pmin(T,C)  # Censored observation time
> delta = T < C  # Event indicator
> 
> A = matrix(c(1, 0, 0, 0), nrow=2)
> B = matrix(c(0, 0, 0, 1), nrow=2)
> my.covariance = list(A, B)
> fit = coxme(Surv(time, delta) ~ Z + (1 + Z | Group), varlist = my.covariance)
Error in coxme(Surv(time, delta) ~ Z + (1 + Z | Group), varlist = my.covariance) : 
  In random term 1: Mlist cannot have both covariates and grouping

person Breaking Waves    schedule 27.03.2014    source источник
comment
Я не знаю, как это сделать, но несколько месяцев назад я использовал coxmepackage, и у меня возник вопрос, на который нельзя было ответить здесь. Поэтому я отправил электронное письмо Терри Терно, разработчику пакета, и он ответил мне довольно быстро и полезно, так что, возможно, спросите его!   -  person Vincent    schedule 06.05.2014


Ответы (1)


После долгих поисков в Google и экспериментов я еще раз взглянул на одну из виньеток для пакета coxme. Я нашел один из возможных ответов на проблему в последнем пункте раздела 3, в котором говорится:

"- По умолчанию предполагается полная ковариационная матрица. Модель 2 показывает, что простой способ указать независимость - это поместить эффекты в отдельные термины."

Итак, в примере данных для получения независимых случайных эффектов можно использовать следующую команду:

> fit = coxme(Surv(time, delta) ~ Z + (1 | Group) + (Z | Group))
> fit$vcoef
$Group
Intercept 
0.1181417 

$Group
        Z 
0.2822648 

Показаны только дисперсии случайного пересечения и случайного наклона, поскольку предполагается, что корреляция равна нулю.

Возможная альтернатива - использовать функцию phmm() в пакете phmm. Этот метод даст несколько иное соответствие, поскольку оценка основана на полном правдоподобии.

person Breaking Waves    schedule 04.04.2014