Forecasting observations (t+1) using spatio-temporal areal unit modelling
564 views
Skip to first unread message
Thomas
Jul 20, 2018, 3:57:01 PM7/20/18
to R-inla discussion group
Hi,
I am after some help with forecasting, month+1 after fitting a spatio-temporal areal unit model using the formula and data below;
formula <- y ~ 1 +
f(area, model="bym", graph =adj) +
f(month,model='ar1') +
f(month1,model='iid') +
mod <- inla(formula,family="poisson",data=data, control.compute = list(dic = TRUE)))
############################################################
* not true data
| y | area | month | month1 | area.int | month.int | fitted |
| 142 | 1 | 1 | 1 | 1 | 1 | 142.6 |
| 128 | 2 | 1 | 1 | 2 | 1 | 129.7 |
| 268 | 3 | 1 | 1 | 3 | 1 | 268.2 |
| 156 | 4 | 1 | 1 | 4 | 1 | 157.0 |
| … | … | … | … | ... | ... | |
| 46 | 1 | 2 | 2 | 1 | 2 | 48.3 |
| 287 | 2 | 2 | 2 | 2 | 2 | 286.5 |
| 125 |
3 | 2 | 2 | 3 | 2 | 126.2 |
| 408 | 4 | 2 | 2 | 4 | 2 | 407.6 |
###########################################################
As you can see, the fitted values from the model are close to the true observations.
Random Effect Output:
| Hyperparameters | mean | sd | 0.025quant | 0.5quant | 0.975quant | mode |
| Precision for area (iid) |
0.88 | 0.05 | 0.78 | 0.88 | 0.98 | 0.88 |
| Precision for area (spatial component) | 1718.00 | 1819.00 | 122.99 | 1173.00 | 6580.00 | 338.33 |
| Precision for month | 21.54 | 8.42 | 9.06 | 20.27 | 41.73 | 17.75 |
| Rho for month | 0.62 | 0.14 | 0.32 | 0.63 | 0.85 | 0.66 |
| Precision for month1 | 17990.00 | 18310.00 | 1246.90 | 12550.00 | 66570.00 | 3432.81 |
| Precision for area.int | 2.08 | 0.04 | 2.00 | 2.08 | 2.17 | 2.08 |
| GroupRho for areaint | 0.08 | 0.02 | 0.05 | 0.08 | 0.11 | 0.09 |
###########################################################
For each area, I have counts from the month 01-12 of 2017. So I want to predict the counts for each area for 01/18.
Now I can't reveal the true data for confidentiality reasons, all I can say is that the model is at a regional scale, so covers a large area.
I have tried following the FAQ - How can I do predictions using INLA? Which I tested on for month 12, by simply setting y[i] = NA for all the areas in month 12 I want to predict NOT just the locations I want to predict. This gave values far from the actual values of all areas for month 12, in the region of y=5 (mean).
Please could some one explain how I can predict counts for month 12 in space&time.
Finn Lindgren
Jul 20, 2018, 4:11:26 PM7/20/18
to Thomas, R-inla discussion group
Did you set the "link" parameter which tells inla to apply the link function also to NA "response" predictions in fitted.values?
See
(Important: the posterior predictions are for the model parameter controlling the mean, not the counts themselves; predictive count distributions require an extra step, typically using inla.posterior.sample; inlabru::predict can do it semi-automatically, but inlabru doesn't quite support all models yet)
Finn
--
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-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-discussion-group@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
--
Finn Lindgren
Jul 20, 2018, 4:13:04 PM7/20/18
to Thomas, R-inla discussion group
Ah, just noticed that you did mention the FAQ. But since you didn't mention link (link=1 should be sufficient in your case), please double check that.
Finn
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion-group+unsubscr...@googlegroups.com.
To post to this group, send email to r-inla-discussion-group@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
--
Thomas
Jul 20, 2018, 4:25:21 PM7/20/18
to R-inla discussion group
Hi Finn,
Firstly thank you for such a quick response!
No I did not set the link function, that was just thinking about that... so i'll try link=1, is that the correct assumption.
Also thank you for directing me in the right direction. I have come across inla.posterior.sample() and from my understanding, control.compute( config=TRUE) needs to be called in the INLA call.
Would 1000 samples be enough for 7896 observations? 658(areas)*12(months)?
Please could you direct me to maybe a link or some sample code in how to predict month+1 using the inla.posterior.sample() function, as there seems to not be a lot of examples online.
Thank you for your support, it's very appreciated,
Thomas.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
Finn Lindgren
Jul 20, 2018, 4:34:23 PM7/20/18
to Thomas, R-inla discussion group
On 20 July 2018 at 21:25, Thomas <tast...@gmail.com> wrote:
Hi Finn,Firstly thank you for such a quick response!No I did not set the link function, that was just thinking about that... so i'll try link=1, is that the correct assumption.
Correct. That means "apply the link function from the first likelihood to all the elements", and since you only have one likelihood, that should do it.
Also thank you for directing me in the right direction. I have come across inla.posterior.sample() and from my understanding, control.compute( config=TRUE) needs to be called in the INLA call.
Correct.
Would 1000 samples be enough for 7896 observations? 658(areas)*12(months)?
The number of posterior samples required has very little to do with the amount of input data; inla.posterior.sample generates independent samples, so the number of samples depends mostly on how many samples are needed to estimate the desired quantities with the desired precision (e.g. 10 is too small to reliably estimate a variance, but 1000 is quite decent).
Please could you direct me to maybe a link or some sample code in how to predict month+1 using the inla.posterior.sample() function, as there seems to not be a lot of examples online.
I don't have any links handy, sorry; hopefully someone else will chime in. But in your case, your model should really be set up in the way you already have it, including month+1 in the model.
It's then "just" a matter of extracting the relevant part of each posterior sample, and compute whatever quantity you're interested in. For example, to compute P(N<=n | y),
P(N <= n | y) = E( P(N <= n | \eta, y) ) \approx 1/S \sum_{s=1}^S P(N <= n | \eta^{[s]}, y)
where \eta is the linear predictor, and \eta^{[s]} is each sample of \eta.
Another option is to sample from the counts.
Finn
--
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-group+unsubscr...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
--
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-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-discussion-group@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
Thomas
Jul 20, 2018, 4:52:22 PM7/20/18
to R-inla discussion group
Hi Finn,
Thanks for your quick response again and clarifying a few things for me.
Ok, now I understand the required data structure and the theoretical underpinning of what's required for forecasting month+1 using a spatio-temporal areal unit model.
Yes, hopefully someone very kind will chirp in with some sample code they already have / can contribute to.
Regards,
Tom.
On Friday, 20 July 2018 21:34:23 UTC+1, Finn Lindgren wrote:
On 20 July 2018 at 21:25, Thomas <tast...@gmail.com> wrote:Hi Finn,Firstly thank you for such a quick response!No I did not set the link function, that was just thinking about that... so i'll try link=1, is that the correct assumption.Correct. That means "apply the link function from the first likelihood to all the elements", and since you only have one likelihood, that should do it.Also thank you for directing me in the right direction. I have come across inla.posterior.sample() and from my understanding, control.compute( config=TRUE) needs to be called in the INLA call.Correct.Would 1000 samples be enough for 7896 observations? 658(areas)*12(months)?The number of posterior samples required has very little to do with the amount of input data; inla.posterior.sample generates independent samples, so the number of samples depends mostly on how many samples are needed to estimate the desired quantities with the desired precision (e.g. 10 is too small to reliably estimate a variance, but 1000 is quite decent).Please could you direct me to maybe a link or some sample code in how to predict month+1 using the inla.posterior.sample() function, as there seems to not be a lot of examples online.I don't have any links handy, sorry; hopefully someone else will chime in. But in your case, your model should really be set up in the way you already have it, including month+1 in the model.It's then "just" a matter of extracting the relevant part of each posterior sample, and compute whatever quantity you're interested in. For example, to compute P(N<=n | y),P(N <= n | y) = E( P(N <= n | \eta, y) ) \approx 1/S \sum_{s=1}^S P(N <= n | \eta^{[s]}, y)where \eta is the linear predictor, and \eta^{[s]} is each sample of \eta.Another option is to sample from the counts.Finn
--
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-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
--
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-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
Haakon Bakka
Jul 25, 2018, 6:26:17 AM7/25/18
to Thomas, R-inla discussion group
Hi,
The only code I have for posterior sampling is explained here:
I suggest you add a fake NA-observation year for the year that you want to predict.
I assume that someone has done something similar to what you want in the past, so maybe someone has more relevant code.
Kind regards,
Haakon
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion-group+unsubscr...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
--
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-group+unsubscr...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
--
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-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-discussion-group@googlegroups.com.
Thomas
Jul 25, 2018, 4:52:57 PM7/25/18
to R-inla discussion group
Hi Haakon,
Firstly thank you for your reply, I found your website very useful and it certainly makes the use of the inla.posterior.sample() meaningful in terms of seeing the structure of what this function does!
# define n samplesn.samples = 1000# generate 1000 posterior samples from all of its components from the INLA modelsamples = inla.posterior.sample(n.samples, result = train)
# create a generic code to see the hidden part of the INLA result object & start contents = train$misc$configs$contentscontentsPreviously I have added fake NA observations for the months that I want to predict and simply "used control.predictor = list(link = 1)))". This gives an output but not as accurate as using a similar method to your (since it is the inverse-link applied to the linear predictor).
I have followed steps 1-5 in your tutorial ok but I am slightly confused after this point.
I have two points below:
1.0 For my requirement, I understand that I will need to build the linear predictor sample by adding together both the fixed effects and random effects in my model, which requires accounting for each latent element (and their start index).
Since I ran link=1 (which obtained values for the defined NA observations I wanted to predict, when I called INLA), technically I could just just run the below code for predicting on an observation already known? (an NA observation in the dataset, with a value predicted by INLA) Correct me if i'm wrong. From: 7 Predicting for a different experiment on an old plate
our.experiment = list(x1 = 0, x2 = 0, plate = 7)
nr = 57
s = samples[[nr]]$latent
## beta1 * 0
f.x1.0 = s[44, , drop=F] * our.experiment$x1
## beta1 * 1
f.x2.1 = s[45, , drop=F] * our.experiment$x2
## f(plate = 7)
# - the same plate as before
f.plate.7 = s[28, , drop=F]
## The intercept
int = s[43, , drop=F]
predictor.our.experiment = drop(f.x1.0 + f.x2.1 + f.plate.7 + int) # see below for questionsamples.link = inla.link.invlogit(predictor.our.experiment
2.0
To get data values, by transforming the predictor sample, you use:
"We now need to transform through the link function and the data sampling. We know we use the logit transform..."
samples.link7 = inla.link.invlogit(samples.pred7)
I haven't called in INLA() control.family() and can't call it from:
res1$.args$control.family$control.link$model Firstly, this has to be called to obtain data sampling? Secondly, since my model uses Poisson likelihood, I guess I would need to use log?
control.family1 = list(control.link=list(model="logit") == control.family1 = list(control.link=list(model="log")I would be grateful for any assistance.
Tom.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
--
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-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
--
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-group+unsub...@googlegroups.com.
Haakon Bakka
Jul 26, 2018, 2:49:39 AM7/26/18
to Thomas, R-inla discussion group
If you have a log-link, then this transformation, to get from the linear predictor to the data scale, is the exponential ("exp").
Haakon
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion-group+unsubscr...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
--
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-group+unsubscr...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
--
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-group+unsubscr...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
--
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-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-discussion-group@googlegroups.com.
Thomas
Jul 27, 2018, 6:59:10 AM7/27/18
to R-inla discussion group
Hi Haakon,
It seems I was confused between the family and family.control (I guess from working late at night!), what you've wrote makes sense.
Regards,
Tom.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
--
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-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
--
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-group+unsub...@googlegroups.com.
To post to this group, send email to r-inla-disc...@googlegroups.com.
Visit this group at https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.
--
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-group+unsub...@googlegroups.com.
Reply all
Reply to author
Forward
0 new messages