随机斜率 Cox 比例风险
Random slopes Cox Proportional Hazards
我一直在尝试使用 coxme
为模型中的每个协变量提取随机斜率。
library (coxme)
Start <- runif(5000, 1985, 2015)
Stop <- Start + runif(5000, 2, 10)
S <- data.frame (
X1 <- runif(5000, 5.0, 7.5),
X2 <- runif(5000, 5.0, 7.5),
D <- rbinom(5000, 1, 0.8),
Letters <- sample(LETTERS, 5000, replace = TRUE),
Start <- Start,
Stop <- Stop
)
S_ind1 <- Surv (time = S$Start, time2 = S$Stop, event = S$D)
a <- coxme (S_ind1 ~ X1 + X2 + (X1 + X2|Letters), data = S)
我得到的是:
Error in gchol(kfun(theta, varlist, vparm, ntheta, ncoef)) :
NA/NaN/Inf in foreign function call (arg 5)
In addition: Warning messages:
1: In sqrt(xvar * zvar) : NaNs produced
2: In sqrt(xvar * zvar) : NaNs produced
使用我自己的数据时,我经常得到:
Error in coxme.fit(X, Y, strats, offset, init, control, weights = weights, :
'Calloc' could not allocate memory (56076596 of 8 bytes)
是否可以使用 coxme
来包含随机坡度?
如果没有,是否有使用其他包的其他替代方案?
coxme
包的作者 Terry Therneau 通过电子邮件回答 - 他让我post 在这里。
下面是我对您的示例的重写,删除了 Surv 间接寻址并在 data.frame 调用中使用了“=”(我有点惊讶 <- 在该上下文中有效) ,并添加 set.seed 以便示例可重现。
library (coxme)
set.seed(1953)
time1 <- runif(5000, 1985, 2015)
time2 <- time1 + runif(5000, 2, 10)
test <- data.frame (
x1 = runif(5000, 5.0, 7.5),
x2 = runif(5000, 5.0, 7.5),
death = rbinom(5000, 1, 0.8),
letters = sample(LETTERS, 5000, replace = TRUE),
time1 = time1,
time2 = time2)
fit1 <- coxme(Surv(time1, time2, death) ~ x1 + x2 + (1|letters), data=test)
fit2 <- coxme(Surv(time1, time2, death) ~ x1 + x2 + (1+x1 | letters), test)
fit3 <- coxme(Surv(time1, time2, death) ~ x1 + x2 + (1+x2 | letters), test)
fit4 <- coxme(Surv(time1, time2, death) ~ x1 + x2 + (1+ x1 + x2 | letters),
data=test, vinit= c(1e-6, 1e-8, 1e-8))
*1。所有模型都可以工作,直到 fit4.
我发现你的模型令人担忧,因为它有一个随机斜率但没有随机截距,就像所有通过原点的回归让我担心的一样:我很难解释结果。尽管 lme 默认包含截距项,但 coxme 没有。
我曾希望 fit4 能奏效,也许有了更好的初始估计它就会奏效。 coxme 的底层代码是我在所有生存工作中遇到的最难的最大化问题,从最大化器很容易迷路并且永远找不到它的方式的意义上来说很难。这是一个有时需要手动操作的函数,通过有限的迭代计数 and/or 开始估计。我希望不是这样,我有一些长期计划通过添加一个替代的基于 MCMC 的最大化器来改进这一点,理论上它永远不会丢失但以更长的计算时间为代价。
如果任何方差太接近于零,则 sqrt() 消息往往会作为舍入误差的函数出现。当然,在您的测试用例中,实际 MLE 的方差为 0。发生这种情况时,我通常会通过使用一系列固定方差(vfixed 参数)进行拟合来直接检查零方差。如果可能性不变或随着方差值达到 1e-6 或更小而增加,那么我假设 MLE 为零并从模型中删除该随机项。
特里 T.*
我一直在尝试使用 coxme
为模型中的每个协变量提取随机斜率。
library (coxme)
Start <- runif(5000, 1985, 2015)
Stop <- Start + runif(5000, 2, 10)
S <- data.frame (
X1 <- runif(5000, 5.0, 7.5),
X2 <- runif(5000, 5.0, 7.5),
D <- rbinom(5000, 1, 0.8),
Letters <- sample(LETTERS, 5000, replace = TRUE),
Start <- Start,
Stop <- Stop
)
S_ind1 <- Surv (time = S$Start, time2 = S$Stop, event = S$D)
a <- coxme (S_ind1 ~ X1 + X2 + (X1 + X2|Letters), data = S)
我得到的是:
Error in gchol(kfun(theta, varlist, vparm, ntheta, ncoef)) :
NA/NaN/Inf in foreign function call (arg 5)
In addition: Warning messages:
1: In sqrt(xvar * zvar) : NaNs produced
2: In sqrt(xvar * zvar) : NaNs produced
使用我自己的数据时,我经常得到:
Error in coxme.fit(X, Y, strats, offset, init, control, weights = weights, :
'Calloc' could not allocate memory (56076596 of 8 bytes)
是否可以使用 coxme
来包含随机坡度?
如果没有,是否有使用其他包的其他替代方案?
coxme
包的作者 Terry Therneau 通过电子邮件回答 - 他让我post 在这里。
下面是我对您的示例的重写,删除了 Surv 间接寻址并在 data.frame 调用中使用了“=”(我有点惊讶 <- 在该上下文中有效) ,并添加 set.seed 以便示例可重现。
library (coxme)
set.seed(1953)
time1 <- runif(5000, 1985, 2015)
time2 <- time1 + runif(5000, 2, 10)
test <- data.frame (
x1 = runif(5000, 5.0, 7.5),
x2 = runif(5000, 5.0, 7.5),
death = rbinom(5000, 1, 0.8),
letters = sample(LETTERS, 5000, replace = TRUE),
time1 = time1,
time2 = time2)
fit1 <- coxme(Surv(time1, time2, death) ~ x1 + x2 + (1|letters), data=test)
fit2 <- coxme(Surv(time1, time2, death) ~ x1 + x2 + (1+x1 | letters), test)
fit3 <- coxme(Surv(time1, time2, death) ~ x1 + x2 + (1+x2 | letters), test)
fit4 <- coxme(Surv(time1, time2, death) ~ x1 + x2 + (1+ x1 + x2 | letters),
data=test, vinit= c(1e-6, 1e-8, 1e-8))
*1。所有模型都可以工作,直到 fit4.
我发现你的模型令人担忧,因为它有一个随机斜率但没有随机截距,就像所有通过原点的回归让我担心的一样:我很难解释结果。尽管 lme 默认包含截距项,但 coxme 没有。
我曾希望 fit4 能奏效,也许有了更好的初始估计它就会奏效。 coxme 的底层代码是我在所有生存工作中遇到的最难的最大化问题,从最大化器很容易迷路并且永远找不到它的方式的意义上来说很难。这是一个有时需要手动操作的函数,通过有限的迭代计数 and/or 开始估计。我希望不是这样,我有一些长期计划通过添加一个替代的基于 MCMC 的最大化器来改进这一点,理论上它永远不会丢失但以更长的计算时间为代价。
如果任何方差太接近于零,则 sqrt() 消息往往会作为舍入误差的函数出现。当然,在您的测试用例中,实际 MLE 的方差为 0。发生这种情况时,我通常会通过使用一系列固定方差(vfixed 参数)进行拟合来直接检查零方差。如果可能性不变或随着方差值达到 1e-6 或更小而增加,那么我假设 MLE 为零并从模型中删除该随机项。 特里 T.*