inlabru-org / inlabru

inlabru
https://inlabru-org.github.io/inlabru/
73 stars 21 forks source link

Problems with the group specification #221

Open JoseCansec opened 6 months ago

JoseCansec commented 6 months ago

Dear Finn, we are currently trying to fit a separable spatio-temporal model which includes a grouping variable (id cell) within the temporal model (hour) specification. The main problem is that the amount of time it takes to fit the model, sometimes giving a inla.core crashed type of error. We are not sure if we are specyfing the control.group vairbable correctly as we want to include the repeated measurments within a cell (xy) during each hour. I share our model specification bellow:

Spatial model

cmp.d.sth2 = corden ~ 1 + myspde(main = coordinates, model = Hosts.spde.d)+ f2(hour.2,model="ou",n=n.hour,group=id.cell,ngroup=n.cell.pos,control.group =list(model="iid"))

Fit de model

bru.d.sth2 = bru(cmp.d.sth2, family = "gamma", data=pos2.0.sp, options=list(num.threads="10", control.inla=list(int.strategy="eb", strategy="gaussian"), control.compute = list(cpo = TRUE, dic = TRUE, waic = TRUE)))

Thank you in advance for the help and we look forward to your answer. Jose

finnlindgren commented 6 months ago

Hard to say without more information. What's the data structure, e.g. head(pos2.0.sp) and str(pos2.0.sp)? Also, try

bru.d.sth2 = bru(cmp.d.sth2, family = "gamma", data=pos2.0.sp,
  options=list(num.threads="10", control.inla=list(int.strategy="eb", strategy="gaussian"),
    control.compute = list(cpo = TRUE, dic = TRUE, waic = TRUE),
    bru_run = FALSE))
summary(bru.d.sth2, verbose = TRUE)

since that might show me some useful information as well.

JoseCansec commented 6 months ago

Dear Finn

The structure of my data and the first five rows are these:

head(pos.15m);str(pos.15m) id.cell xcell ycell hour corden hour.2 1 439 3298 439 3298 1 1101.3370 1 2 440 3298 440 3298 1 871.2178 1 3 440 3299 440 3299 1 3204.5860 1 4 441 3299 441 3299 1 1684.4960 1 5 441 3299 441 3299 2 1021.3860 2 6 441 3300 441 3300 2 1757.6060 2 'data.frame': 629 obs. of 6 variables: $ id.cell: Factor w/ 187 levels "420 3303","420 3304",..: 116 124 125 132 132 133 134 135 143 133 ... $ xcell : num 439 440 440 441 441 441 441 441 442 441 ... $ ycell : num 3298 3298 3299 3299 3299 ... $ hour : num 1 1 1 1 2 2 2 2 3 3 ... $ corden : num 1101 871 3205 1684 1021 ... $ hour.2 : num 1 1 1 1 2 2 2 2 3 3 ...

The summary of one run (we are doing 1000 simulations) with verbose=TRUE is the following:

summary(bru.d.sth2, verbose=TRUE)inlabru version: 2.10.0 INLA version: 23.04.24 Components: Label: myspde Type: main = spde, group = exchangeable, replicate = iid Input: main = coordinates, group = 1L, replicate = 1L Map: pipe = multi(main = fm_mesh_2d, group = index, replicate = index) -> scale INLA formula: ~ . + f(myspde, model = BRU_myspde_main_model, ngroup = 1, nrep = 1, values = BRU_myspde_values) Label: f2 Type: main = ou, group = iid, replicate = iid Input: main = hour.2, group = id.cell, replicate = 1L Map: pipe = multi(main = fm_mesh_1d, group = index, replicate = index) -> scale INLA formula: ~ . + f(f2, model = BRU_f2_main_model, n = 168, group = f2.group, ngroup = 187, control.group = list(model = "iid"), nrep = 1, values = BRU_f2_values) Label: Intercept Type: main = linear, group = exchangeable, replicate = iid Input: main = 1, group = 1L, replicate = 1L Map: pipe = multi(main = linear, group = index, replicate = index) -> scale INLA formula: ~ . + f(Intercept, model = BRU_Intercept_main_model, ngroup = 1, nrep = 1, values = BRU_Intercept_values) Likelihoods: Family: 'gamma' Data class: 'SpatialPointsDataFrame' Predictor: corden ~ . Time used: Pre = 5.53, Running = 374, Post = 1.14, Total = 381 Fixed effects: mean sd 0.025quant 0.5quant 0.975quant mode kld Intercept 7.285 0.393 6.514 7.285 8.056 7.285 0

Random effects: Name Model myspde SPDE2 model f2 Ornstein-Uhlenbeck model

Model hyperparameters: mean sd 0.025quant 0.5quant 0.975quant mode Precision parameter for the Gamma observations 9.20 5.46e-01 8.091 9.214 10.23 9.307 Range for myspde 45.39 3.22e+01 14.307 36.603 133.00 24.765 Stdev for myspde 0.78 4.48e-01 0.299 0.665 1.99 0.487 Precision for f2 10891.82 1.55e+04 438.171 5997.062 51320.93 1075.249 Phi for f2 32.09 7.30e+02 0.027 1.328 178.78 0.032

Deviance Information Criterion (DIC) ...............: 9820.53 Deviance Information Criterion (DIC, saturated) ....: 676.40 Effective number of parameters .....................: 22.11

Watanabe-Akaike information criterion (WAIC) ...: 9824.56 Effective number of parameters .................: 24.58

Marginal log-Likelihood: -4936.38 CPO, PIT is computed Posterior summaries for the linear predictor and the fitted values are computed

(Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

We would like to know if the id.cell group specification is correct given that we need to consider the repeated measurments in each cell for each hour.

Thanks in advance.

Jose

El jue, 21 dic 2023 a las 10:27, Finn Lindgren @.***>) escribió:

Hard to say without more information. What's the data structure, e.g. head(pos2.0.sp) and str(pos2.0.sp)? Also, try

bru.d.sth2 = bru(cmp.d.sth2, family = "gamma", data=pos2.0.sp, options=list(num.threads="10", control.inla=list(int.strategy="eb", strategy="gaussian"), control.compute = list(cpo = TRUE, dic = TRUE, waic = TRUE), bru_run = FALSE)) summary(bru.d.sth2, verbose = TRUE)

since that might show me some useful information as well.

— Reply to this email directly, view it on GitHub https://github.com/inlabru-org/inlabru/issues/221#issuecomment-1865930401, or unsubscribe https://github.com/notifications/unsubscribe-auth/ARDE42V2CG57XUPPJYQSBRDYKP6G5AVCNFSM6AAAAABA4VMIMGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQNRVHEZTANBQGE . You are receiving this because you authored the thread.Message ID: @.***>

finnlindgren commented 5 months ago

Hi Jose, I'm now "back in action" after holidays and exam marking. In your info above, you gave a different data variable, pos.15m, than the one used in the actual estimation, pos2.0.sp. Since the code uses sp features, I'll need to see the pos2.0.sp information, as pos.15m doesn't have any geometry information. Also print(Hosts.spde.d$mesh) would be useful. The f2 component has around 31k elements, but I couldn't see the spde model size.