Recently I was looking at a Kaggle dataset on air quality and the relationships between the concentrations of contaminants and environmental variables such as wind speed and direction, temperature, humidity, time of day and time of year for 2017. The idea was to develop a state space model for forecasting the level of particulate matter in the atmosphere. Using this model we’ll know if it safe to go outside or best to shut the windows and wait for the smog to pass. I also used Bayesian methods to estimate the probability of the the next day being over the threshold of air pollution standards.
Data cleaning
Firstly we’ll import the dataset and change the column names to something a little more usable.
suppressPackageStartupMessages(library(lubridate))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(dlm))
suppressPackageStartupMessages(library(RColorBrewer))
suppressPackageStartupMessages(library(corrplot))
suppressPackageStartupMessages(library(mice))
suppressPackageStartupMessages(library(gridExtra))
suppressPackageStartupMessages(library(forecast))
suppressPackageStartupMessages(library(MCMCpack))
# import data set
air <- read.csv("Earlwood_Air_Data_17_18.csv")
# cleaning up the columns of this dataset because they're horrible
colnames(air) <- c("date", "time", "wind_dir", "temp", "wind_speed", "no", "no2", "co", "ozone", "ozone4hr", "pm10", "pm2.5", "humidity", "sd1", "co8hr")
# set colours
colRamp <- colorRampPalette(c("purple2", "green2"))(100)
myCol <- rgb(22/255, 207/255, 193/255)
myCol2 <- rgb(239/255, 185/255, 63/255)
head(air)
## date time wind_dir temp wind_speed no no2 co ozone ozone4hr pm10 ## 1 01/01/2017 01:00 152.3 22.6 0.4 0 0.4 NA 2.0 2.1 23.6 ## 2 01/01/2017 02:00 134.0 22.6 0.3 NA NA NA NA 2.2 21.0 ## 3 01/01/2017 03:00 132.2 22.6 0.3 0 0.6 NA 1.7 2.0 20.0 ## 4 01/01/2017 04:00 125.8 22.7 0.2 0 0.5 NA 1.7 1.8 21.4 ## 5 01/01/2017 05:00 107.9 22.8 0.6 0 0.3 NA 2.1 1.8 21.5 ## 6 01/01/2017 06:00 107.5 23.0 0.5 0 0.3 NA 2.1 1.9 23.5 ## pm2.5 humidity sd1 co8hr ## 1 7.0 87.2 49.01 NA ## 2 6.6 87.2 46.56 NA ## 3 7.2 87.0 47.40 NA ## 4 7.1 87.2 53.70 NA ## 5 4.3 86.8 41.65 NA ## 6 8.6 85.7 40.73 NA
On inspection there appears to be quite a few missing values. These need to be treated in a sensible way otherwise we will end up with not so sensible results. We’ll use the MICE package for investigating and treating the missing values.
md.pattern(air)
## date time humidity temp wind_dir wind_speed sd1 pm10 pm2.5 no no2 ## 7670 1 1 1 1 1 1 1 1 1 1 1 ## 3 1 1 1 1 1 1 1 1 1 1 1 ## 57 1 1 1 1 1 1 1 1 1 1 1 ## 24 1 1 1 1 1 1 1 0 1 1 1 ## 108 1 1 1 1 1 1 1 1 0 1 1 ## 27 1 1 1 1 1 1 1 1 1 0 0 ## 401 1 1 1 1 1 1 1 1 1 1 1 ## 3 1 1 1 1 1 1 1 0 1 1 1 ## 19 1 1 1 1 1 1 1 0 0 1 1 ## 346 1 1 1 1 1 1 1 1 1 0 0 ## 3 1 1 1 1 1 1 1 1 1 0 0 ## 2 1 1 1 1 1 1 1 0 1 0 0 ## 4 1 1 1 1 1 1 1 1 0 1 1 ## 2 1 1 1 1 1 1 1 0 0 1 1 ## 1 1 1 1 1 0 0 0 1 1 1 1 ## 70 1 1 1 1 1 1 1 1 1 0 0 ## 3 1 1 1 1 1 1 1 0 1 0 0 ## 1 1 1 1 1 1 1 1 0 0 1 1 ## 1 1 1 1 1 0 0 0 0 1 1 1 ## 1 1 1 1 1 1 1 1 0 1 0 0 ## 1 1 1 1 1 1 1 1 1 0 0 0 ## 6 1 1 1 1 1 1 1 0 0 0 0 ## 1 1 1 1 1 1 1 1 0 0 0 0 ## 6 1 1 1 1 1 1 1 0 0 0 0 ## 1 1 1 1 1 0 0 0 1 1 0 0 ## 1 1 1 0 0 1 1 1 0 1 0 0 ## 1 1 1 1 0 0 0 0 0 0 1 1 ## 4 1 1 0 0 0 0 0 0 0 0 0 ## 17 1 1 0 0 0 0 0 0 0 0 0 ## 0 0 22 23 25 25 25 92 170 489 489 ## ozone4hr ozone co co8hr ## 7670 1 1 0 0 2 ## 3 1 0 0 0 3 ## 57 0 1 0 0 3 ## 24 1 1 0 0 3 ## 108 1 1 0 0 3 ## 27 1 1 0 0 4 ## 401 0 0 0 0 4 ## 3 0 1 0 0 4 ## 19 1 1 0 0 4 ## 346 1 0 0 0 5 ## 3 0 1 0 0 5 ## 2 1 1 0 0 5 ## 4 0 0 0 0 5 ## 2 0 1 0 0 5 ## 1 1 1 0 0 5 ## 70 0 0 0 0 6 ## 3 1 0 0 0 6 ## 1 0 0 0 0 6 ## 1 1 1 0 0 6 ## 1 0 0 0 0 7 ## 1 0 0 0 0 7 ## 6 1 0 0 0 7 ## 1 0 1 0 0 7 ## 6 0 0 0 0 8 ## 1 1 0 0 0 8 ## 1 0 0 0 0 9 ## 1 1 0 0 0 9 ## 4 1 0 0 0 12 ## 17 0 0 0 0 13 ## 568 866 8784 8784 20362
The “co” and “co8hr” columns are completely empty so these will be removed. The ozone and NO columns are also fairly patchy, thsee will be treated using multiple imutation however we’ll look at the correlation structure of the full dataset first. If we look at the correlation structure after imputation there will dependencies in the data. Given there are only 10% of the values missing for ozone it isn’t likely to be a problem but it’s good practice to be mindful of the structure of your data.
# remove co and co8hr
remove_list <- c("co", "co8hr")
air <- air[,!(match.fun("%in%")(colnames(air), remove_list))]
# retain the complete records
air_complete <- air[complete.cases(air),]
dim(air_complete)
## [1] 7670 13
# correlation plot
cor.mat <- cor(air_complete[,-(1:2)])
par(mfrow=c(1,2))
corrplot(cor.mat, method="ellipse", tl.col="black", tl.cex=0.7)
corrplot(cor.mat, method="number", tl.col="black", number.cex=0.7, tl.cex=0.7)
The correlation structure in this data doesn’t seem to be particularly strong, the strongest being wind speed and sd1 at -0.63 and ozone and temperature at 0.62. We’re mostly interested in the correlations with 
and 
which aren’t strong at all.
Using a multiple imputation technique we’ll treat the missing values.
airimp <- mice(air, meth="pmm", maxit = 5, m = 3)
## ## iter imp variable ## 1 1 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1 ## 1 2 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1 ## 1 3 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1 ## 2 1 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1 ## 2 2 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1 ## 2 3 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1 ## 3 1 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1 ## 3 2 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1 ## 3 3 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1 ## 4 1 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1 ## 4 2 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1 ## 4 3 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1 ## 5 1 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1 ## 5 2 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1 ## 5 3 wind_dir temp wind_speed no no2 ozone ozone4hr pm10 pm2.5 humidity sd1
densityplot(airimp)
The plots above show the density of the imputed values (red) with respect to the observedv(blue). If the assumption is the data is missing at random or missing completely at random we would expect to see the densities to be similar. In most part we do however, wind direction and humidity are a little different to the observed. This may require further investigation to find a more suitable imputation method. The other variables seem ok. We’ll come back to address wind direction and humidity.
# keeping all the imputed records
airc <- complete(airimp)
Fix the date and time into a POSIXct structure.
# fixing the date and time fields
airc$date <- as.Date(airc$date, format="%d/%m/%Y")
airc$dt <- ymd_hms(paste0(airc$date, " ", airc$time, ":00"))
Data exploration
Firstly it is important to see how the variables are distributed. These will essentially be those plotted above but it will be good to see them in slightly more detail. The ozone 4 hr window variable will be dropped now that the finer level ozone (1 hr) has been imputed.
var_list <- colnames(airc)[!(match.fun("%in%")(colnames(airc), c("date", "time", "ozone4hr", "dt")))]
density_plots <- function(var) ggplot(airc, aes_string(x = var)) + geom_histogram(fill = "darkcyan", col = "grey20")
plot_list <- lapply(var_list, density_plots)
grid.arrange(grobs=plot_list, nrow=3, ncol = 4, as.table=FALSE)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
We want to visualise this data as a time series and identify any interesting patterns over time. Given this data only extends for 1 year seasonal trends will not be able to be picked up however, there may be day to day cycles we can observe such as temperature fluctuations. Since ozone and temperature showed the highest correlation I’ll start my search there.
grid.arrange(ggplot(airc, aes(x=dt, y=temp)) + geom_line(),
ggplot(airc, aes(x=dt, y=ozone)) + geom_line(), nrow=2)
There is a lot of detail in these plots. To reduce the granularity to gain a better picture I could aggregate to the day level by either taking the max or mean of the values for each day, or simply choose a smaller window to observe such as a single month. The first option will eliminate much of the noise in the data but it will also eliminate the daily cycle of temperatures increasing throughout the day until midday, then decreasing. I expect to see this same pattern in the ozone data since the two variables show some positive correlation, abliet weak. Given the hottest day last year occurred in February I’ll focus my attention there (as good of a reason as any).
# create a crude measure of light intensity based on the hour of the day for plotting purposes
airc$sunlight <- 1-cos(2*pi/24*hour(airc$dt))
colRampSun <- colorRampPalette(c("darkblue", "orange2", "yellow3"))(100)
# select feb and plot
feb <- match.fun("==")(month(airc$dt), 2)
grid.arrange(
ggplot(airc[feb,], aes(x=dt, y=temp, col=sunlight)) + geom_line() + scale_color_gradientn(colours=colRampSun) + theme_dark(),
ggplot(airc[feb,], aes(x=dt, y=ozone, col=sunlight)) + geom_line() + scale_color_gradientn(colours=colRampSun) + theme_dark(),
nrow=2)
When viewing an individual month it is easy to see the relationship between these two features. The colour scale indicates the time of day at a crude level where the darkest colour occurs at midnight and the lightest at midday. This all makes intuitive sense, as the sun goes down the temperatures drop.
Air pollution
A quick Google search returned particulate matter 10 micrometers () or less in diameter and 2.5 micrometers (PM2.5) or less in diameter are standardised meassures of air pollution which form environmental policies, for example the Australian governemnt has set national ambient air quality standard of 50 micrograms for PM10 averaged over a 24 hour period. A recent study has shown that an increase of 10 
on the previous day is associated with a 0.6% increase in the risk of mortality. The air quality index (AQI) is a standardised measure to advise on the health risk of outdoor activities and calculated as
![\[ \text{AQI}=100\frac{\text{pm}_{10}}{50} \]](https://gradientdescending.com/wp-content/ql-cache/quicklatex.com-9ec6398aa4b2fe480592b0b6f1e7791f_l3.png "Rendered by QuickLaTeX.com")
This provides us with some good direction. We will build a state space model to forecast the levels for the next day and for any given day compute the probability of experiencing an average of greater than 50 PM10 levels (or an AQI > 1) using Bayesian simulation.
First of all lets look at the series of for 2017 and for Feb.
grid.arrange(
ggplot(airc, aes(x=dt, y=pm10, col=sunlight)) + geom_line() + theme_dark() + scale_color_gradientn(colours = colRampSun) + labs(title="PM10 for 2017"),
ggplot(airc[feb,], aes(x=dt, y=pm10, col=sunlight)) + geom_line() + theme_dark() + scale_color_gradientn(colours = colRampSun) + labs(title="PM10 for February"),
nrow=2)
There are some large spikes scattered throughout the year which are greater than 50 although these values are recorded at the hour. To control some of these spikes it may be worth taking the square root transform as there appears to be some skewness. There are also some negative values in the data which cause some problem. In this instance we’ll replace those values with the smallest non-zero value in the data (also, I expect that these are incorrect values anyway and should probably imputed with something more reasonable).
min_pm10 <- min(airc$pm10[airc$pm10 > 0])
airc <- tbl_df(airc) %>%
mutate(sq_pm10 = ifelse(pm10 <= 0, sqrt(min_pm10), sqrt(pm10)))
## Warning in sqrt(pm10): NaNs produced
grid.arrange(
ggplot(airc, aes(x=pm10)) + geom_histogram(fill="darkcyan", col = "grey20") + labs(title="PM10"),
ggplot(airc, aes(x=sq_pm10)) + geom_histogram(fill="darkcyan", col = "grey20") + labs(title="log_PM10"),
nrow=1
)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`. ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
To be consistent with the air pollution standards KPI these need to be averaged over a 24 hours period. A 24 hr moving average will be applied to the series.
air_day <- tbl_df(airc) %>%
mutate(pm10_ma = ma(pm10, 24),
sq_pm10_ma = ma(sq_pm10, 24),
temp_ma = ma(temp, 24),
humidity_ma = ma(humidity, 24),
wind_dir_ma = ma(wind_dir, 24),
wind_speed_ma = ma(wind_speed, 24),
ozone_ma = ma(ozone, 24),
no2_ma = ma(no2, 24))
grid.arrange(
ggplot(airc, aes(x=dt, y=sq_pm10, col=sunlight)) + geom_line() + theme_dark() + scale_color_gradientn(colours = colRampSun) + labs(title="sqrt(PM10) for 2017"),
ggplot(airc[feb,], aes(x=dt, y=sq_pm10, col=sunlight)) + geom_line() + theme_dark() + scale_color_gradientn(colours = colRampSun) + labs(title="sqrt(PM10) for February"),
nrow=2
)
grid.arrange(
ggplot(air_day, aes(x=dt, y=sq_pm10_ma, col=sunlight)) + geom_line() + theme_dark() + scale_color_gradientn(colours = colRampSun) + labs(title="sqrt(PM10) for 2017 (moving average over 24 hours)"),
ggplot(air_day[feb,], aes(x=dt, y=sq_pm10_ma, col=sunlight)) + geom_line() + theme_dark() + scale_color_gradientn(colours = colRampSun) + labs(title="sqrt(PM10) for February (moving average over 24 hrs)"),
nrow=2)
## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
## Warning: Removed 24 rows containing missing values (geom_path).
Observing the series for all of 2017 it still appears to be quite noisey without any clear cyclic patter With a 24 hour moving average applied. Only 1 point of the PM10 levels excede 50 (or 7 on the square root scale). It would be intersting to find out what caused this spike, a bushfire/burn off perhaps?.
State space model
Given the data is quite noisey and doesn’t have any clear trend a suitable model will be a local linear model with a regression component using temperature, wind speed, wind direction and humidty. A Kalman filter will be applied to the series to smooth out the noise and to forecast for the next day(s).
The raw data has a clear daily cycle in which case a seasonal component would be suitable. Given the data points are sampled hourly this would mean 25 seasonal factors to estimate and that would take a while to run. For the time being the model will be fit without seasonal components. Given there is temperature has this cyclic component and is incorporated in the regression part of the model this is still coming into play to some degree.
# function for fitting the state space model
pm10_ssm <- function(dates, nAhead, par, mle_opt){
# date ids
dates_id_m <- match.fun("%in%")(air_day$date, dates) # only the dates for modelling
dates_id_f <- match.fun("%in%")(air_day$date, seq(min(dates), max(dates)+nAhead, by = "day")) # the dates for modelling + the forecasting dates
# set up the x matrix for the regression component
air_x <- air_day[dates_id_m, c("humidity", "temp", "wind_speed", "wind_dir")]
air_x_f <- air_day[dates_id_f, c("humidity", "temp", "wind_speed", "wind_dir")]
air_y <- air_day$sq_pm10[dates_id_m]
air_y_f <- c(air_day$sq_pm10[dates_id_m], rep(NA, sum(dates_id_f)-sum(dates_id_m)))
air_y_t <- air_day$sq_pm10[dates_id_f]
# state space model
ssm <- function(parm, X) dlmModPoly(1, dV=parm[1]) + dlmModReg(X=X, dV=parm[2], dW=parm[3:7])
# estimate the paramters using MLE. can be slow for large amounts of data
if(mle_opt){
mle <- dlmMLE(y=air_y, parm=par, X=air_x, build=ssm)
fit <- ssm(parm=mle$par, X=air_x_f)
par <- mle$par
}else{
fit <- fit <- ssm(parm=par, X=air_x_f)
}
# fit model and apply kalman filter
filtered <- dlmFilter(air_y_f, fit)
# create plot
df <- data.frame(dt = c(air_day$dt[dates_id_f], air_day$dt[dates_id_f]),
pm10 = c(air_y_t, filtered$f),
series = c(rep("observed", sum(dates_id_f)), rep("filtered", sum(dates_id_f))))
g <- ggplot(df, aes(x=dt, y=pm10, col=series)) +
geom_line() +
scale_color_manual(values = c("orange2", "darkblue")) +
geom_vline(xintercept = ymd_hms(paste(max(dates), "23:00:00")), col = "red") +
labs(x = "Date", y = "sqrt PM10 levels", title = "Kalman filter of PM10 levels + forecast from red line") +
scale_y_continuous(limits = c(0, max(df$pm10))) +
theme_dark()
return(list(filtered = filtered, par = par, fit = fit, g = g))
}
mod_output <- pm10_ssm(dates = seq(as.Date("2017-02-01"), as.Date("2017-02-28"), by = "day"),
nAhead = 3,
par = c(1, 1, 1, 0, 0.83, 1, 0),
mle_opt = TRUE)
mod_output$g
The above graph has been forecast using the observations for February and forecasting the first 3 days in March. The observations look very similar to what you’d expect from a random walk and this is essentially what the SSM picks up on. With systems like this there is a lot of uncertainty in the data and going back further in time won’t add anything to the model, instead it’s the short term observations that will drive most of the forecasting power.
mod_output <- pm10_ssm(dates = seq(as.Date("2017-02-26"), as.Date("2017-02-28"), by = "day"),
nAhead = 3,
par = rep(1, 7),
mle_opt = TRUE)
mod_output$g
While it’s still higher than the observerd it is picking up on the features in the series a little more.
mod_output <- pm10_ssm(dates = seq(as.Date("2017-03-17"), as.Date("2017-03-19"), by = "day"),
nAhead = 3,
par = rep(1, 7),
mle_opt = TRUE)
mod_output$g
Plotting a random time period shows that the model is definitely picking up on some features in the data. This appears to be doing quite well forecasting the next 3 days. While the first 1.5 days are fairly flat it does pick up on the dip near the end of the 21st. This will be information will be gained from the regression component of the model. A key feature it doesn’t pick up on is the large spike in the evening of the 22nd. This appears to be an outlier and would require a far more sophisticated model than the one here. Ordinarily outliers are removed or treated but that largely depends on the context. Here we would want to ensure we are capturing the likelihood of observing such spikes given we are talking about peoples health. This could be a separate model on it’s own.
If a 24 hr moving average is applied to the forecasted series I expect the forecasted and observed would be quite close. Feel free to take this function and plot for more time frames.
Bayesian Estimation
When there is a lot of uncertianty in complex systems Bayesian approaches work quite well. A Bayesian framework can be applied to state space models however, for demonstration a simple Bayesian regression model will be fit on the data using the MCMCpack package. The previous days environmental observations will be used to estimate the posterior distribution of levels for the following day. This will enable the calculation of the probability a level greater than 
will be observed in the following day. For simplicity the dataset will be averaged to the daily level, that way we don’t have to untangle the correlations by taking the moving average. This is a good approximation though, once built the observed moving average values can be applied to the model to prodice the predictive posterior distribution as the variances should be approximately the same.
air_daily_avg <- tbl_df(airc) %>%
mutate(yr = year(dt), mth = month(dt), day = day(dt)) %>%
mutate(dt = as.Date(paste0(yr, "-", mth, "-", day), format="%Y-%m-%d")) %>%
group_by(dt) %>%
summarise(temp = mean(temp),
humidity = mean(humidity),
wind_speed = mean(wind_speed),
wind_dir = mean(wind_dir),
ozone = mean(ozone),
no2 = mean(no2),
sq_pm10 = sqrt(mean(pm10)))
On the averaged dataset set the dependent variable to be at lag 1 i.e. one day ahead.
# set sq_pm10 to be one day ahead
#air_daily_avg <- air_day
air_daily_avg$lagged_sq_pm10 <- c(air_daily_avg$sq_pm10[2:nrow(air_daily_avg)], NA)
The model formula will include temp, wind speed and humidity. Wind direction wasn’t as strong as the other 3. There are some interesting relationships between ozone and no2 however I want to keep this model reduced to simple environmental observations at this time.
# fit bayesian model
form <- as.formula(lagged_sq_pm10 ~ temp + wind_speed + humidity)
mcmc_mod <- MCMCregress(form, data = air_daily_avg)
summary(mcmc_mod)
## ## Iterations = 1001:11000 ## Thinning interval = 1 ## Number of chains = 1 ## Sample size per chain = 10000 ## ## 1. Empirical mean and standard deviation for each variable, ## plus standard error of the mean: ## ## Mean SD Naive SE Time-series SE ## (Intercept) 4.67532 0.310424 3.104e-03 3.104e-03 ## temp 0.04864 0.008019 8.019e-05 8.019e-05 ## wind_speed -0.37186 0.083304 8.330e-04 8.154e-04 ## humidity -0.01323 0.003473 3.473e-05 3.473e-05 ## sigma2 0.50790 0.037907 3.791e-04 3.907e-04 ## ## 2. Quantiles for each variable: ## ## 2.5% 25% 50% 75% 97.5% ## (Intercept) 4.07870 4.46487 4.67709 4.88181 5.301783 ## temp 0.03287 0.04327 0.04863 0.05401 0.064450 ## wind_speed -0.53596 -0.42589 -0.37163 -0.31657 -0.210449 ## humidity -0.02024 -0.01551 -0.01325 -0.01090 -0.006509 ## sigma2 0.43915 0.48128 0.50591 0.53225 0.587171
par(mfrow=c(2,2))
plot(mcmc_mod, auto.layout=FALSE, smooth=TRUE)
All environmental variables are quite strong predictors. The posterior distributions of the parameters are now used to simulate the predictive posterior distribution of for the following day. To test this a day will be randomly chosen and treated as a new observation.
mcmc_x <- as.matrix(mcmc_mod)[,1:(length(form)+1)]
sigma <- as.matrix(mcmc_mod)[,(length(form)+2)]
# choose a day and pretend it's new data
day <- sample(1:nrow(air_daily_avg), 1)
newdata <- as.matrix(cbind(1, as.data.frame(air_daily_avg[day, c("temp", "wind_speed", "humidity")])))
# simulate
z <- mcmc_x %*% t(newdata) + rnorm(10000, 0, sigma)
# plot
ggplot(data.frame(pred_sq_pm10 = z[,1]), aes(x=pred_sq_pm10)) + geom_histogram(fill="darkcyan", col = "black") +
geom_vline(xintercept=air_daily_avg$lagged_sq_pm10[day], col = "red")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
mean(z > sqrt(50))
## [1] 0
newdata
## 1 temp wind_speed humidity ## [1,] 1 23.69583 1.383333 69.55417
The predictive posterior distribution is 
where 
4.3877499 and 
0.5172252. The probability that we will obsereve levels greater than with the given observations is 0. Out of 10000 simulations there were 0 that were greater than the threshold, effectively saying it’s extrememly unlikely. This isn’t too surprising. The data shows only one instances where the levels exceded the threshold and this could be an outlier.
For demonstration lets create “the perfect storm”.
newdata <- as.matrix(cbind(1, data.frame(temp = 37, wind_speed = 0, humidity = 25)))
# simulate
z <- mcmc_x %*% t(newdata) + rnorm(10000, 0, sigma)
# plot
ggplot(data.frame(pred_sq_pm10 = z[,1]), aes(x=pred_sq_pm10)) + geom_histogram(fill="darkcyan", col = "black")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
mean(z > sqrt(50))
## [1] 0.0586
newdata
## 1 temp wind_speed humidity ## [1,] 1 37 0 25
Under these values there is a 5.9% chance of observing levels greater than . I’m sure with more observations and a more complex model this will change to be a higher probability however based on 2017 data this is what we’re seeing. This Bayesian model could be improved by applying suitable priors to better reflect the prior beliefs about the system.
Conclusion
We have fit a basic state space model to forecast the next few days of levels based on historical observations and other environmental observations and a Bayesian model to estimate the probablity of observing 
on 24 hr averaged data. This is an initial basic analysis of the data. Next steps include:
- Applying to PM2.5
- Fitting a seasonal component to the SSM
- Incorporating more covariates into the SSM and Bayesian model
- Adding suitable priors to the Bayesain model
- Better imputation of missing values
- Sourcing more historical data
