---
title: "BayesGP: COVID-19 Example"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{BayesGP: COVID-19 Example}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
oldpar <- par(no.readonly = TRUE)  # Save current settings 
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>", fig.height=3, fig.width=5, margins=TRUE
)
knitr::knit_hooks$set(margins = function(before, options, envir) {
  if (!before) return()
  graphics::par(mar = c(1.5 + 0.9, 1.5 + 0.9, 0.2, 0.2), mgp = c(1.45, 0.45, 0), cex = 1.25, bty='n')
})
```


```{r setup}
library(BayesGP)
```

# COVID-19 Example

## Data and Model

We will illustrate the use of BayesGP using the `covid_canada` dataset, which contains the daily death count of COVID-19 in Canada.

```{r}
head(covid_canada)
```

For simplicity, let's consider the following model:
\[Y_i|\lambda_i \sim \text{Poisson}(\lambda_i)\]
\[\log(\lambda_i) = \mathbf{x}_i^T\boldsymbol{\beta} + f(t_i)\]
where $\mathbf{x}_i$ denotes the fixed effect of weekdays, and $f$ is an unknown function to be inferred.

To make inference of the unknown function $f$, we use the $\text{IWP}_3(\sigma)$ model:
\[\frac{\partial^p{f}(t)}{\partial t^p} = \sigma \xi(t),\]
with the boundary (initial) conditions that $\frac{\partial^q{f}(0)}{\partial t^q} = 0$ for all $0\leq q <p$.
Here $\xi(t)$ is the standard Gaussian white noise process, or can be viewed as the distributional derivative of the standard Brownian motion.

## Inference

To fit the above model using the BayesGP package, we simply do the following:
```{r warning=FALSE}
fit_result <- model_fit(new_deaths ~ weekdays1 + weekdays2 + weekdays3 + weekdays4 + weekdays5 + weekdays6 +
                          f(smoothing_var = t, model = "IWP", order = 3, k = 100, sd.prior = list(prior = "exp", param = list(u = 0.02, alpha = 0.5), h = 1)), 
                        data = covid_canada, method = "aghq", family = "Poisson")
```

We can take a look at the posterior summary of this model:
```{r}
summary(fit_result)
```


We can also see the inferred function $f$:
```{r warning=FALSE}
plot(fit_result)
```

We can use the `predict` function to obtain the posterior summary of $f$ or its derivative at `new_data`:

For the function $f$:
```{r warning=FALSE}
predict_f <- predict(fit_result, variable = "t", newdata = data.frame(t = seq(from = 605, to = 615, by = 0.1)))
matplot(x = predict_f[,1], y = predict_f[,2:4], type = 'l', ylab = "f", xlab = "t", col = c("grey", "grey", "black"), lty = c("dashed", "dashed", "solid"))
```


For the first derivative:
```{r warning=FALSE}
predict_f1st <- predict(fit_result, variable = "t", newdata = data.frame(t = seq(from = 605, to = 615, by = 0.1)), deriv = 1)
matplot(x = predict_f1st[,1], y = predict_f1st[,2:4], type = 'l', ylab = "f'", xlab = "t", col = c("grey", "grey", "black"), lty = c("dashed", "dashed", "solid"))
```


For the second derivative:
```{r warning=FALSE}
predict_f2nd <- predict(fit_result, variable = "t", newdata = data.frame(t = seq(from = 605, to = 617, by = 0.1)), deriv = 2)
matplot(x = predict_f2nd[,1], y = predict_f2nd[,2:4], type = 'l', ylab = "f''", xlab = "t", col = c("grey", "grey", "black"), lty = c("dashed", "dashed", "solid"))
par(oldpar)
```