---
title: "Marginal Effects of Elements in RBGlass1 Data"
author: "Matthew Harris <<me@me.com>>"
date: "`r Sys.Date()`"
output:
  tufte::tufte_html: default
liftr:
  from: "rocker/r-base:latest"
  maintainer: "Matthew Harris"
  email: "me@me.com"
  cran:
    - tufte
    - tidyverse
    - archdata
    - ggeffects
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


## Libraries
```{r libraries, message=FALSE, warning=FALSE}
library("tidyverse")
library("archdata")
library("ggeffects")
```


## Data
These data are from the `archdata` package and details the concentrations for 11 major and minor elements in 105 Romano-British waste glass specimens from two furnace sites, Leicester and Mancetter. 

```{r}
data(RBGlass1)
rbg <- RBGlass1 
summary(rbg)
```

## Exploratory Analysis
Across the eleven elements in this data set, it appears the `Al`, `Fe`, and `P` have pretty good discrimination between the two furnaces.

```{r plot_hist, fig.height=6, fig.width=8, message=FALSE, warning=FALSE, paged.print=FALSE}
ggplot(gather(rbg,element,value,-Site), aes(x=value, fill=Site)) +
  geom_histogram(bins = 15, position = "fill") +
  facet_wrap(~element, scales = "free") +
  theme_bw()
```

## Model
Fit a binomial generalized linear model to these data with the `glm` function.
```{r fit_glm}
glm1 <- glm(Site ~ ., data = rbg, family=binomial(link='logit'))
summary(glm1)
```

## Evaluate
Compute the marginal effect for six of the eleven elements
```{r marginals}
marginals <- tibble(elements = c("Fe","P","Ti","Al","Pb","Mn"), 
                    fit = list(glm1)) %>%
  mutate(marginal = purrr::map2(fit, elements, ggpredict)) %>%
  select(-fit) %>%
  unnest()
```

## Plot

```{r plot_marginals, fig.width=8, fig.height=6}
ggplot(marginals, aes(x = x, y = predicted, fill = elements)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .25) +
  theme_bw() +
  facet_wrap(~elements, scales = "free")

```

## Session information

The R session information for compiling this document is shown below.

```{r session}
sessionInfo()
```