Introduction

In this project I wanted to show some methods for handling spatial data given for a set of points. I wanted to show that it is possible to estimate the values from point to the whole country (here: Poland). I hope that my contribution will help others to learn from the data given by points about the patterns for bigger teritory.

Air pollution - why this topic?

There is more and more discussion around the world about air pollution and pollution-related ailments. One of the most dangerous particulate matter which causes lung cancer, among other illnesses, is PM2.5.

One of the challenges with analysing air pollution and PM2.5 is the data availability. Air measurement stations are usually scattered and capture air quality only at specific locations. With my project I want to suggest an approach to use point-based data also for calculating and visualising air pollution estimates for the bigger area like the whole country.

What methods did I use?

Statistical tools I employ include Gaussian Process calibration and estimation, bandwith estimation for point patterns, spatial interpolation, Generalized Additive Model, Gaussian Process and OLS.

Speaking of packages, I use a dozen of them for spatial analysis and visualisation - sf, sp, maptools, raster, osm, ggplot2, ggspatial, spatstat, and others. For dealing with time-series I use lubridate; tidyr and dplyr for less messy code; and stringr for text manipulation.

Where does data come from?

I utilize publicly available data provided by the Chief Inspectorate For Environmental Protection in Poland (https://powietrze.gios.gov.pl/pjp/archives). Shapefiles with administrative boundaries within Poland come from GIS Support (https://gis-support.pl/baza-wiedzy-2/dane-do-pobrania/granice-administracyjne/).

Packages and introductory data

# Check your working directory
getwd()
## [1] "C:/Users/tymme/Studia_nauka_niechronione/Seminars_R_ML/Air_pollution_Poland"
# Loading necessary libraries and packages
library(sf) # spatial data frames
library(ggplot2) # plotting
library(viridis) # colourful maps
library(stringr) #strings operations
library(dplyr) # handling the data easier
library(tidyr) # less messy code + lengthening/widening data
library(openxlsx) # Opening excel files, changing data types
library(sp) # for spatial polygons maps
library(spdep) # for spatial statistics
library(spatialreg) # for spatial regressions
library(raster) # For nice maps  with osm background
library(ggspatial) # For nice maps  with osm background
library(lubridate) # For operations on dates
library(mgcv) # for modelling (e.g. Marginal Likelihood, Generalized Cross Validation)
library(epm) # for some spatial calculations
library(glue) # for "f-strings" in R
library(broom) # transforms model outputs into tibble
library(lava) # for %ni% operator
library(factoextra) # for clustering
library(spatstat) # for density estimation
library(maptools) # point pattern processing
library(ggthemes) # for nicer themes in ggplot

Reading shapefiles

With spatial analysis it is crucial to load appropiate shapefiles. They are available free of charge for countries, continents, administrative regions, rivers and much more. Here I read boundries of the whole country (Poland), “województwa” - regions and “powiaty” - counties. After reading I need to unify projections.

Polska <- read_sf("Polska/polska.shp")
WOJ <- st_read("wojewodztwa/wojewodztwa.shp", stringsAsFactors = FALSE)
## Reading layer `wojewodztwa' from data source 
##   `C:\Users\tymme\Studia_nauka_niechronione\Seminars_R_ML\Air_pollution_Poland\wojewodztwa\wojewodztwa.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 16 features and 29 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 171677.6 ymin: 133223.7 xmax: 861895.7 ymax: 774923.7
## Projected CRS: ETRS89 / Poland CS92
WOJ$jpt_nazwa_<-iconv(WOJ$jpt_nazwa_, "WINDOWS-1250", "UTF-8") # introducing Polish characters
WOJ <- WOJ[,c(5,6,30)] # leaving only necessary columns

POW <- st_read("powiaty/powiaty.shp", stringsAsFactors = F)
## Reading layer `powiaty' from data source 
##   `C:\Users\tymme\Studia_nauka_niechronione\Seminars_R_ML\Air_pollution_Poland\powiaty\powiaty.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 380 features and 29 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 171677.6 ymin: 133223.7 xmax: 861895.7 ymax: 774923.7
## Projected CRS: ETRS89 / Poland CS92
POW$jpt_nazwa_ <- iconv(POW$jpt_nazwa_, "WINDOWS-1250", "UTF-8") # introducing Polish characters
POW <- POW[,c(5,6,30)] # leaving only necessary columns

# Unifying projections
WOJ <- st_transform(WOJ, 4326)
POW <- st_transform(POW, 4326)
Polska <- st_transform(Polska, 4326)

Loading measurement stations

Now I read measuring stations (whichever particles they measure). I want to gain some insights from location of those stations into possible pollution distribution. Hypothetically, where worse air is, the people are more concerned about the problem and therefore more stations should be in operation.

stacje <- read.xlsx("Data/Metadane oraz kody stacji i stanowisk pomiarowych.xlsx", sheet = 1)

# Transform excel time format to universal date format
stacje$Data.uruchomienia <- convertToDateTime(stacje$Data.uruchomienia, origin = "1900-01-01")
stacje$Data.zamknięcia <- convertToDateTime(stacje$Data.zamknięcia, origin = "1900-01-01")

# Leaving only stations still in operation
stacje_czynne <- stacje %>%
  filter(is.na(Data.zamknięcia)==T)
#Creating df, that I will use later (after loading pollution data)
pm25 <- stacje_czynne
# Transforming df to sf format
stacje_czynne_sf <- st_as_sf(stacje_czynne,coords = c(colnames(stacje_czynne)[15], colnames(stacje_czynne)[14]) ,crs=4326)

Reading pollution data

Now I read PM2.5 data from all measuring stations in Poland. Data I use is calculated for every 24 hours. For my analysis I aggregate those values for each month, calculating means. Among many upsides, it allows me to overcome the problem of missing records - e.g. when some stations used to break down for a couple of days.

pollution_pm25 <- read.xlsx("Data/Wyniki pomiarów z 2022 roku/2022_PM25_24g.xlsx")

# Renaming columns and deleting redundant rows
colnames(pollution_pm25) <- pollution_pm25[1,]
colnames(pollution_pm25)[1] <- "Time"
pollution_pm25 <- pollution_pm25[c(6:nrow(pollution_pm25)),]
rownames(pollution_pm25) <- 1:nrow(pollution_pm25)

# Change data format from character to numeric
for (i in 1:ncol(pollution_pm25)){
  pollution_pm25[,i] <- as.numeric(pollution_pm25[,i])
}

# Change excel time format to normal time format
pollution_pm25[,1] <- convertToDateTime(pollution_pm25[,1], origin = "1900-01-01")

# Extracting from which month the observation is (package lubridate)
pollution_pm25$month <- month(ymd(pollution_pm25$Time))

# Counting monthly averages of pollution
pollution_pm25_monthly <- pollution_pm25[c(1:12),]

for ( i in 1:12){
  monthly <- pollution_pm25 %>% filter(month==i) %>%
    summarise(across(all_of(colnames(pollution_pm25)),
                     ~mean(.x,na.rm = T)))
  for (j in 2:ncol(pollution_pm25_monthly)){
    pollution_pm25_monthly[i,j] <- monthly[nrow(monthly),j]
  }
}

# Creating df with pollution values
# Adding columns for months
pm25$January <- 0
pm25$February <- 0
pm25$March <- 0
pm25$April <- 0
pm25$May <- 0
pm25$June <- 0
pm25$July <- 0
pm25$August <- 0
pm25$September <- 0
pm25$October <- 0
pm25$November <- 0
pm25$December <- 0

# Subsetting all stations in operation
pm25 <- pm25[pm25[,2] %in% colnames(pollution_pm25_monthly),]

# Copying values to pm25 dataframe
for (i in 1:nrow(pm25)){
  name_pm25 <- pm25[i,2]
  for (j in 2:ncol(pollution_pm25_monthly)){
    name_pollution <- colnames(pollution_pm25_monthly)[j]
    if (name_pm25==name_pollution){
      for ( k in 1:12){
        pm25[i,k+15] <- pollution_pm25_monthly[k,j]
      }
    }
  }
}

# Creating sf with pollution values
pm25_sf <- st_as_sf(pm25, coords = c(colnames(pm25)[15], colnames(pm25)[14]) ,crs=4326)

Measurement stations plotting

Here I plot measurement stations as a simple point pattern. To every plot I add graticules and borders of the regions.

All stations

At first I present all stations that measure any air features and were active in 2022.

ggplot() +
  geom_sf(data=WOJ) +
  geom_point(aes(x=WGS84.λ.E, y=WGS84.φ.N),
             data=stacje_czynne, color="darkred", size=1.5)+
  scale_fill_viridis(option = "H")+
  ggtitle(label = "Measurement stations in Poland")+
  xlab(label="")+ylab(label = "")+
  theme_gdocs()

Only stations active for PM2.5 measures

Here I plot locations of only these stations, which were measuring PM2.5. I exclude stations measuring other air features like PM10, benzene or nitrogen dioxide.

ggplot()+
  geom_sf(data=WOJ) +
  geom_point(aes(x=WGS84.λ.E, y=WGS84.φ.N,
                 col=January),
             data=pm25, )+
  scale_color_viridis(discrete=F)+
  ggtitle(label = "Stations actively measuring PM2.5")+
  guides(color=guide_legend(title = "Pollution in January"))+
  xlab(label="")+ylab(label = "")+
  theme_gdocs()

Map of measuring stations with osm background

Again, I show here only stations measuring PM2.5.

ggplot(pm25_sf) + 
  annotation_map_tile(type="osm",zoom=7) + 
  geom_sf(col='red',size=1.5)
## Zoom: 7

Plotting measurement stations - smoothed

In this section, I use smoothing to calculate so-called “stations intensity”. The idea is that the area with higher “stations intensity” is in the proximity of numerous stations and close to them. The area with low “stations intensity” is far away to even a single station.

Here I use various classes not used before: ppp and raster. At first I change projection to obtain cartographic projection. (My idea with bandwith estimation does not work with geographic estimation.) Then, I calculate bandwidths using four different estimations methods. The outcomes of each of them will be displayed below.

At first I adjust projections.

Polska <- st_transform(Polska, crs=2180)
stacje_czynne_sf <- st_transform(stacje_czynne_sf, crs=2180)
pm25_sf <- st_transform(pm25_sf, crs=2180)

Then I create the window necessary for ppp class.

winn_stations <- as.owin(as_Spatial(Polska))

Below I save all working stations in ppp class.

xy_stations <- st_coordinates(stacje_czynne_sf)  # retrieve coordinates
PP_stations <- as.ppp(xy_stations, winn_stations)
PP_stations <- as.ppp(PP_stations)

Then I save to ppp class only working stations that actively measure PM2.5.

xy_stations_2 <- st_coordinates(pm25_sf)  # retrieve coordinates
PP_stations_2 <- as.ppp(xy_stations_2, winn_stations)
PP_stations_2 <- as.ppp(PP_stations_2)

Now I estimate the bandwith/degree of smoothing. At first I use all stations, then those measuring PM2.5. (Warning! calculating takes some time.)

scott_stations <- bw.scott(PP_stations) # Scott’s rule
ppl_stations <- bw.ppl(PP_stations) # Likelihood cross validation
diggle_stations <- bw.diggle(PP_stations) # Diggle’s approach
CvL_stations <- bw.CvL(PP_stations) # Cronie and van Lieshout’s criterio

scott_stations_2 <- bw.scott(PP_stations_2) # Scott’s rule
ppl_stations_2 <- bw.ppl(PP_stations_2) # Likelihood cross validation
diggle_stations_2 <- bw.diggle(PP_stations_2) # Diggle’s approach
CvL_stations_2 <- bw.CvL(PP_stations_2) # Cronie and van Lieshout’s criterio

Calculating density - all stations

I use Jones Diggle Estimator for edge correction. It is a bit slower, but more accurate than the default density option in “density” command.

Scott’s rule

d_scott_dig_stations <- density(PP_stations,
                                sigma = scott_stations,
                                positive=TRUE,
                                leaveoneout=FALSE,
                                diggle=TRUE)
plot(d_scott_dig_stations,
     main="Intensity of all active stations")

Likelihood Cross Validation

d_ppl_dig_stations <- density(PP_stations,
                          sigma = ppl_stations,
                          leaveoneout=FALSE,
                          positive=TRUE,
                          diggle=TRUE)
plot(d_ppl_dig_stations,
     main="Intensity of all active stations")

Diggle’s approach to smoothing bandwidth

d_diggle_dig_stations <- density(PP_stations,
                          sigma = diggle_stations,
                          leaveoneout=FALSE,
                          positive=TRUE,
                          diggle=TRUE)
plot(d_diggle_dig_stations,
     main="Intensity of all active stations")

Cronie and van Lieshout’s criterion

d_cvl_dig_stations <- density(PP_stations,
                       sigma = CvL_stations,
                       leaveoneout=FALSE,
                       positive=TRUE,
                       diggle=TRUE)
plot(d_cvl_dig_stations,
     main="Intensity of all active stations")

Calculating density - stations in operation

Again, I use Jones Diggle Estimator for edge correction. It is a bit slower, but more accurate than default density option.

Scott’s rule

d_scott_dig_stations_2 <- density(PP_stations_2,
                                sigma = scott_stations_2,
                                positive=TRUE,
                                leaveoneout=FALSE,
                                diggle=TRUE)
plot(d_scott_dig_stations_2,
     main="Intensity of stations measuring Pm2.5")

Likelihood Cross Validation

d_ppl_dig_stations_2 <- density(PP_stations_2,
                          sigma = ppl_stations_2,
                          leaveoneout=FALSE,
                          positive=TRUE,
                          diggle=TRUE)
plot(d_ppl_dig_stations_2,
     main="Intensity of stations measuring Pm2.5")

Diggle’s approach to smoothing bandwidth

d_diggle_dig_stations_2 <- density(PP_stations_2,
                          sigma = diggle_stations_2,
                          leaveoneout=FALSE,
                          positive=TRUE,
                          diggle=TRUE)
plot(d_diggle_dig_stations_2,
     main="Intensity of stations measuring Pm2.5")

Cronie and van Lieshout’s criterion

d_cvl_dig_stations_2 <- density(PP_stations_2,
                       sigma = CvL_stations_2,
                       leaveoneout=FALSE,
                       positive=TRUE,
                       diggle=TRUE)
plot(d_cvl_dig_stations_2,
     main="Intensity of stations measuring Pm2.5")

Further exploring stations intensity

The plot above shows the density of stations. The intuition suggests that the more dense area may be 1) more populated or 2) more polluted. To check what is the underlying cause, I’ll check both density and pollution. I use just one of the bandwidths (Likelihood Cross Validation for all stations), but one can try any of them. I map estimations on powiats for better visualisation + for further modelling. Firstly, I convert ppp class to raster and plot the map.

stations_intensity <- raster(d_ppl_dig_stations)
plot(stations_intensity,
     xaxt = "n",
     yaxt = "n")

Then, I unify projections of “POW” and “stations_intensity”. The latter should not have any projection yet. Let’s check it

projection(stations_intensity) # check projection
## [1] NA

Indeed, the projection is “NA”. Therefore, I need to set it. The projection should be a cartographic one, not geographic. So, I use so-called ‘PUWG 1992’ or in longer version - ’ Państwowy Układ Współrzędnych Geodezyjnych 1992’. Its EPSG code is 2180.

After unifying projections I will convert “stations_intensity” to sf class, using POW object as a benchmark. As the effect, I will obtain a mean value of the intensity in every powiat.

projection(stations_intensity) <- st_crs(2180)$proj4string
POW <- st_transform(POW, crs=st_crs(2180)$proj4string)
# Convert raster to sf class
stations_intensity2 <- rasterToGrid(stations_intensity, target = POW, crop = T, fun = 'mean')

Now I will plot the results

plot(stations_intensity2,
     main="Stations intensity by powiat")

ggplot()+
  geom_sf(data=stations_intensity2, mapping=aes(geometry=geometry, fill=vals))+
  scale_fill_viridis(option = 'C')+
  theme_gdocs()+
  ggtitle(label = "Stations intensity by powiat")

At last I will add intensity values to the POW object in order to easier handle all the data. Importantly, both objects have the same order of geometries, so the assignment requires only one line of code. To check the correctness I will make a plot using POW object. It should be identical to the one above.

POW$stations <- stations_intensity2$vals 
ggplot()+
  geom_sf(data=POW, mapping=aes(geometry=geometry, fill=stations))+
  scale_fill_viridis(option = 'C')+
  theme_gdocs()+
  ggtitle(label = "Stations intensity by powiat")

Checking population

Here I visually check whether population plot and stations intensity plot are similar. At first I load the relevant data.

population <- read.csv("Data/population_density_2022.csv", sep = ";", dec = ",")

Then, I ensure that TERYT codes, that I use to merge dataframes would be camparable across datasets (written in the same format, without redundant values like 0 at the beginning and so on). At first I transform character code to numeric format in both dataframes

# 1. Ensuring the code in of the same length in every cell
for (i in 1:nrow(population)){
  if (str_count(population[i,1], pattern = "[:digit:]")==6){
    population[i,1]=paste('0',population[i,1], sep='')
  }
}

# 2. Removing "000" from the end
population$Kod=sub(population$Kod, pattern="000", replacement="")

# 3. Saving a value as a number
population$Kod <- as.numeric(population$Kod)
POW$jpt_kod_je <- as.numeric(POW$jpt_kod_je)

# 4. Removing unncessary rows from population dataframe and renaming them
population <- population[,c(1,3)]
colnames(population) <- c("jpt_kod_je","pop_density")

# 5. Merging dataframes
POW <- left_join(POW, population, by = join_by(jpt_kod_je))

Finally, I can map population density for every powiat in Poland.

ggplot()+
  geom_sf(data=POW, mapping=aes(geometry=geometry, fill=stations))+
  scale_fill_viridis(option = 'C')+
  theme_gdocs()+
  ggtitle(label = "Stations intensity by powiat")

ggplot()+
  geom_sf(data=POW, mapping=aes(geometry=geometry, fill=log(pop_density)))+
  scale_fill_viridis(option = 'C')+
  theme_gdocs()+
  ggtitle(label = "Population density (log)")

It looks like the population density is valid for Silesia and Warsaw. However, it is not enough to explain the whole pattern of stations. I need to check pollution of PM 2.5

Model fitting - one selected month

The point of this section is to examine where in Poland pollution of PM2.5 is the highest. To do so, I fit a model using Gausian Process and point-measured pollution values from a selected month (here I use January). As a result, I will obtain a proxy for real pollution in given geographic entity.

Firstly, I transform the sf object containing the locations and pollution values from active stations measuring PM2.5. I need to save it to a tibble, because it is a format required by gam function, which I will utilize to estimate Gausian Process.

pm25_sf <- st_transform(pm25_sf, crs=4326)
pm25_tibble <- pm25_sf %>% st_coordinates() %>% 
  as_tibble() %>%
  bind_cols(pm25_sf %>% 
              st_drop_geometry())

Choosing k paramenter for modelling - via the AIC

Then, I should decide which model will be the most suitable for estimating. The parameter which I can manipulate is k, that controlls the level of smoothing. I’ll choose k, for which AIC of the model is the lowest.

res <- tibble(k=5:55) %>% 
  rowwise() %>% 
  mutate(
    AIC_k =
      AIC(
        gam(
          as.formula(glue(
            "January~s(X,Y,k={k},m=2,bs='gp')")),
          data=pm25_tibble)))
# And now plot the results
ggplot(res,aes(x=k,y=AIC_k)) + 
  geom_point(col='darkred') +
  geom_line()+
  theme_gdocs()

From the plot above I see, that the lowest AIC is somewhere between 35 and 40. So, I repeat the computation and plot again.

res2 <- tibble(k=35:40) %>% 
  rowwise() %>% 
  mutate(
    AIC_k =
      AIC(
        gam(
          as.formula(glue(
            "January~s(X,Y,k={k},m=2,bs='gp')")),
          data=pm25_tibble)))
# And the plot
ggplot(res2,aes(x=k,y=AIC_k)) + 
  geom_point(col='darkred') +
  geom_line()+
  theme_gdocs()

The lowest AIC is for k=38

Fitting the Gaussian Process model

Once I calculated the optimal k, I can fit a model. I define a smoothing function by setting marginal penalty to 2 (m=2), but one can chenge it if wishes to.

gp <- gam(January~s(X,Y,k=38,m=2,bs='gp'),data=pm25_tibble) 

gp
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## January ~ s(X, Y, k = 38, m = 2, bs = "gp")
## 
## Estimated degrees of freedom:
## 11  total = 12.02 
## 
## GCV score: 16.92855

Visualise The Gaussian Process using powiats - only January

Having fitted the Gaussian Process model for pollution in January, I want to visualise it properly. I will switch the projection back to WGS84 and then assign the values predicted with gaussian process to the POW object. The value for every powiat will be equal to the value at the centroid for this powiat.

POW <- st_transform(POW, crs = 4326)
powiats <- POW %>% st_centroid() %>% 
  st_coordinates() %>% as_tibble() %>% 
  predict(gp,.)
## Warning: st_centroid assumes attributes are constant over geometries
POW$predicted_gp <- powiats

Now it’s time for plotting the results.

ggplot(POW,aes(fill=predicted_gp,col=predicted_gp)) +
  geom_sf() +
  scale_fill_viridis_c()+
  guides(fill=guide_legend(title = "Pollution (μg)"),
         color="none")+
  ggtitle(label = "Pollution of PM 2,5 in air",
          subtitle = "January 2022")+
  theme_gdocs()

The pot clearly indicates that the best quality of air (in terms of PM2.5 pollution) is by the sea, near the border with Germany in the west and at the north-east corner of Poland. However, one should keep in mind, that the calculation here is based only on the values for one selected month - January. Since pollution patterns may change through the year, it is advised to treat the outcome of this model with caution. Additionally, the result is calculated based on the limited number of stations. It would be wise to check standard errors of my estimates then.

Estimating Standard Errors of GP Values

In this section I want to examine where my estimates are the most plausible and where the chances of me being wrong are higher. Below I calculate standard errors for a centroid of every powiat and plot them. I add the location of every station for the reference.

powiats_se <- POW %>% st_centroid() %>% 
  st_coordinates() %>% as_tibble()
## Warning: st_centroid assumes attributes are constant over geometries
powiats_se <- powiats_se %>%
  predict(gp,.,se.fit = TRUE) %>% .$se.fit %>%
  mutate(powiats_se,se=.)
POW$predicted_gp_se <- powiats_se$se
ggplot(POW) +
  geom_sf(aes(fill=predicted_gp_se,col=predicted_gp_se)) +
  scale_color_viridis_c(trans='sqrt') +
  scale_fill_viridis_c(trans='sqrt') +
  geom_sf(data=pm25_sf,col='red',mapping=NULL)+
  theme_gdocs()

Embedding GPs into all-year model

Here I estimate the model similar to what I did above. Here, however, I use the data from the whole year, adding monthly effects to account for the fixed seasonal effect (I assume it is constant across different parts of Poland).

I need to create a sf object with all months at once. Hence, I need to lengthen my data.

pm25_sf_long <- pivot_longer(pm25_sf, 
                             cols = January:December, 
                             names_to = "Month",
                             values_to = "Pollution")

Similarly as before, I transform the sf object to a tibble. This time I use the data for all months.

pm25_tibble_long <- pm25_sf_long %>% st_coordinates() %>% 
  as_tibble() %>%
  bind_cols(pm25_sf_long %>% 
              st_drop_geometry())

And I fit the second Gaussian Process model. Parameters k and m remains unchanged.

gp2 <- gam(Pollution~s(X,Y,k=38,m=2,bs='gp')+Month+0,
           data=pm25_tibble_long)

Plotting model GP2 for January

In the second model I obtained the values for all months. However, for easier comparison with the previous plot, I will start with showing pollution estimates for January.

powiats_January <- POW %>% st_centroid() %>% 
  st_coordinates() %>% as_tibble() %>% 
  mutate(Month='January')%>%
  predict(gp2,.)
## Warning: st_centroid assumes attributes are constant over geometries
POW$predicted_January <- powiats_January

And now plotting…

ggplot(POW,aes(fill=predicted_January,col=predicted_January)) +
  geom_sf() +
  scale_fill_viridis(option = "H")+
  guides(fill=guide_legend(title = "Pollution (μg)"),
         color="none")+
  ggtitle(label = "Pollution of PM 2,5 in air",
          subtitle = "January 2022")+
  theme_gdocs()

Monthly effects in GP2 model

The model for all months calculated the fixed monthly effects. I extract them to check duting what time of a year pollution of PM2.5 is on average the worst.

months <- c('January','February','March','April','May','June',
            'July','August','September','October','November','December')
effects <- tidy(gp2, parametric=TRUE) %>% 
  mutate(Month = str_remove(term,"Month") %>%
           factor(levels=months,ordered=TRUE)) %>% 
  arrange(Month) 
ggplot(effects,aes(x=Month,y=estimate)) + 
  geom_col()+
  theme_gdocs()

There is little pollution April-September. The peak from this model is in March.

Embedding GPs into other model 2 -> 12 processes at once

Here I run 12 separate Gaussian Processes for every month each

pm25_tibble_long <- pm25_tibble_long %>%  # gam wants factors not character variables
  mutate(Month=factor(Month,levels=months))

gp3 <- gam(Pollution~s(X,Y,k=38,m=2,bs='gp',by=Month)+0,data=pm25_tibble_long)

Having fitted the third model I want to compare its AIC with the previous one.

AIC(gp2,gp3)
##           df      AIC
## gp2 44.39189 3585.067
## gp3 25.00000 5844.293

From results it looks as though the simpler (second) model is more plausible then the third one. I can look at gp3 anyway, to obtain new insights from the data.

powiats3 <- NULL
for (m in months) {
  powiats3 <- powiats3 %>% 
    bind_rows(POW %>% 
                mutate(Month=m))
}
powiats3 <- powiats3 %>% 
  mutate(Month=factor(Month,levels=months)) 

powiats3 <- powiats3 %>% # it takes some time
  st_centroid() %>%
  st_coordinates() %>% 
  as_tibble() %>% bind_cols(powiats3,.) %>%
  mutate(pred=predict(gp3,newdata=.))
## Warning: st_centroid assumes attributes are constant over geometries
ggplot(powiats3,aes(fill=pred,col=pred)) +
  geom_sf(show.legend = FALSE) +
  scale_fill_viridis(option = "H") + 
  facet_wrap(~Month,ncol=4)+
  theme_gdocs()

From the plots above I see, that in May-October the PM2,5 pollution is minimal.

Smoothing Seasonal Effects

Here I treat months as ‘pseudo-spatial’ areal entities. I assume that months are correlated with the neighbouring one (eg. pollution in January is correlated with pollution in both February and December). At first I create a list reflecting the correlation.

m_nabs <- list(January=c(12,2),February=c(1,3),March=c(2,4),April=c(3,5),
               May=c(4,6),June=c(5,7),July=c(6,8),August=c(7,9),
               September=c(8,10),October=c(9,11),November=c(10,12),December=c(11,1))

Then I fits the fourth model

gp4 <- gam(Pollution~s(Month,bs='mrf',xt=list(nb=m_nabs),k=4)+
             s(X,Y,k=38,m=2,bs='gp'),data=pm25_tibble_long)

And compare the AIC with model number 2.

AIC(gp2,gp4)
##           df      AIC
## gp2 44.39189 3585.067
## gp4 30.17615 4275.409

AIC still suggests that the second model is better but let’s look at the monthly pattern anyway. It may provide some new information about pollution patterns.

At first let’s focus on a randomly choosen powiat and check monthly effects there.

predset <- tibble(jpt_kod_je=1812,
                  Month=factor(months,
                               levels=months))
# Unifying data type
predset$jpt_kod_je <- as.numeric(predset$jpt_kod_je)
POW$jpt_kod_je <- as.numeric(POW$jpt_kod_je)

predset <- POW %>% st_centroid() %>% 
  st_coordinates() %>% as_tibble() %>% 
  bind_cols(POW) %>% 
  inner_join(predset) %>% 
  mutate(pred=predict(gp4,newdata=.))
## Warning: st_centroid assumes attributes are constant over geometries
## Joining with `by = join_by(jpt_kod_je)`
# Plotting
ggplot(predset,aes(x=Month,y=pred)) +
  geom_col()+
  theme_gdocs()

Simple spatial models and summary

Up to this point I presented different ways how to model values known for a point and project them to the whole country. I also presented how save those estimates for a smaller entity - powiat. As I’ve pointed it makes it easier to understand maps. For instance, it is more informative to see that Pm2.5 pollution is high in powiat wodzisławski than to say that it is high in the South near the border with Czechia.

Here I will show that saving values for each entity allows also to create a regression model. The data from statistical offices are usually available for some administrative entity: powiat, county, region and so on.

Stations intensity examination

At first I’ll create a simple OLS model explaining “stations intensity” by population density and estimated pollution.

model = lm(data=POW, formula = stations~predicted_gp+log(pop_density))
summary(model)
## 
## Call:
## lm(formula = stations ~ predicted_gp + log(pop_density), data = POW)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.517e-09 -3.409e-10 -3.754e-11  2.689e-10  1.946e-09 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -2.480e-09  2.010e-10  -12.34   <2e-16 ***
## predicted_gp      1.279e-10  1.084e-11   11.80   <2e-16 ***
## log(pop_density)  2.650e-10  2.426e-11   10.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.676e-10 on 377 degrees of freedom
## Multiple R-squared:  0.4765, Adjusted R-squared:  0.4737 
## F-statistic: 171.6 on 2 and 377 DF,  p-value: < 2.2e-16

It seems that the location of stations are determined by the population density (where more people live, there are more people wanting to have an air quality station nearby) and pollution itself (people in polluted areas want to know how polluted air is).

Further one may like to run spatial models like Spatial Durbin Model, Spatial Lag Model, Geographically Weighted Regression or others.

Lung cancer ailments analysis

Here I experiment with lung cancer analysis. I don’t want to analyse the determinants of lung cancer (which is mostly caused by smoking and other similar behaviours) but to offer additional dimension of how to use processed point data.

At first I add a column with names matching to powiats’ names from cancer data

POW$jpt_nazwa_nowa <- sub(POW$jpt_nazwa_, pattern="powiat ", replacement="")

Then I read lung cancer ailment data. The data comes from the Polish National Cancer Registry https://onkologia.org.pl/pl/raporty

cancer <- read.csv("Data/Krajowy_rejestr_nowotworów/ailments_map_2019_powiaty.csv")

Now I add TERYT code to powiat’s, because names can overlap (cancer data offers only names, not unique TERYT code.) I create a vector with duplicate names:

duplicate_names <- c("bielski","brzeski","grodziski",
                     "krośnieński","nowodworski","opolski",
                     "ostrowski","średzki","świdnicki","tomaszowski")

Then, I create a dataframe with TERYT codes for all powiats.

TERYT_POW <- POW %>% st_drop_geometry()

Now I create a dataframe with TERYT code if it’s powiat with unique name, force NA otherwise.

cancer2 <- cancer %>% filter(Region %ni% duplicate_names)
cancer2 <- left_join(cancer2,TERYT_POW, join_by(Region==jpt_nazwa_nowa))
cancer2 <- left_join(cancer, cancer2, by = c("Region", "Współczynnik.surowy"))

At last I must manually place TERYT codes based on Współczynnik.surowy in df and on the map on https://onkologia.org.pl/pl/raporty. If they match, I fill the dataframe with corresponding TERYT code.

cancer2[7,3] <- 1405
cancer2[178,3] <- 3005
cancer2[9,3] <- 0618
cancer2[276, 3] <- 1016
cancer2[11,3] <- 1202
cancer2[16,3] <- 1601
cancer2[46,3] <- 1807
cancer2[128,3] <- 0802
cancer2[52,3] <- 3025
cancer2[366,3] <- 0218
cancer2[76,3] <- 2003
cancer2[265,3] <- 2402
cancer2[328,3] <- 1609
cancer2[121,3] <- 0612
cancer2[157,3] <- 1416
cancer2[209,3] <- 3017
cancer2[162,3] <- 2210
cancer2[376,3] <- 1414
cancer2[169,3] <- 0617
cancer2[206,3] <- 0219

Finally, I trim the dataframe and merge cancer data to POW

cancer2 <- cancer2[,2:4]
cancer2[is.na(cancer2$jpt_nazwa_)==T,"jpt_nazwa_"]=paste("powiat",cancer2[is.na(cancer2$jpt_nazwa_)==T,"Region"],paste="")
cancer2$jpt_kod_je <- as.numeric(cancer2$jpt_kod_je)
POW$jpt_kod_je <- as.numeric(POW$jpt_kod_je)
POW <- left_join(POW, cancer2, by = join_by(jpt_kod_je))
POW <- POW[,-10]
colnames(POW)[2] <- "jpt_nazwa_"

It is time then to run a simple OLS model.

model2 <- lm(POW, formula=Współczynnik.surowy~predicted_January+log(pop_density))
summary(model2)
## 
## Call:
## lm(formula = Współczynnik.surowy ~ predicted_January + log(pop_density), 
##     data = POW)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.279 -12.443  -1.702  10.193 109.782 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        63.3822     6.8047   9.315  < 2e-16 ***
## predicted_January  -1.3022     0.3845  -3.387 0.000781 ***
## log(pop_density)    3.6700     0.7828   4.688 3.86e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.07 on 377 degrees of freedom
## Multiple R-squared:  0.0655, Adjusted R-squared:  0.06055 
## F-statistic: 13.21 on 2 and 377 DF,  p-value: 2.843e-06

The outcome is surprising - the coefficient of pollution is negative. It means that model omits some important information like healthcare availability. Also, such a simple regression disregards possible influence of population density on pollution. However, this OLS regression is clearly shows other application of data obtained from points.

The end

What I hoped to show above is that it is possible to obtain data for every administrative entity, having only values for specific points. It is possible to gain insights just from the location of the points, and of course from values at those points.

There are various ways of estimating and predicting values like pollution, therefore every decent researcher should carefully account for the use of one or other method.