search

mclements / rstpm2

An R package for generalised survival models
28 stars 11 forks source link

Cannot estimate models with frailty terms #5

Closed ellessenne closed 7 years ago

ellessenne commented 7 years ago

Dear Mark, I have been trying to use rstpm2 to estimate models with shared frailty terms:

library(rstpm2)
data(brcancer)
fit <- stpm2(Surv(rectime, censrec == 1) ~ hormon, df = 5, data = brcancer, frailty = TRUE, cluster = "x4", RandDist = "LogN")

The estimation fails, as it gets to a bad Hessian that cannot be inverted:

Error in solve.default(hessian) : 
Lapack routine dgesv: system is exactly singular: U[2,2] = 0

I get the same error if I use a Gamma-distributed frailty. I tried comparing with Stata's stmixed:

webuse brcancer
stset rectime censrec
stmixed hormon || x4: , dist(fpm) df(5)

and this worked. Any insight into this, or tip for solving it? I tried:

My session info:

> sessionInfo()
R version 3.3.3 (2017-03-06)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=English_United Kingdom.1252 
[2] LC_CTYPE=English_United Kingdom.1252   
[3] LC_MONETARY=English_United Kingdom.1252
[4] LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.1252    

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
[1] rstpm2_1.3.4    survival_2.41-3

loaded via a namespace (and not attached):
 [1] Matrix_1.2-8      tools_3.3.3       mgcv_1.8-17       Rcpp_0.12.10     
 [5] bbmle_1.0.18      fastGHQuad_0.2    nlme_3.1-131      grid_3.3.3       
 [9] numDeriv_2016.8-1 stats4_3.3.3      lattice_0.20-35

Many thanks in advance for your feedback! Best wishes,

Alessandro

mclements commented 7 years ago

Thank you for this report.

brcancer$x4 has only 3 levels and would be more appropriately fitted as a fixed effect. I am surprised that stmixed works: how did you interpret the answer that it gives?

Does this address your question? How's your time in Leicester?

Sincerely, Mark.

On 04/12/2017 11:57 AM, Alessandro Gasparini wrote:

Dear Mark, I have been trying to use rstpm2 to estimate models with shared frailty terms:

library(rstpm2) data(brcancer) fit <- stpm2(Surv(rectime, censrec == 1) ~ hormon, df = 5, data = brcancer, frailty = TRUE, cluster = "x4", RandDist = "LogN")

The estimation fails, as it gets to a bad Hessian that cannot be inverted:

Error in solve.default(hessian) : Lapack routine dgesv: system is exactly singular: U[2,2] = 0

I get the same error if I use a Gamma-distributed frailty. I tried comparing with Stata's stmixed:

webuse brcancer stset rectime censrec stmixed hormon || x4: , dist(fpm) df(5)

and this worked. Any insight into this, or tip for solving it? I tried:

My session info:

sessionInfo() R version 3.3.3 (2017-03-06) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 7 x64 (build 7601) Service Pack 1

locale: [1] LC_COLLATE=English_United Kingdom.1252 [2] LC_CTYPE=English_United Kingdom.1252 [3] LC_MONETARY=English_United Kingdom.1252 [4] LC_NUMERIC=C [5] LC_TIME=English_United Kingdom.1252

attached base packages: [1] splines stats graphics grDevices utils datasets methods [8] base

other attached packages: [1] rstpm2_1.3.4 survival_2.41-3

loaded via a namespace (and not attached): [1] Matrix_1.2-8 tools_3.3.3 mgcv_1.8-17 Rcpp_0.12.10 [5] bbmle_1.0.18 fastGHQuad_0.2 nlme_3.1-131 grid_3.3.3 [9] numDeriv_2016.8-1 stats4_3.3.3 lattice_0.20-35

Many thanks in advance for your feedback! Best wishes,

Alessandro

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/mclements/rstpm2/issues/5, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AAiBCRH6wpJeywiNDQNmykk3p3R2QU0Kks5rvKAcgaJpZM4M7NF-.

ellessenne commented 7 years ago

Many thanks for your feedback Mark - also, everything is going well here in Leicester!

Anyway, the breast cancer dataset was more of an example than a practical application: I am currently simulating data under different scenarios to compare the performance of a set of models - including Royston-Parmar models with frailties. I tried also with the catheter dataset (data(kidney, package = "survival")), which includes repeated events (infections) for each individual, using an individual-level frailty - and had no luck there either. Any way I could dig into it or give you more information about what triggers the error? Best wishes,

Alessandro

xinliu commented 7 years ago

Hi Alessandro,

Please try it again.

data(brcancer, package = "rstpm2") summary(stpm2(Surv(rectime, censrec == 1) ~ hormon, df = 5, data = brcancer, cluster = brcancer$x4, RandDist = "Gamma"))

data(kidney, package = "survival") summary(stpm2(Surv(time, status == 1) ~ ns(age, df = 4) + sex, df = 4, cluster = kidney$id, data = kidney, RandDist = "LogN"))

summary(stpm2(Surv(time, status == 1) ~ age + sex, df = 4, cluster = kidney$id, data = kidney, RandDist = "LogN"))

summary(stpm2(Surv(time, status == 1) ~ age + sex, df = 4, cluster = kidney$id, data = kidney, RandDist = "Gamma"))

Xing-Rong

ellessenne commented 7 years ago

That's great Xing-Rong, thanks a lot! I would kindly suggest adding your first example to the documentation for the stpm2 function, it could be useful to avoid confusion as if I run:

library(rstpm2)
data(brcancer, package = "rstpm2")
summary(stpm2(Surv(rectime, censrec == 1) ~ hormon, df = 5, data = brcancer, cluster = x4, RandDist = "Gamma"))

it doesn't work anymore (cannot find x4 in the current environment). I was expecting the value passed to cluster to be evaluated within the environment defined by the data argument, or it being a string processed using standard evaluation - hence my confusion. Anyway, thanks a lot for clarifying this! Best wishes,

Alessandro