glm_binomial.Rmd

---
title: "Generalised Linear Models"
subtitle: "Logistic regression"
author: "Francisco Rodríguez-Sánchez"
institute: "https://frodriguezsanchez.net"
aspectratio: 43  # use 169 for wide format
fontsize: 10pt
output: 
  binb::metropolis:
    keep_tex: no
    incremental: yes
    fig_caption: no
    pandoc_args: ["--lua-filter=hideslide.lua"]
urlcolor: blue
linkcolor: blue
header-includes:
  - \definecolor{shadecolor}{RGB}{230,230,230}
  # - \setbeamercolor{frametitle}{bg=black}
---



```{r knitr_setup, include=FALSE, cache=FALSE}

library("knitr")

### Chunk options ###

## Text results
opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE, size = 'tiny')

## Code decoration
opts_chunk$set(tidy = FALSE, comment = NA, highlight = TRUE, prompt = FALSE, crop = TRUE)

# ## Cache
# opts_chunk$set(cache = TRUE, cache.path = "knitr_output/cache/")

# ## Plots
# opts_chunk$set(fig.path = "knitr_output/figures/")
opts_chunk$set(fig.align = 'center', out.width = '90%')

### Hooks ###
## Crop plot margins
knit_hooks$set(crop = hook_pdfcrop)

## Reduce font size
## use tinycode = TRUE as chunk option to reduce code font size
# see http://stackoverflow.com/a/39961605
knit_hooks$set(tinycode = function(before, options, envir) {
  if (before) return(paste0("\n \\", options$size, "\n\n"))
  else return("\n\n \\normalsize \n")
  })

```


## Q: Survival of passengers on the Titanic ~ Class

Read `titanic_long.csv` dataset and fit linear model (survival ~ class).

\vspace{5mm}

```{r prepare_titanic_data, echo=FALSE, eval=FALSE}
titanic <- read.table("http://www.amstat.org/publications/jse/datasets/titanic.dat.txt")
names(titanic) <- c("class", "age", "sex", "survived")
titanic$class <- factor(titanic$class, labels = c("crew", "first", "second", "third"))
titanic$age <- factor(titanic$age, labels = c("child", "adult"))
titanic$sex <- factor(titanic$sex, labels = c("female", "male"))
write.csv(titanic, file = "data/titanic_long.csv", row.names=FALSE, quote=FALSE)
```

```{r read_titanic, echo=FALSE}
titanic <- read.csv("data/titanic_long.csv")
head(titanic)
```


## Quiz: Did passenger class influence survival?

https://pollev.com/franciscorod726


## Let's check linear model:

\footnotesize

```{r titanic_lm, echo=TRUE}
m5 <- lm(survived ~ class, data = titanic)
library("easystats")
check_model(m5)
```



## Weird residuals!

```{r titanic_lm_resid, echo=FALSE}
hist(resid(m5))
```


## What if your residuals are clearly non-normal \newline or variance not constant (heteroscedasticity)?

\Large

- Binary variables (0/1)

- Counts (0, 1, 2, 3, ...)

- Categories ("small", "medium", "large"...)

\vspace{5mm}

- **Generalised Linear Models to the rescue!**


## Generalised Linear Models

1. **Response variable** - distribution `family`
    + Bernouilli - Binomial
    + Poisson
    + Gamma
    + etc

\vspace{3mm}

2. **Predictors** (continuous or categorical)

\vspace{3mm}

3. **Link function**
    + Gaussian: identity
    + Binomial: logit, probit
    + Poisson: log...
    + See [`family`](http://www.rdocumentation.org/packages/stats/functions/family).



## The modelling process

```{r out.width="40%"}
include_graphics("images/modeling_process.png")
```

Bolker 2008



## Bernouilli - Binomial distribution (Logistic regression) 

Response variable: **Yes/No** (e.g. survival, sex, presence/absence)

Canonical link function: `logit` (*log odds*), but others possible (see `family`)

$$
\begin{aligned} 
logit(p) = \log \left( \dfrac {p} {1-p}\right) \\ 
\end{aligned} 
$$

Then

$$
\begin{aligned} 
logit(P(alive)) = a + bx \\  
P(alive) = invlogit(a + bx) = \dfrac {e^{a+bx}} {1+e^{a+bx}} \\  
\end{aligned} 
$$

## Where is the variance?

In a Gaussian GLM

$$
y \sim Normal(\mu, \sigma)
$$

In a Binomial GLM

$$
y \sim Binomial(n, p)
$$
`n` = number of trials

`p` = probability of success

$$
{Var}(y)=np(1-p)
$$
(maximum variance when `p` around 0.5)



# Back to survival of Titanic passengers 

## How many survived in each class?

\vspace{3mm}

```{r echo=TRUE}
table(titanic$class, titanic$survived)
```



## How many survived in each class? (*dplyr*)

\vspace{3mm}

```{r titanic_dplyr, echo=c(-1)}
library("dplyr")
titanic %>%
  group_by(class, survived) %>%
  summarise(count = n())
```




## Data visualisation (mosaic plot)

```{r titanic_eda, echo=TRUE, out.width = "100%"}
plot(factor(survived) ~ factor(class), data = titanic)
```

## Mosaic plots (ggplot2)

\footnotesize

```{r echo=TRUE, out.width="80%"}
library("ggmosaic")
ggplot(titanic) +
  geom_mosaic(aes(x = product(survived, class))) +
  labs(x = "", y = "Survived")
```



## Fitting GLMs in R: `glm`

```{r echo=TRUE}
tit.glm <- glm(survived ~ class, 
               data = titanic, 
               family = binomial)
```

\vspace{3mm}

which corresponds to 

$$
\begin{aligned}  
logit(P(survival)_{i}) = a + b \cdot class_{i} \\  
logit(P(survival)_{i}) = a + b_{first} + c_{second} + d_{third} \\
\end{aligned}  
$$



## Interpreting binomial GLM

\scriptsize
```{r titanic_glm, echo=1}
tit.glm <- glm(survived ~ class, data = titanic, family = binomial)
summary(tit.glm)
```


## Binomial GLM estimates are in `logit` scale!

We need to **back-transform** (apply *inverse logit*):

  - Manually: `plogis`

  - Automatically: `easystats`, etc.
  

::: hide :::

## Interpreting logistic regression output (`effects` pkg)

```{r tit_glm_effects, echo=TRUE}
library("effects")
allEffects(tit.glm)
```


## Interpreting logistic regression output (`effects` pkg)

Including confidence intervals:

\footnotesize

```{r tit_glm_effects2, echo=TRUE}
summary(allEffects(tit.glm))
```

:::

## Interpreting logistic regression output (`easystats`)

```{r echo=TRUE}
library("easystats")   # 'modelbased' pkg
estimate_means(tit.glm)
```


## Analysing differences among factor levels (class)

\footnotesize

```{r echo=TRUE}
estimate_contrasts(tit.glm)
```


## Pseudo R-squared for GLMs

```{r echo=TRUE}
library("easystats")   # 'performance' pkg
r2(tit.glm)
```

\vspace{5mm}

But there are caveats (e.g. see [here](https://stats.stackexchange.com/questions/3559/which-pseudo-r2-measure-is-the-one-to-report-for-logistic-regression-cox-s) and [here](http://data.library.virginia.edu/is-r-squared-useless/))




## Presenting model results

```{r echo=TRUE}
kable(xtable::xtable(tit.glm), digits = 2)
```



## Presenting model results

\footnotesize

```{r echo=TRUE}
library("modelsummary")
modelsummary(tit.glm, output = "markdown")
```


## Visualising model: `effects` package

```{r effects_plot, echo=TRUE}
plot(allEffects(tit.glm))
```


## Visualising model: `visreg` package

```{r echo=2}
library("visreg")
visreg(tit.glm, scale = "response", rug = FALSE)
```



## Visualising model: `sjPlot` package

```{r echo=4, out.width="70%", eval=FALSE}
library(ggplot2)
library(sjPlot)
theme_set(theme_minimal(base_size = 16))
sjPlot::plot_model(tit.glm, type = "eff")
```


## Visualising model: `easystats` (`see` package)

```{r echo=2}
library("easystats")
plot(parameters(tit.glm), show_intercept = TRUE)
```


## Model without intercept

```{r echo=TRUE, out.width="70%"}
no.intercept <- glm(survived ~ class - 1, family = binomial, data = titanic)
plot(parameters(no.intercept))
```


# Model checking

## `plot(model)` not very useful with binomial GLM

```{r tit_glm_check, echo=2, tinycode = TRUE}
layout(matrix(1:4, nrow = 2))
plot(tit.glm)
dev.off()
```


::: hide :::

## Binned residual plots for logistic regression

```{r echo=TRUE, out.width="80%", eval=FALSE}
plot(binned_residuals(tit.glm))
```



```{r binnedplot, eval=FALSE}
predvals <- predict(tit.glm, type="response")
arm::binnedplot(predvals, titanic$survived - predvals)
```

:::

## check_model (easystats)

```{r echo=TRUE, eval=FALSE}
check_model(tit.glm)
```

::: hide :::

## Posterior predictive checking

Simulate data from fitted model (**yrep**) and compare with observed data (**y**)

```{r out.width="70%", echo = 2, tinycode = TRUE}
library("easystats")
check_predictions(tit.glm)
```


```{r eval=FALSE}
library(bayesplot)
sims <- simulate(tit.glm, nsim = 100)
ppc_bars(titanic$survived, yrep = t(as.matrix(sims)))
```

:::


## Residual diagnostics with DHARMa

\footnotesize

```{r echo=TRUE}
library("DHARMa")
simulateResiduals(tit.glm, plot = TRUE)
```

See https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html


## Calibration plot

Compares predicted vs observed probabilities (grouped by quantiles)

\footnotesize

```{r echo=TRUE, out.width="60%"}
library("predtools")
titanic$surv.pred <- predict(tit.glm, type = "response")
calibration_plot(data = titanic, obs = "survived", pred = "surv.pred", 
                 x_lim = c(0,1), y_lim = c(0,1)) 
```

---

Passenger class was important, but lots of unexplained variation

```{r echo=FALSE}
visreg(tit.glm, scale = "response")
```


---

The goal is not to test whether the model's assumptions are "true", because all models are false.

Rather, the goal is to assess exactly **how the model fails to describe the data**, as a path towards **model comprehension, revision, and improvement**.

Richard McElreath. _Statistical Rethinking_



## Recapitulating


1. **Visualise data**

\vspace{2mm}

2. **Fit model**: `glm`. Don't forget to specify `family`!

\vspace{2mm}

3. **Examine model**: `summary`

\vspace{2mm}

4. **Back-transform parameters** from *logit* into probability scale (`estimate_means`)

\vspace{2mm}

5. **Plot model**: `visreg`, ...

\vspace{2mm}

6. **Check model**: `check_model`, `DHARMa::simulateResiduals`, `calibration_plot`





# Q: Did men have higher survival than women?

## Quiz

https://pollev.com/franciscorod726


## First, visualise data

```{r tit_sex_eda, out.width="100%"}
plot(factor(survived) ~ as.factor(sex), data = titanic)
```

## Fit model

\footnotesize

```{r tit_sex, echo=FALSE}
tit.sex <- glm(survived ~ sex, data = titanic, family = binomial)
summary(tit.sex)
```


## Model interpretation

\footnotesize

```{r tit_sex_effects, echo=TRUE}
estimate_means(tit.sex)
```

```{r echo=FALSE, out.width="50%"}
visreg(tit.sex, scale = "response", rug = FALSE)
```





## Model checking

```{r echo=TRUE}
simulateResiduals(tit.sex, plot = TRUE)
```


# Q: Did women have higher survival because they travelled more in first class?

## Did women have higher survival because they travelled more in first class?

Sex is a confounder

```{r echo=FALSE, out.width="70%"}
library("dagitty")
g1 <- dagitty("dag {
              Class -> Survival
              Sex -> Survival
              Sex -> Class
              }")
plot(g1)
```



## Let's look at the data

\footnotesize

```{r tit_women, echo=TRUE}
table(titanic$class, titanic$survived, titanic$sex)
```




## Quiz

https://pollev.com/franciscorod726


## Fit additive model with both factors 

\footnotesize

```{r echo=FALSE}
tit.sex.class.add <- glm(survived ~ class + sex, family = binomial, data = titanic)
summary(tit.sex.class.add)
```


## Plot additive model

```{r echo=2}
par(mfcol = c(1, 2))
visreg(tit.sex.class.add, scale = "response", rug = FALSE)
dev.off()
```



## Fit model with the interaction of both factors

\footnotesize

```{r tit_sex_class, echo=FALSE}
tit.sex.class.int <- glm(survived ~ class * sex, family = binomial, data = titanic)
summary(tit.sex.class.int)
```


## Women had higher survival than men, even within the same class

```{r tit_sex_class_effects, echo=TRUE}
estimate_means(tit.sex.class.int)
```


## Women had higher survival than men, even within the same class

\footnotesize

```{r tit_sex_class_effects2, echo=TRUE}
visreg(tit.sex.class.int, by = "sex", xvar = "class", scale = "response", rug = FALSE)
```




## Visualising model (`sjPlot`)

```{r echo=TRUE, out.width="80%"}
library("sjPlot")
plot_model(tit.sex.class.int, type = "int")
```


## Comparing models

\scriptsize

```{r echo=T}
library("easystats")   # 'performance' pkg
compare_performance(tit.sex.class.add, tit.sex.class.int)
```


## Comparing parameters

\footnotesize

```{r echo=T}
compare_parameters(tit.sex.class.add, tit.sex.class.int)
```




```{r echo=FALSE, eval = FALSE}
## Calibration plot
titanic$surv.pred <- predict(tit.sex.class.int, type = "response")
calibration_plot(data = titanic, obs = "survived", pred = "surv.pred", 
                 x_lim = c(0,1), y_lim = c(0,1), nTiles = 10)
```






## Extra exercises:

\Large

Is survival related to age?

Are age effects dependent on sex?



# Logistic regression for proportion data


## Read Titanic data in different format

Read `titanic_prop.csv` data.

\vspace{3mm}

```{r read_tit_short, echo = FALSE}
tit.prop <- read.csv("data/titanic_prop.csv")
head(tit.prop)
```

\vspace{3mm}

These are the same data, but summarized (see `Freq` variable).


## Use `cbind(n.success, n.failures)` as response

\footnotesize

```{r binom_prop, echo=1}
prop.glm <- glm(cbind(Yes, No) ~ Class, data = tit.prop, family = binomial)
summary(prop.glm)
```


::: hide :::

## Effects

```{r prop_glm_effects, echo=FALSE}
allEffects(prop.glm)
```

\vspace{5mm}

**Compare with former model based on binary data:**
\vspace{3mm}

```{r comp, echo=FALSE}
allEffects(tit.glm)
```

\vspace{3mm}

Same results!

:::


## Survival probability by class

```{r echo=TRUE}
estimate_means(prop.glm)
```





# Logistic regression with continuous predictors


----

Example dataset: [GDP and infant mortality](http://vincentarelbundock.github.io/Rdatasets/doc/car/UN.html)

\vspace{3mm}

Read `UN_GDP_infantmortality.csv`.

\vspace{3mm}

```{r read_gdp, echo = FALSE}
#gdp <- read.csv("http://vincentarelbundock.github.io/Rdatasets/csv/car/UN.csv")
gdp <- read.csv("data/UN_GDP_infantmortality.csv")
names(gdp) <- c("country", "mortality", "gdp")
summary(gdp)
```


## Q: Is infant mortality related to GDP?

https://pollev.com/franciscorod726


## Visualising data

```{r gdp_eda}
plot(mortality ~ gdp, data = gdp, main = "Infant mortality (per 1000 births)")
```


## Fit model

\footnotesize

```{r gdp_glm, echo=1}
gdp.glm <- glm(cbind(mortality, 1000 - mortality) ~ gdp, 
               data = gdp, family = binomial)
summary(gdp.glm)
```


::: hide :::

## Effects

```{r gdp_effects, echo=T}
allEffects(gdp.glm)
```

## Effects plot

```{r gdp_effectsplot}
plot(allEffects(gdp.glm))
```

:::



## Plot model using visreg:

```{r gdp_visreg, echo=c(2,3)}
library(visreg)
visreg(gdp.glm, scale = "response")
points(mortality/1000 ~ gdp, data = gdp)
```


## Residuals diagnostics with DHARMa

```{r echo=TRUE}
simulateResiduals(gdp.glm, plot = TRUE)
```




# Overdispersion

---

**Overdispersion**:

more variation in the data than assumed by statistical model

$$
{Var}(y)=np(1-p)
$$


## Testing for overdispersion (DHARMa)

```{r echo = TRUE}
simres <- simulateResiduals(gdp.glm, refit = TRUE)
testDispersion(simres, plot = FALSE)
```


---

`quasibinomial` allows us to model overdispersed binomial data


## Overdispersion in logistic regression with proportion data

\footnotesize 

```{r logreg_overdisp, echo=1}
gdp.overdisp <- glm(cbind(mortality, 1000 - mortality) ~ gdp, 
               data = gdp, family = quasibinomial)
summary(gdp.overdisp)
```


## Mean estimates do not change after accounting for overdispersion

But standard errors (uncertainty) do!

\scriptsize

```{r logreg_overdisp2, echo=TRUE}
parameters(gdp.overdisp)
```

\vspace{5mm}

```{r echo=T}
parameters(gdp.glm)
```



## But standard errors (uncertainty) do!

:::::::::::::: {.columns align=center}

::: {.column width="50%"}
```{r overdisp_plot1, echo=FALSE, out.width="100%"}
visreg(gdp.glm, scale = "response", main = "Binomial")
points(mortality/1000 ~ gdp, data = gdp, pch = 20)
```
:::

::: {.column width="50%" }
```{r overdisp_plot2, echo=FALSE, out.width="100%"}
visreg(gdp.overdisp, scale = "response", main = "Quasibinomial")
points(mortality/1000 ~ gdp, data = gdp, pch = 20)
```
:::
::::::::::::::



# Think about the shape of relationships


## Think about the shape of relationships

Not everything has to be linear...

```{r}
visreg(gdp.glm, ylab = "Mortality (logit scale)")
```


## Residuals show non-linear pattern

```{r echo=FALSE}
library(ggResidpanel)
resid_panel(gdp.glm)
```


## Calibration plot shows non-linear pattern

```{r echo=FALSE}
gdp.na <- na.omit(gdp)
gdp.na$fit <- fitted(gdp.glm)
plot(gdp.na$fit, gdp.na$mortality, 
     xlab = "Probability of mortality (predicted)",
     ylab = "Mortality (observed)")
```


## Calibration plot shows non-linear pattern

```{r}
gdp.na$mort.pred <- predict(gdp.glm, type = "response")
gdp.na$mort.prob <- gdp.na$mortality/1000
calibration_plot(data = gdp.na, obs = "mort.prob", pred = "mort.pred", 
                 x_lim = c(0,0.1), y_lim = c(0,0.1), nTiles = 10)
```


## Trying polynomial predictor (GDP + GDP^2^)

:::::::::::::: {.columns align=center}

::: {.column width="50%"}
```{r echo=FALSE, out.width="100%"}
visreg(gdp.overdisp, main = "Mortality ~ GDP", ylab = "Mortality (logit scale)")
```
:::

::: {.column width="50%" }
```{r echo=FALSE, out.width="100%"}
gdp.overdisp2 <- glm(cbind(mortality, 1000 - mortality) ~ gdp + I(gdp*gdp), 
               data = gdp, family = quasibinomial)
visreg(gdp.overdisp2, main = "Mortality ~ GDP + GDP^2", ylab = "Mortality (logit scale)")
```
:::
::::::::::::::






## Think about the shape of relationships

:::::::::::::: {.columns align=center}

::: {.column width="50%"}
```{r echo=FALSE, out.width="100%"}
visreg(gdp.overdisp, main = "Mortality ~ GDP", scale = "response", ylab = "Mortality")
points(mortality/1000 ~ gdp, data = gdp, pch = 20)
```
:::

::: {.column width="50%" }
```{r echo=FALSE, out.width="100%"}
gdp.overdisp2 <- glm(cbind(mortality, 1000 - mortality) ~ gdp + I(gdp*gdp), 
               data = gdp, family = quasibinomial)
visreg(gdp.overdisp2, main = "Mortality ~ GDP + GDP^2", scale = "response", ylab = "Mortality")
points(mortality/1000 ~ gdp, data = gdp, pch = 20)
```
:::
::::::::::::::





## Trying log(GDP)

:::::::::::::: {.columns align=center}

::: {.column width="50%"}
```{r echo=FALSE, out.width="100%"}
visreg(gdp.overdisp, main = "Mortality ~ GDP", ylab = "Mortality (logit scale)")
```
:::

::: {.column width="50%" }
```{r echo=FALSE, out.width="100%"}
gdp.overdisp2 <- glm(cbind(mortality, 1000 - mortality) ~ log(gdp), 
               data = gdp, family = quasibinomial)
visreg(gdp.overdisp2, main = "Mortality ~ log(GDP)", ylab = "Mortality (logit scale)")
```
:::
::::::::::::::


## Trying log(GDP)

:::::::::::::: {.columns align=center}

::: {.column width="50%"}
```{r echo=FALSE, out.width="100%"}
visreg(gdp.overdisp, main = "Mortality ~ GDP", scale = "response", ylab = "Mortality")
points(mortality/1000 ~ gdp, data = gdp, pch = 20)
```
:::

::: {.column width="50%" }
```{r echo=FALSE, out.width="100%"}
gdp.overdisp2 <- glm(cbind(mortality, 1000 - mortality) ~ log(gdp), 
               data = gdp, family = quasibinomial)
visreg(gdp.overdisp2, main = "Mortality ~ log(GDP)", scale = "response", ylab = "Mortality")
points(mortality/1000 ~ gdp, data = gdp, pch = 20)
```
:::
::::::::::::::





```{r eval=FALSE, echo=FALSE}
## Trying Poisson
m <- glm(mortality ~ log(gdp), data = gdp, family = quasipoisson)
summary(m)
visreg(m, scale = "response")
points(mortality ~ gdp, data = gdp)

gdp.na <- na.omit(gdp)
gdp.na$fit <- fitted(m)
plot(gdp.na$fit, gdp.na$mortality, 
     xlab = "Probability of mortality (predicted)",
     ylab = "Mortality (observed)")
```




## More examples

- `moth.csv`: Probability of moth predation on trunk trees depending on morph (light/dark) and distance to Liverpool ([Bishop 1972](https://doi.org/10.2307/3513))

- `seedset.csv`: Comparing seed set among plants (Data from [Harder et al. 2011](https://datadryad.org/resource/doi:10.5061/dryad.0vf86nb1.2))

- `soccer.csv`: Probability of scoring penalty depending on goalkeeper's team being ahead, behind or tied ([Roskes et al 2011](https://doi.org/10.1177%2F0956797611418677))


# Moth predation

## The industrial revolution and evolution of dark morphs 

```{r out.width="90%"}
include_graphics("images/moths.jpg")
```

## The data

```{r echo = 1}
moth <- read.csv("data/moth.csv")
head(moth)
```

## Creating new variable: REMAIN

```{r echo = 1}
moth$REMAIN <- moth$PLACED - moth$REMOVED
head(moth)
```

## Did some morph have higher predation overall?

\footnotesize

```{r}
pred.morph <- glm(cbind(REMOVED, REMAIN) ~ MORPH, data = moth, family = binomial)
summary(pred.morph)
```

## Did some morph have higher predation overall?

```{r}
plot(allEffects(pred.morph))
```


## Did predation increase farther from city centre?

\footnotesize

```{r }
pred.dist <- glm(cbind(REMOVED, REMAIN) ~ DISTANCE, data = moth, family = binomial)
summary(pred.dist)
```

## Did predation increase farther from city centre?

```{r }
plot(allEffects(pred.dist))
```

## Did dark morph have lower predation in city & light have lower predation in countryside?

\footnotesize

```{r }
pred.int <- glm(cbind(REMOVED, REMAIN) ~ MORPH * DISTANCE, data = moth, family = binomial)
summary(pred.int)
```


## Did dark morph have lower predation in city & light have lower predation in countryside?

```{r }
plot(allEffects(pred.int))
```

## Model check

```{r echo = TRUE}
simulateResiduals(pred.int, plot = TRUE)
```


# Seed set among plants

## Seed set among plants

```{r echo=FALSE}
include_graphics("images/tomato.jpg")
```


## Seed set among plants

```{r}
seed <- readr::read_csv("data/seedset.csv")
head(seed)
seed$plant <- as.factor(seed$plant)
```

## Questions:

https://pollev.com/franciscorod726
\vspace{5mm}

- Is seed set related to proportion of outcross pollen (pcmass)?

\vspace{5mm}

- Which plant had lower seed set?





## Number of seeds vs Number of ovules

```{r}
plot(seeds ~ ovulecnt, data = seed)
```


## Number of seeds vs Proportion outcross pollen

```{r}
plot(seeds ~ pcmass, data = seed)
```

## Seed set across plants

```{r echo=FALSE}
seedm <- glm(cbind(seeds, ovulecnt - seeds) ~ plant, data = seed, family = binomial)
#summary(seedm)
plot(allEffects(seedm))
```


## Seed set ~ outcross pollen

```{r echo=FALSE}
seedm <- glm(cbind(seeds, ovulecnt - seeds) ~ plant + pcmass, data = seed, family = binomial)
#summary(seedm)
plot(allEffects(seedm))
```


# Probability of scoring penalty

## Data on penalty shots

```{r echo = TRUE}
soccer <- read.csv("data/soccer.csv")
soccer
```

Does probability of scoring penalty depends on match situation?

https://pollev.com/franciscorod726


## Probability of scoring depending on match situation

```{r echo=FALSE}
soccer.mod <- glm(cbind(Scored, Nshots - Scored) ~ GoalkeeperTeam, data = soccer, family = binomial)
visreg(soccer.mod, scale = "response", 
       ylab = "Probability of scoring")
```