---
title: "Lung cancer risk in Pennsylvania, USA, 2002"
output: flexdashboard::flex_dashboard
---
  
  
```{r, echo=FALSE, results='hide'}
library(SpatialEpi)
library(ggplot2)
library(leaflet)
library(DT)
library(spdep)
library(sp)
library(dplyr)
library(INLA)

# Data and map
data(pennLC)
map <- pennLC$spatial.polygon

# Observed cases
d <- group_by(pennLC$data, county) %>% summarize(Y = sum(cases))
# Expected cases
pennLC$data <- pennLC$data[order(pennLC$data$county, 
                                 pennLC$data$race,
                                 pennLC$data$gender,
                                 pennLC$data$age), ]
E <- expected(population = pennLC$data$population,
              cases = pennLC$data$cases,
              n.strata = 16)
d$E <- E[match(d$county, unique(pennLC$data$county))]
# Smoker proportions
d <- merge(d, pennLC$smoking, by = "county")
# SMRs
d$SMR <- d$Y/d$E

# Add data to map
rownames(d) <- d$county
map <- SpatialPolygonsDataFrame(map, d, match.ID = TRUE)

# INLA
nb <- poly2nb(map)
nb2INLA("map.adj", nb)
g <- inla.read.graph(filename = "map.adj")
map$re_u <- 1:nrow(map@data)
map$re_v <- 1:nrow(map@data)
formula <- Y ~ smoking + f(re_u, model = "besag", graph = g, scale.model = TRUE) + f(re_v, model = "iid")
res <- inla(formula, family = "poisson", data = map@data,
            E = E, control.predictor = list(compute = TRUE))

# Add results to map
map$RR <- res$summary.fitted.values[, "mean"]
map$LL <- res$summary.fitted.values[, "0.025quant"]
map$UL <- res$summary.fitted.values[, "0.975quant"]

# Plot coefficient covariate
# ggplot object plot_covariate
marginal <- inla.smarginal(res$marginals.fixed$smoking)
marginal <- data.frame(marginal)
plot_covariate <- ggplot(marginal, aes(x = x, y = y)) +
  geom_line() +
  labs(x = expression(beta[1]), y = "Density") +
  geom_vline(xintercept = 0, col = "blue") +
  theme_bw()

# Map risk
# leaflet object leaflet_map
pal <- colorNumeric(palette = "YlOrRd", domain = map$RR)

labels <- sprintf("<strong> %s </strong> <br/> Observed: %s <br/> Expected: %s <br/>
                  Smokers proportion: %s <br/>SMR: %s <br/>RR: %s (%s, %s)",
                  map$county, map$Y,  round(map$E, 2),  map$smoking, round(map$SMR, 2),
                  round(map$RR, 2), round(map$LL, 2), round(map$UL, 2)) %>% lapply(htmltools::HTML)

leaflet_map <- leaflet(map) %>% addTiles() %>%
  addPolygons(color = "grey", weight = 1, fillColor = ~pal(RR), fillOpacity = 0.5,
              highlightOptions = highlightOptions(weight = 4), label = labels, 
              labelOptions = labelOptions(style = list("font-weight" = "normal", padding = "3px 8px"),
                                          textsize = "15px", direction = "auto")) %>%
  addLegend(pal = pal, values = ~RR, opacity = 0.5, title = "RR", position = "bottomright")
```

  
Column {data-width=600}
-------------------------------------
  
### Map
  
```{r}
leaflet_map
```

Column {data-width=400}
-------------------------------------
  
### Smoking effect
  
```{r}
plot_covariate
```   

### Table

```{r}
datatable(map@data)
```