Skip to content

Commit

Permalink
update GLM slides
Browse files Browse the repository at this point in the history
  • Loading branch information
Pakillo committed Nov 27, 2023
1 parent 5f027ae commit 7530ff0
Show file tree
Hide file tree
Showing 4 changed files with 63 additions and 46 deletions.
19 changes: 10 additions & 9 deletions glm_binomial.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -227,8 +227,8 @@ plot(factor(survived) ~ factor(class), data = titanic)

\footnotesize

```{r echo=2:3}
library(ggmosaic)
```{r echo=TRUE, out.width="80%"}
library("ggmosaic")
ggplot(titanic) +
geom_mosaic(aes(x = product(survived, class))) +
labs(x = "", y = "Survived")
Expand Down Expand Up @@ -566,16 +566,15 @@ simulateResiduals(tit.sex, plot = TRUE)

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

```{r echo=FALSE}
Sex is a confounder

```{r echo=FALSE, out.width="70%"}
library("dagitty")
g1 <- dagitty("dag {
Class -> Survival
Sex -> Survival
Sex -> Class
}")
coordinates(g1) <- list(
x = c(Class = 1, Sex = 1, Survival = 2),
y = c(Class = 0, Sex = 2, Survival = 1)
)
plot(g1)
```

Expand Down Expand Up @@ -629,14 +628,16 @@ summary(tit.sex.class.int)

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

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


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

```{r tit_sex_class_effects2, echo=FALSE}
\footnotesize

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

Expand Down
Binary file modified glm_binomial.pdf
Binary file not shown.
90 changes: 53 additions & 37 deletions glm_count.Rmd
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
---
title: "GLM for count data: Poisson regression"
title: "GLM for count data"
author: "Francisco Rodríguez-Sánchez"
institute: "https://frodriguezsanchez.net"
aspectratio: 43 # use 169 for wide format
Expand Down Expand Up @@ -149,31 +149,6 @@ seedl.glm <- glm(count ~ light, data = seedl, family = poisson)
summary(seedl.glm)
```

\normalsize


## Parameter estimates are in log scale!

Parameter estimates (log scale):
```{r poisson_params, echo=T}
coef(seedl.glm)[1]
```

\vspace{5mm}

**We need to back-transform**: apply the inverse of the logarithm

```{r echo=T}
exp(coef(seedl.glm)[1])
```


## Using effects package

```{r echo=2}
library(effects)
allEffects(seedl.glm)
```


## Estimated distribution of the **slope** parameter
Expand Down Expand Up @@ -203,15 +178,41 @@ plot(simulate_parameters(seedl.glm)) +



## So what's the relationship between Nseedlings and light?

```{r poisson_effects, echo=2}
#allEffects(seedl.glm)
plot(allEffects(seedl.glm))
## Parameter estimates are in log scale!

Parameter estimates (log scale):
```{r poisson_params, echo=T}
coef(seedl.glm)[1]
```

\vspace{5mm}

**We need to back-transform**: apply the inverse of the logarithm

```{r echo=T}
exp(coef(seedl.glm)[1])
```


## Using easystats

\footnotesize

```{r message = FALSE, echo=TRUE}
library("easystats")
parameters(seedl.glm)
parameters(seedl.glm, exponentiate = TRUE)
```

## How Nseedlings decrease with light

```{r}
estimate_relation(seedl.glm)
```


## Using visreg
## Visualising how Nseedlings decrease with light

```{r poisson_visreg, echo = 2:3}
library(visreg)
Expand All @@ -220,14 +221,16 @@ points(count ~ light, data = seedl, pch = 20)
```

::: hide :::

## Using sjPlot

```{r echo=4}
```{r echo=4, eval=FALSE}
library(sjPlot)
library(ggplot2)
theme_set(theme_minimal(base_size = 16))
sjPlot::plot_model(seedl.glm, type = "eff", show.data = TRUE)
```

:::


Expand Down Expand Up @@ -374,9 +377,9 @@ summary(seedl.overdisp)

\footnotesize

```{r poisson_overdisp2, echo=TRUE}
allEffects(seedl.overdisp)
allEffects(seedl.glm)
```{r poisson_overdisp2, echo=TRUE, message = FALSE}
parameters(seedl.overdisp)
parameters(seedl.glm)
```


Expand All @@ -386,13 +389,13 @@ allEffects(seedl.glm)
:::::::::::::: {.columns align=center}
::: {.column width="50%"}
```{r pois_overdisp_eff1, echo=FALSE}
plot(allEffects(seedl.overdisp))
visreg(seedl.overdisp, scale = "response")
```
:::

::: {.column width="50%" }
```{r pois_overdisp_eff2, echo=FALSE}
plot(allEffects(seedl.glm))
visreg(seedl.glm, scale = "response")
```
:::
::::::::::::::
Expand Down Expand Up @@ -475,6 +478,19 @@ predict(seedl.glm, newdata = new.lights, type = "response", se.fit = TRUE)
```


## Prediction (easystats)

\footnotesize

```{r echo=TRUE}
new.lights <- data.frame(light = c(10, 90))
estimate_expectation(seedl.glm, data = new.lights)
```

```{r echo=TRUE}
estimate_prediction(seedl.glm, data = new.lights)
```



## Poisson GLM: more examples
Expand Down
Binary file modified glm_count.pdf
Binary file not shown.

0 comments on commit 7530ff0

Please sign in to comment.