seasonal latent effects

185 views
Skip to first unread message

Rustam Mussa

unread,
Mar 5, 2022, 8:53:21 AM3/5/22
to R-inla discussion group

Dear INLA experts,

 

For seasonal latent effects, the periodicity of the seasonal effect is set with argument season.length. When you do this, does the periodicity depend on a unit of observation (data points)? For example, if I have a variable “date” which is measured year by year, as 2000, 2001, 2002… 2020, and I want to assume that the seasonal effect covers one year, then should I go as f(date, model=”seasonal”, season.length=1)? And should it be season.length=5, if I decide that the periodicity is 5 years?

 

Sincerely,

Rustam 

Helpdesk

unread,
Mar 5, 2022, 1:53:09 PM3/5/22
to Rustam Mussa, R-inla discussion group
you can check the example in

inla.doc("seasonal")

the period is related to the indexing set, which is easiest defined as
integers, like 1,2,3...

the season.length=5 is that is the cycle length.

The easiest way to check that you got it right, it to simulate some easy
to estimate data, and check that you get it back and all is looking
good.

see also the GMRF-book at the end of chapter 3.




-----Original Message-----
From: Rustam Mussa <rmuss...@gmail.com>
To: R-inla discussion group <r-inla-disc...@googlegroups.com>
Subject: [r-inla] seasonal latent effects
Date: Sat, 5 Mar 2022 05:53:21 -0800 (PST)

Dear INLA experts,
 
For seasonal latent effects, the periodicity of the seasonal effect is
set with argument season.length.When you do this, does the periodicity
depend on a unit of observation (data points)? For example, if I have a
variable “date” which is measured year by year, as 2000, 2001, 2002…
2020, and I want to assume that the seasonal effect covers one year,
then should I go as f(date, model=”seasonal”, season.length=1)? And
should it be season.length=5, if I decide that the periodicity is 5
years?
 
Sincerely,
Rustam 
--
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/3bf3929f-ef25-4350-81e3-98358533bff3n%40googlegroups.com.

Angela Fanelli

unread,
Jan 20, 2025, 5:51:17 AM1/20/25
to R-inla discussion group
Hi Rustam, 
Thank you for your answer. 
I created a reproducible code to make sure I implement it properly in INLA. 

Håvard Rue suggested using NA observations formulation. 
However I still do not understand how to check and assess the effect of seasonality. 
I prefer not to switch to a rw1 model. 

Also once I make predictions I am not sure if I have to convert the posterior fitted distribution from the logistic scale to the scale of the data using a custom function. I've been using the following function:

myfun <- function(x) { exp(x) / (1 + exp(x)) }

I saw that people do a transformation when predicting using  lincomb  but I am not sure if it is necessary if I extract directly the predictions .. they look fine without a transformation.

Best wishes 
ReproducibleModel.R

Helpdesk (Haavard Rue)

unread,
Jan 20, 2025, 4:03:00 PM1/20/25
to Angela Fanelli, R-inla discussion group
for a NA observation, you need to set the argument 'link' to 'i' to say that
this NA is related to the i'th likelihood. Otherwise, you can map it yourself,
using the function inla.tmarginal and inla.zmarginal


n = 100
z = runif(n)
eta = 1 + 0.1*z
N = 10


## logit
p = exp(eta)/(1+exp(eta))
y = rbinom(n, size = N, prob = p)
y[(n-10):n] = NA
r = inla(y ~ 1 + z, data = data.frame(y, z), family = "binomial", Ntrials =
rep(N, n),
control.family = list(control.link = list(model = "logit")),
control.predictor = list(compute=TRUE))

rr = inla(y ~ 1 + z, data = data.frame(y, z), family = "binomial", Ntrials =
rep(N, n),
control.family = list(control.link = list(model = "logit")),
control.predictor = list(compute=TRUE, link=1))

rrr = inla(y ~ 1 + z, data = data.frame(y, z), family = "binomial", Ntrials =
rep(N, n),
control.family = list(control.link = list(model = "logit")),
control.predictor = list(compute=TRUE, link=NA))

tail(cbind(r$summary.fitted.values$mean, rr$summary.fitted.values$mean,
rrr$summary.fitted.values$mean))
Håvard Rue
he...@r-inla.org
Reply all
Reply to author
Forward
0 new messages