Issue specifying link for prediction -- nbinomial distributional assumption

18 views
Skip to first unread message

Nicholas Arisco

unread,
Jul 16, 2024, 2:29:42 PM (5 days ago) Jul 16
to R-inla discussion group
Hello,

I am working on a spatiotemporal model that uses a set of scenario values to predict a subset of years of the outcome data, with a negative binomial distributional assumption. The outcome to be predicted is set to NA for the years of interest. The model runs fine, but the summary.linear.predictor values for the prediction period are not being transformed from the summary.fitted.values output. The other outcome data are being correctly transformed.

My model code is as follow:

m0 <- cases ~ 1 +
  scale(area_mined / area) +
  scale(month_control_cosine) +
  scale(defrst_minus_mining_km2 / area) +
  scale(population) +
  scale(area) +
  f(nome.2, model='bym', graph=mun.adj, scale.model=T, param=c(0.001,0.001)) +
  f(nome.1, model='iid') +
  f(yn, model='iid') +
  f(ymn, model='iid') +
  f(inla.group(total_precipitation_sum), model='ar1') +
  f(inla.group(temperature_2m_max), model='ar', order=2) +
  f(ONI, model='ar', order=2)

m0_resultspre2017 <- inla(m0, family="nbinomial", data=d1_full3,
                          control.predictor(list(compute=TRUE,link=1)),
                          control.compute(list(dic=TRUE, waic = TRUE, cpo = TRUE, config = TRUE)),
                          control.inla(list(tolerance=0.01)),
                          verbose = T, safe = T)

Any advice would be greatly appreciated.

Sincerely,
Nick

Helpdesk (Haavard Rue)

unread,
Jul 17, 2024, 12:43:45 PM (4 days ago) Jul 17
to Nicholas Arisco, R-inla discussion group
I think there are several issues with the code below.

On Tue, 2024-07-16 at 11:25 -0700, Nicholas Arisco wrote:
> Hello,
>
> I am working on a spatiotemporal model that uses a set of scenario values to
> predict a subset of years of the outcome data, with a negative binomial
> distributional assumption. The outcome to be predicted is set to NA for the
> years of interest. The model runs fine, but the summary.linear.predictor
> values for the prediction period are not being transformed from the
> summary.fitted.values output. The other outcome data are being correctly
> transformed.
>
> My model code is as follow:
>
> m0 <- cases ~ 1 +
>   scale(area_mined / area) +
>   scale(month_control_cosine) +
>   scale(defrst_minus_mining_km2 / area) +
>   scale(population) +
>   scale(area) +


>   f(nome.2, model='bym', graph=mun.adj, scale.model=T, param=c(0.001,0.001)) +

There is better model than 'bym', called 'bym2' that you might use instead. see
the doc and the relevant publication



>   f(nome.1, model='iid') +
>   f(yn, model='iid') +
>   f(ymn, model='iid') +
>   f(inla.group(total_precipitation_sum), model='ar1') +
>   f(inla.group(temperature_2m_max), model='ar', order=2) +
>   f(ONI, model='ar', order=2)
>
> m0_resultspre2017 <- inla(m0, family="nbinomial", data=d1_full3,
>                           control.predictor(list(compute=TRUE,link=1)),
>                           control.compute(list(dic=TRUE, waic = TRUE, cpo =
> TRUE, config = TRUE)),
>                           control.inla(list(tolerance=0.01)),
>                           verbose = T, safe = T)
>
> Any advice would be greatly appreciated.
>
> Sincerely,
> Nick
> --
> You received this message because you are subscribed to the Google Groups "R-
> inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to r-inla-discussion...@googlegroups.com.
> To view this discussion on the web, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/d200e416-66f1-4374-9190-c97246809ffdn%40googlegroups.com
> .

--
Håvard Rue
he...@r-inla.org

Helpdesk (Haavard Rue)

unread,
Jul 17, 2024, 1:06:06 PM (4 days ago) Jul 17
to Nicholas Arisco, R-inla discussion group
[that was sent to early... I'll try again]

I think there are several issues with the code below.

On Tue, 2024-07-16 at 11:25 -0700, Nicholas Arisco wrote:
> Hello,
>
> I am working on a spatiotemporal model that uses a set of scenario values to
> predict a subset of years of the outcome data, with a negative binomial
> distributional assumption. The outcome to be predicted is set to NA for the
> years of interest. The model runs fine, but the summary.linear.predictor
> values for the prediction period are not being transformed from the
> summary.fitted.values output. The other outcome data are being correctly
> transformed.
>
> My model code is as follow:
>
> m0 <- cases ~ 1 +
>   scale(area_mined / area) +
>   scale(month_control_cosine) +
>   scale(defrst_minus_mining_km2 / area) +
>   scale(population) +
>   scale(area) +


>   f(nome.2, model='bym', graph=mun.adj, scale.model=T, param=c(0.001,0.001)) +

There is better model than 'bym', called 'bym2' that you might use instead. see
the doc and the relevant publication

https://journals.sagepub.com/doi/10.1177/0962280216660421

>   f(inla.group(total_precipitation_sum), model='ar1') +

the arguments to 'ar1', should be integers, 1, 2,...,n, but this is not
enforeced. so if

f(idx, model="ar1")

and idx=c(1.1, 3.3, 4.4), these are mapped into time=1,2,3 even though the
distance inbetween is not the same. so the use of 'inla.group()' is somewhat
questionable for the ar1 and similar for 'ar'.

see the attached runme3.R
>


> m0_resultspre2017 <- inla(m0, family="nbinomial", data=d1_full3,
>                           control.predictor(list(compute=TRUE,link=1)),
>                           control.compute(list(dic=TRUE, waic = TRUE, cpo =
> TRUE, config = TRUE)),
>                           control.inla(list(tolerance=0.01)),
>                           verbose = T, safe = T)


the 'control.compute()' function is not to be confused with the
'control.compute' argument. its ment to auto-completion when writing, not for
use in code, you can replace it (and the others control.xxx similarly), as

control.compute = control.compute(dic=TRUE, waic = TRUE, cpo = TRUE, config =
TRUE)

or simply

control.compute = list(dic=TRUE, waic = TRUE, cpo = TRUE, config = TRUE)




please retry, and let me know if there are still issues

best
H




--
Håvard Rue
he...@r-inla.org
runme3.R

Arisco, Nicholas

unread,
Jul 17, 2024, 4:15:36 PM (4 days ago) Jul 17
to Helpdesk, R-inla discussion group
Thank you, Havard, this is very helpful. You have successfully fixed my issue!
Reply all
Reply to author
Forward
0 new messages