2014-03-27 51 views
2

我使用函数coxme()在coxme软件包中拟合R中的混合效果Cox模型。在我的模型中,我有一个审查生存时间$ X $,一个单一的协变量$ Z $和一个分组变量$ Group $。有一个随机的截距和一个随机的斜率,即我拟合了模型$ \ lambda(t | Z,b_0,b_1)= \ lambda_0(t)e^{\ beta Z + b_0 + b_1 Z} $。在混合效果中指定变异结构R中的Cox模型

我想指定随机效应的协方差矩阵的结构;特别是我想指出$ 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 
+0

我不知道该怎么做,但几个月前我用了'coxme'package,但有一个问题在这里无法解答。因此,我发送了一封电子邮件给Terry Therneau,这个软件包的开发人员,他回答我的速度很快,很有用,所以也许问他吧! – Vincent

回答

1

多google搜索和尝试后,我花了二看one of the vignettescoxme包。我在第3节的最后一点找到了一个可能的答案,它说

“ - 通过默认,假设一个完整的协方差矩阵。模型2显示指定独立性的简单方法是将以单独的方式效果。“在数据示例

所以,下面的命令可以被用来获得独立的随机效应:

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

$Group 
     Z 
0.2822648 

只有随机截距和随机斜率的方差被示出,因为所述相关性假定为零。

可能的替代方案是使用phmm包中的phmm()函数。由于估计是基于完全可能性的,所以这种方法会给出略微不同的拟合。

相关问题