INLA crashed with larger data, but worked in small data

[LC Chien]

Hi,

I tried to write a program for fitting a Cox PH model with my own data including over 8000 cancer patients.

The INLA crashed with the following message (it is very long, so I only take a short part of it):

        file: smtp-taucs.c  hgid: 9aeb760f40f0  date: Mon Oct 03 17:46:48 2011 +0200

        Function: GMRFLib_factorise_sparse_matrix_TAUCS(), Line: 771, Thread: 0

        Fail to factorize Q. I will try to fix it...

        file: smtp-taucs.c  hgid: 9aeb760f40f0  date: Mon Oct 03 17:46:48 2011 +0200

        Function: GMRFLib_factorise_sparse_matrix_TAUCS(), Line: 771, Thread: 0

        Fail to factorize Q. I will try to fix it...

        file: smtp-taucs.c  hgid: 9aeb760f40f0  date: Mon Oct 03 17:46:48 2011 +0200

        Function: GMRFLib_factorise_sparse_matrix_TAUCS(), Line: 771, Thread: 0

        Fail to factorize Q. I will try to fix it...

       This application has requested the Runtime to terminate it in an unusual way.

       Please contact the application's support team for more information.

However, when I reduced the sample size to 1000, it works. Hope some people can help this out. Thanks!!

(ps) My program looks like:

formula1=inla.surv(time,status1)~sex+white+black+f(agedx,model="rw2")

model1 = inla(formula1, family="coxph",

  data=crc, verbose=TRUE, quantiles = c(0.025, 0.1, 0.9, 0.975),

  control.hazard=list(model="rw1",n.intervals=20,param=c(1,0.001)))

[Finn Lindgren]

Can you please try adding the option
  control.inla(diagonal=100)
in your inla call (also try values other than 100, such as 10 and 10000) and see if that helps?  The option adjusts the Hessian in the optimisation in a way that may it more robust.

/Finn

[Daniel Simpson]

To further what Finn said, here is an example of this code being used in a different context (spatial extremes).  The first run essentially gets good starting values for the proper INLA method.  The options strategy="gaussian" and int.strategy="eb" will produce a cheap Gaussian approximation and marginalise things using empirical Bayes methods.  This is very fast.  (also - don't worry about control.predictor, it is needed for the specific model, but isnt' important to what I'm describing.)

starting.value = inla(formula,

            data=c(str.data$data, list(spde=spde_m)),

            family="gev",

            control.inla = list(diagonal = 100, strategy = "gaussian", int.strategy = "eb"),

            control.predictor=list(A=str.data$A),

            verbose=TRUE)

The starting value can then be used by calling the control.mode option.  This tells INLA where to start looking for the mode and restart=TRUE tells it to keep improving the approximation.

result = inla(formula,

            data=c(str.data$data, list(spde=spde_m)),

            family="gev",

            control.mode = list(result = starting.value, restart = TRUE),

            control.predictor=list(A=str.data$A),

            verbose=TRUE)

If this doesn't work, with diagonal =100, try it with a larger value and gradually "cool" it.  For example try the sequence diagonal = 1000, diagonal = 100, diagonal = 10.

[LC Chien]

I tried but the program still crashed. I am wondering whether the large sample size is really a concern. If this is a overall problem in any data set with large samples, that will be a huge limitation in INLA.  :(

[LC Chien]

I am sorry that I don't really understand the approach of obtaining a better starting value. I will learn it later. But if someone is interested in my data, I can share it with my program for people who can have a test to figure the problem out. I am so frustrated of using MCMC to estimate models due to its time-consumption. INLA is a big hope that I can save my time to do more research. I hope the problem can be eventually solved.

[Finn Lindgren]

The data size by itself is not the issue.  I run models with almost 125000 data points, with a model graph of about 850000 without running into this problem. (I do need lots of memory though).  Did you also try the more detailed suggestions from Daniel?

Can you tell me the internal size of the latent model graph (stated just before the start of the optimisation when using the verbose=TRUE option to inla)?

/Finn

[LC Chien]

Here I am updating the progress of this question. Now it has been solved. Thanks for Finn's great advice!

Just as Finn mentioned, this program should include a control.inla option and modify the number of diagonal. It looks like:

crc <- read.table("detroit_inla.dat", header=T, sep="")

crc <- crc[order(crc$time),]

## Conventional Cox PH model ##

formula1=inla.surv(time,status1)~sex+white+black+f(agedx,model="rw2")

## Geosurvival model ##

formula2=inla.surv(time,status1)~sex+white+black+f(agedx,model="rw2")+f(tract2,model="besag",graph.file="c:/lchien/crc1/data/map/detroit_inla.gra")

model1 = inla(formula1, family="coxph",

                        data=crc, verbose=TRUE, quantiles = c(0.025, 0.1, 0.9, 0.975),

                        control.inla=list(diagonal=10000),

                        control.hazard=list(model="rw1",n.intervals=20,param=c(1,0.001)))

model2 = inla(formula2, family="coxph",

                        data=crc, verbose=TRUE, quantiles = c(0.025, 0.1, 0.9, 0.975),

                        control.inla=list(diagonal=10000),

                        control.hazard=list(model="rw1",n.intervals=20,param=c(1,0.001)))

summary(model1)

summary(model2)

They spent 86 seconds and 216 seconds, respectively! I also tried diagonal = 10 and diagonal = 1000, both of them worked as well.

An additional update is that a problem appeared in the spatial function of the geosurvival model. I found the graph file doesn't work while using 0 to be the first code of the location. I used to read a statement in one of the tutorial documentations which mentioned that the graph file also works with starting location code = 0. So far I haven't figure out why it doesn't work in my graph file.  After writing an additional SAS program to convert the graph file from starting location = 0 to starting location = 1, the INLA() finally ran well. 

Thanks all advice in these days!

[Håvard Rue]

On Thu, 2011-11-03 at 08:53 -0700, LC Chien wrote:

>
> model1 = inla(formula1, family="coxph",
> data=crc, verbose=TRUE, quantiles = c(0.025,
> 0.1, 0.9, 0.975),
> control.inla=list(diagonal=10000),
>
> control.hazard=list(model="rw1",n.intervals=20,param=c(1,0.001)))

I just want to add the not that the 'diagonal=10000' option really
*change* the model, as its add the constant 'diagonal' to all the
diagonal elements of the precision matrix for the latent field. So its
ment as an emergency option really...

it means that, in order to get the the correct results, you'll need to
rerun the model with diagonal=0, starting from the (better) initial
values found as the mode in the first run. so the correct call is

model1 = inla(formula1, family="coxph",
data=crc, verbose=TRUE, quantiles = c(0.025, 0.1,
0.9, 0.975),
control.inla=list(diagonal=10000),

## these options make it runs faster as we're
only interested
## in the mode
control.inla = list(int.strategy = "eb",
strategy = "gaussian"),
##

control.hazard=list(model="rw1",n.intervals=20,param=c(1,0.001)))

## and then we restart it... adding the initial values as the mode of
'model1'.

model1 = inla(formula1, family="coxph",
data=crc, verbose=TRUE, quantiles = c(0.025,
0.1, 0.9, 0.975),

control.inla=list(diagonal=0),
##
## saying, use the mode in model1, and restart
the optimization
control.mode = list( result = model1, restart =
TRUE),
##

control.hazard=list(model="rw1",n.intervals=20,param=c(1,0.001)))

>
> An additional update is that a problem appeared in the spatial
> function of the geosurvival model. I found the graph file doesn't work
> while using 0 to be the first code of the location. I used to read a
> statement in one of the tutorial documentations which mentioned that
> the graph file also works with starting location code = 0. So far I
> haven't figure out why it doesn't work in my graph file. After
> writing an additional SAS program to convert the graph file from
> starting location = 0 to starting location = 1, the INLA() finally ran
> well.

the story about the 0-based graphs: this feature is a left-over
from old times. a graph can be 0-based but the R-code must be
1-based!

to convert a 0-based graph into a 1-based graph, you can do:

inla.write.graph( inla.read.graph("old.graph"), "new.graph")

H

--
Håvard Rue
Department of Mathematical Sciences
Norwegian University of Science and Technology
N-7491 Trondheim, Norway
Voice: +47-7359-3533 URL : http://www.math.ntnu.no/~hrue
Fax : +47-7359-3524 Email: havar...@math.ntnu.no

This message was created in a Microsoft-free computing environment.

[LC Chien]

I was still testing the R codes that you suggested. I understand Håvard's consideration because the result is much different from the result generated by SAS PROC PHREG procedure. So it is really need to rerun the model with the initial values. However, the INLA() crashed again when reruning the program with control.mode option. The error message is the same showing in my original post. I attached my updated program and data, and hope experts here can point out my mistake in the program. I pretty appreciate any advanced help.

[LC Chien]

I think I solved this problem by replacing control.inla=list(diagonal=0) to control.inla=list(diagonal=10), and all estimates were reasonable. However, when I used control.inla=list(diagonal=100) to rerun the model, the estimates became not reasonable again. I don't know why, but I will use control.inla=list(diagonal=10) temporarily.

[LC Chien]

Here is the updated progress of this post. The problem has been solved by including a proper prior for the intercept, so now there is no need to try different initial values and restart the model-fitting.

model1 = inla(formula1, family="coxph",
              data=crc, verbose=TRUE, quantiles = c(0.025,0.1, 0.9, 0.975),

              control.inla=list(diagonal=0),

              control.fixed = list(prec.intercept = 0.1),
              control.hazard=list(model="rw1",n.intervals=20,param=c(1,0.001)))

Thanks for Håvard's guidance.

[liyuji...@gmail.com]

Hi, Could I have your dataset so that i could try this code?  Thanks sincerely!

在 2011年11月3日星期四UTC+1下午4时53分39秒，LC Chien写道：