I have the following model for which i am trying to plot predictions so to analyse the spatial-temporal trends over the years of my data (2012-2014).
mod = gam(TB ~ offset(log(Population)) s(Indigenous, k = 10, bs = "cr") s(Urbanisation, k = 10, bs = "cr")
s(Density, k = 10, bs = "cr") Poverty s(Poor_Sanitation, k = 10, bs = "cr") Unemployment Timeliness
as.factor(Year) Region s(lon, lat),
data = TBdata,
family = nb(link = 'log'))
Currently I have used the following code to plot my model but this doesn't give me temporal trends (just one plot not split by year), i also think this is just plotting my smooth function of lon and lat. If anyone knows how to correctly plot the predictions of this model so i can get the spatial-temporal trends that would be great (ggplot preferred by any method is fine).
my current plot
plot(mod, contour.col = 'white', too.far = 0.10, scheme = 2, rug = F)
output
data snipet
TBdata %>% as.data.frame()
Indigenous Illiteracy Urbanisation Density Poverty Poor_Sanitation Unemployment Timeliness Year TB Population Region lon lat
1 0.33498668 6.354471 84.10664 0.7135083 31.27839 15.269975 5.411208 59.23954 2012 323 559543 11001 -60.67326 -12.1375562
2 6.44677661 8.493550 71.36292 0.7429604 48.59763 29.364664 5.916373 58.06452 2012 15 73193 11002 -64.02859 -9.4315409
3 0.16009349 10.211680 66.96757 0.6638205 42.28298 31.923537 4.084887 88.65248 2012 47 176003 11003 -62.35638 -10.6366765
4 0.54930178 10.294286 68.67525 0.5951172 38.89858 28.534552 5.497259 75.15528 2012 45 294157 11004 -64.09781 -11.6949147
5 0.31773624 11.418378 42.84595 0.6385168 51.57131 48.406288 4.035860 58.33333 2012 17 69841 11005 -61.21610 -13.1534602
6 1.02317144 10.267994 66.49738 0.5946587 39.10929 29.287196 5.578758 69.31818 2012 66 228654 11006 -61.68878 -11.9327327
7 0.82154405 7.811040 82.75737 0.5927965 31.71298 12.124683 6.040940 88.76404 2012 24 136730 11007 -62.56662 -9.6590901
8 0.11125568 11.782911 65.34103 0.5550268 40.84772 29.258613 5.059368 41.17647 2012 5 51890 11008 -62.88370 -11.6058899
9 3.26071252 22.840871 57.55142 1.0104743 60.54176 32.830252 6.722286 50.00000 2012 59 135803 12001 -71.12488 -8.6984193
10 8.27981067 32.395581 51.29996 1.2319246 65.52161 44.536382 4.761273 45.00000 2012 13 76221 12002 -69.65906 -9.6142587
11 7.81641388 26.125845 63.67133 1.0323498 58.63709 35.371653 7.526261 58.33333 2012 13 52651 12003 -67.92047 -9.9827217
12 0.18487825 11.014522 83.02779 0.7966288 42.60303 11.748025 8.647023 31.61616 2012 246 433224 12004 -72.86194 -8.2736825
13 1.23080070 19.107763 66.32175 0.8262261 55.68908 28.966361 5.895377 70.90909 2012 19 60887 12005 -69.15972 -10.7140791
14 50.64622783 23.387271 48.26342 1.4629935 48.07222 45.458143 7.683733 36.09023 2012 48 99826 13001 -65.46983 -4.0204889
15 10.04667257 21.056304 48.48717 1.5299792 72.38631 47.042878 15.439268 0.00000 2012 0 25044 13002 -59.92364 -1.5879241
16 27.62278330 21.747460 58.06581 1.4155145 70.70161 39.524745 6.358736 43.02326 2012 103 231482 13003 -65.13166 -0.3415058
17 4.30755994 33.107526 62.22066 1.1412305 72.10023 40.095659 6.539286 48.36066 2012 42 131695 13004 -65.10241 -7.1406849
18 3.01262696 15.785024 74.05870 1.2235418 66.95602 20.833828 5.363723 30.39216 2012 43 87520 13005 -58.19999 -2.9087320
19 0.94701077 18.122909 60.82228 1.1627655 66.90984 35.127488 9.134779 39.74359 2012 46 156671 13006 -60.24980 -3.3062838
20 0.65529611 4.908563 93.81673 0.8593049 41.11021 4.893957 10.483849 75.53917 2012 1829 2106861 13007 -60.54443 -6.6830480
21 1.35932242 9.200504 47.65380 0.8861087 55.57790 37.740240 10.428850 50.00000 2012 32 55600 13008 -68.71086 -5.9785253
22 0.64001789 9.758551 60.98325 1.1415410 67.23815 30.632589 8.406879 59.20398 2012 51 156384 13009 -67.55678 -1.8679793
23 5.00288446 8.046211 57.06692 1.2392834 73.76746 35.573960 7.824595 45.21452 2012 77 248506 13010 -58.56150 -3.1484377
24 5.98590106 27.766366 58.63355 1.1346604 71.10207 39.876155 8.986358 46.00000 2012 22 49500 13011 -63.22067 -4.1223522
25 7.64428333 29.406475 59.70856 1.0691215 70.59885 42.578621 6.897342 54.79452 2012 12 70911 13012 -68.14551 -8.1089091
26 6.27901221 16.298789 54.29879 1.1055992 65.44075 42.302509 8.442431 40.58824 2012 56 170985 13013 -69.69795 -4.7593017
27 8.39040188 8.206571 89.95105 0.7979057 38.37151 6.075713 7.358698 64.10891 2012 91 334076 14001 -60.49593 0.3182771
28 44.71062619 17.073994 22.33634 1.3535825 57.66587 48.465597 9.245283 53.84615 2012 11 44023 14002 -60.24954 3.4039662
29 6.18822518 18.247396 57.11455 0.9473904 53.79555 32.518013 8.424442 40.00000 2012 7 43635 14003 -61.74950 1.4791942
30 2.83585991 13.326820 50.54670 0.8424786 56.25717 31.095345 6.860730 68.75000 2012 12 47790 14004 -62.22345 3.4860616
31 2.09288322 10.842066 53.12386 1.1715843 70.29388 31.457975 7.813915 46.47887 2012 40 189429 15001 -49.92266 -4.5450924
32 0.63918534 10.374821 60.78618 0.9561001 62.77747 27.277026 7.956556 43.20988 2012 131 492734 15002 -48.68255 -2.5733129
33 1.21355631 17.414926 51.12920 1.2375818 65.28844 36.654216 8.156438 56.09756 2012 14 69092 15003 -51.95188 -7.4381955
34 0.01616554 28.968886 38.94816 1.6750609 77.88347 57.851494 9.060787 54.28571 2012 22 135041 15004 -54.65470 -1.8714902
35 0.03968371 24.890572 42.31165 1.4990517 75.79085 55.581917 7.924278 63.28125 2012 38 209923 15005 -47.77654 -0.8198963
36 0.04052552 13.905926 48.69076 1.1558360 71.62642 38.544809 8.836839 43.43434 2012 36 158333 15006 -49.59232 -7.1006744
37 0.16146379 3.754233 94.76295 0.7778492 40.30134 2.661756 10.683572 62.49457 2012 1880 2185084 15007 -51.01589 -2.1353053
38 0.08202806 9.219135 75.28922 0.8238219 57.50873 14.657248 8.974412 52.50836 2012 185 305405 15008 -50.15950 -6.3840531
39 0.05947901 11.416480 53.22522 0.9225446 67.27691 15.432007 6.932257 37.11656 2012 110 251888 15009 -48.05151 -3.8479350
40 0.08125090 17.971794 56.30766 0.8787671 69.37043 25.782833 7.573189 36.34204 2012 131 390863 15010 -56.87939 -0.1044528
41 0.03864606 13.744725 49.50903 1.2463430 74.00576 42.667051 7.309062 75.55556 2012 101 449907 15011 -49.66102 -5.6414827
42 0.23364324 17.398443 49.76983 1.1178318 68.50080 36.631718 9.384440 55.00000 2012 88 299079 15012 -56.42347 -6.4524857
43 0.23725366 22.500687 38.31865 0.9708972 74.23657 39.255409 5.860432 49.11504 2012 81 441719 15013 -47.18900 -2.0744028
44 4.24519522 15.341410 56.05853 0.8841889 56.77666 32.751180 6.477704 34.27419 2012 72 240165 15014 -50.37440 -0.8876770
45 1.87527384 16.205226 58.56287 0.8116198 54.69351 38.918163 5.465688 56.88073 2012 67 283554 15015 -50.19400 -8.7018861
46 0.52965071 16.844620 67.98471 0.8407914 62.00125 27.628233 9.747115 60.46512 2012 68 340793 15016 -48.03117 -1.3843356
47 0.67886709 16.821733 71.07965 0.8295509 59.87976 19.369363 10.857736 67.46411 2012 71 306994 15017 -49.53206 -2.3657026
48 2.54441144 12.853355 57.68736 0.8028324 50.58523 36.205599 7.418385 57.42574 2012 50 177935 15018 -47.07320 -1.1491093
49 0.51103595 10.931050 75.61964 0.7753638 45.60888 17.048489 10.320373 54.72441 2012 66 271172 15019 -48.49396 -1.3613721
50 0.18999389 13.147962 74.79087 0.7868254 50.54423 19.089217 9.036023 66.66667 2012 85 295238 15020 -49.36448 -0.5230343
51 0.50384846 14.530229 75.98996 0.6910516 50.60517 22.919049 7.678848 57.22222 2012 54 185608 15021 -53.35841 -5.6680040
CodePudding user response：
Producing the smoothed map from a gam is a complex affair, so it's probably best to use the data wrangling that takes place inside plot.gam - helpfully, this actually return data that we can use to plot in whichever way we like.
First we load your workspace and appropriate packages
load("../../Downloads/TBdata.RData")
library(ggplot2)
library(mgcv)
Next we split your data into year groups and carry out your gam on each. We then store the data output of the plot for each of the three years and convert it into a data frame. The resulting data frames are bound together into a single large data frame:
df <- do.call(rbind, lapply(split(TBdata, TBdata$Year), function(x) {
p <- plot(gam(TB ~ offset(log(Population))
s(Indigenous, k = 10, bs = "cr")
s(Urbanisation, k = 10, bs = "cr")
s(Density, k = 10, bs = "cr")
Poverty
s(Poor_Sanitation, k = 10, bs = "cr")
Unemployment
Timeliness
Region
s(lon, lat),
data = x,
family = nb(link = 'log')),
contour.col = 'white', too.far = 0.10, scheme = 2, rug = F)[[5]]
df <- expand.grid(lon = p$x, lat = p$y)
df$TB <- p$fit[,1]
df$Year <- x$Year[1]
df
}))
Having the data in this format makes the plot itself easy using ggplot:
ggplot(df, aes(lon, lat, z = TB, fill = TB))
geom_tile()
geom_contour(color = "white", na.rm = TRUE)
facet_grid(.~Year)
coord_equal()
scale_fill_gradientn(colors = c("red", "yellow", "white"),
values = c(0, 0.9, 1),
na.value = "white")
theme_minimal(base_size = 16)
