```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,
                      message = FALSE,
                      warning = FALSE,
                      dev = "svg")
```

```{r packages, message=FALSE, warning=FALSE}
devtools::load_all() # if using the rproject dowloaded from the slides
# source("utils-glm.R") # if using a standard setup
library(here)
library(tidyr) # for data manipulation
library(dplyr) # for data manipulation
library(ggplot2) # plotting
library(car) # general utilities
library(effects) # for extracting and plotting effects 
library(emmeans) # for marginal means
```

```{r options, include = FALSE}
theme_set(theme_minimal(base_size = 15))
```

```{r script, echo=FALSE}
## Link in Github repo
downloadthis::download_link(
  link = download_link("https://github.com/stat-teaching/SMDA-2023/blob/master/labs/lab11.R"),
  button_label = "Download the R script",
  button_type = "danger",
  has_icon = TRUE,
  icon = "fa fa-save",
  self_contained = FALSE
)
```

```{r loading-data}
dat <- read.csv(here("data", "nwords.csv"))
```

# Overview

This dataset `nwords.csv` represent a developmental psychology study investigating the factors that influence children language development. The dataset includes information about the number of words in a a language task, and some potential predictors: caregiver behavior, socioeconomic status, amount of time spent with the child during a week and the presence of a baby-sitter or not.

- `child`: numeric variable representing the child ID number.
- `timebonding`: numeric variable representing the average hours per week of child-parent activities
- `nwords`: numeric variable representing the number of words produced by the child in a language task.
- `caregiving`: numeric variable representing the caregiver's behavior in a parent-child interaction task, measured on a scale from 0 to 10
- `ses`: categorical variable representing the family socioeconomic status, with three levels: "Low", "Middle", and "High"

1. Importing data and overall check
    - think about factors levels, scale of the numerical variables, etc.
2. Exploratory data analysis of predictors and the relationships between predictors and the number of words
3. Model fitting with `glm()` and `poisson` family starting from additive models and then adding the interactions.
4. Model diagnostic of the chosen model
    - overdispersion
    - residuals
    - outliers and influential points
5. Interpreting the effects and plotting
6. Fit a quasi-poisson and negative binomial version of the chosen model, and check how parameters are affected

# 1. Importing data and overall check

```{r, eval = FALSE}
dat <- read.csv("data/nwords.csv")
```

```{r}
str(dat)
```

Check for `NA` values:

```{r}
count_na(dat) # from utils-glm.R

# equivalent to sapply(dat, function(x) sum(is.na(x)))
```

Everything seems good, we do not have `NA` values.

Let's convert categorical variables into factor setting the appropriate order:

- `ses`: low, middle, high
- `babysitter`: no, yes

```{r}
dat$ses <- factor(dat$ses, levels = c("low", "middle", "high"))
dat$babysitter <- factor(dat$babysitter, levels = c("no", "yes"))

levels(dat$ses)
levels(dat$babysitter)
```

`timebonding` and `caregiving` are the two numerical predictors. Given that we are going to fit main effects and interactions and we want to interpret the intercept and test the interaction on a meaningful point, we create two centered versions of these variables:

```{r}
dat$timebonding0 <- dat$timebonding - mean(dat$timebonding)
dat$caregiving0 <- dat$caregiving - mean(dat$caregiving)
```

Then we can see what happen to the variables:

```{r}
cols <- grepl("timebonding|caregiving", names(dat))
lapply(dat[, cols], function(x) round(c(mean = mean(x), sd = sd(x)), 3))
```

# 2. Exploratory data analysis

```{r}
summary(dat)
```

Let's do some plotting of predictors:

```{r}
par(mfrow = c(2,2))
hist(dat$timebonding)
hist(dat$caregiving)
barplot(table(dat$babysitter))
barplot(table(dat$ses))
```

Comments?

Also the response variable:

```{r}
hist(dat$nwords, probability = TRUE)
lines(density(dat$nwords), col = "salmon", lwd = 2)
```

Let's plot the theoretical distributions, Poisson and Gaussian:

```{r}
m <- mean(dat$nwords)
s <- sd(dat$nwords)
xs <- seq(0, 50, 1)

hist(dat$nwords, probability = TRUE, xlim = c(-10, 50), ylim = c(0, 0.1))
curve(dnorm(x, m, s), add = TRUE, col = "green", lwd = 2)
lines(xs, dpois(xs, m), col = "red", lwd = 2)
lines(xs, dpois(xs, m), col = "red", lwd = 2)
```

Comments?

Let's plot some bivariate distributions:

```{r}
# caregiving ~ timebonding
r <- cor(dat$timebonding, dat$caregiving)
plot(dat$timebonding, dat$caregiving, pch = 19)
abline(lm(dat$caregiving ~ dat$timebonding), col = "red", lwd = 2)
text(30, 2, label = paste("r =", round(r, 2)))
```

```{r}
par(mfrow = c(1,2))
boxplot(caregiving ~ ses, data = dat)
boxplot(timebonding ~ ses, data = dat)
```

```{r}
mosaicplot(table(dat$ses, dat$babysitter))
# or
barplot(table(dat$babysitter, dat$ses), col = c("red", "green"), beside = TRUE)
legend(7,45, legend = c("no", "yes"), fill = c("red", "green"))
```

```{r}
par(mfrow = c(1,2))
boxplot(caregiving ~ babysitter, data = dat)
boxplot(timebonding ~ babysitter, data = dat)
```

Then the relationships with the `nwords` variable:

```{r}
par(mfrow = c(2,2))
plot(dat$timebonding, dat$nwords, pch = 19)
plot(dat$caregiving, dat$nwords, pch = 19)
boxplot(dat$nwords ~ dat$ses)
boxplot(dat$nwords ~ dat$babysitter)
```

Comments?

Finally some interactions plot (with `lm`):

```{r, fig.width=10, fig.height=10}
par(mfrow = c(2,1))
colors <- c(low = "red", middle = "blue", high = "green")
plot(dat$timebonding, dat$nwords, col = colors[dat$ses], pch = 19)
lms <- lapply(split(dat, dat$ses), function(x) lm(nwords ~ timebonding, data = x))
lapply(1:length(lms), function(i) abline(lms[[i]], col = colors[i], lwd = 2))

plot(dat$caregiving, dat$nwords, col = colors[dat$ses], pch = 19)
lms <- lapply(split(dat, dat$ses), function(x) lm(nwords ~ caregiving, data = x))
lapply(1:length(lms), function(i) abline(lms[[i]], col = colors[i], lwd = 2))
```

```{r, fig.width=10, fig.height=10}
par(mfrow = c(2,1))
colors <- c(no = "orange", yes = "purple")

plot(dat$timebonding, dat$nwords, col = colors[dat$babysitter], pch = 19)
lms <- lapply(split(dat, dat$babysitter), function(x) lm(nwords ~ timebonding, data = x))
lapply(1:length(lms), function(i) abline(lms[[i]], col = colors[i], lwd = 2))

plot(dat$caregiving, dat$nwords, col = colors[dat$babysitter], pch = 19)
lms <- lapply(split(dat, dat$babysitter), function(x) lm(nwords ~ caregiving, data = x))
lapply(1:length(lms), function(i) abline(lms[[i]], col = colors[i], lwd = 2))
```

# 3. Model fitting with `glm()` and `poisson`

Let's start by using an additive model:

```{r}
fit <- glm(nwords ~ timebonding + caregiving + babysitter + ses, family = poisson(link = "log"), data = dat)
summary(fit)
```

And always plotting before anything else:

```{r}
plot(allEffects(fit))
```

Comments? How could you describe the results? Something different from the descriptive statistics?

```{r fig.width=10, fig.height=10}
car::residualPlots(fit)
```

Comments? Are we missing something?

```{r}
dat |> 
    ggplot(aes(x = timebonding, y = nwords, color = ses)) +
    geom_point() +
    stat_smooth(method = "glm", method.args = list(family = poisson()), se = FALSE)
```

```{r}
dat |> 
    ggplot(aes(x = caregiving, y = nwords, color = ses)) +
    geom_point() +
    stat_smooth(method = "glm", method.args = list(family = poisson()), se = FALSE)
```

```{r}
dat |> 
    ggplot(aes(x = timebonding, y = nwords, color = babysitter)) +
    geom_point() +
    stat_smooth(method = "glm", method.args = list(family = poisson()), se = FALSE)
```

```{r}
dat |> 
    ggplot(aes(x = caregiving, y = nwords, color = babysitter)) +
    geom_point() +
    stat_smooth(method = "glm", method.args = list(family = poisson()), se = FALSE)
```

Let's add some interactions:

```{r}
fit2 <- glm(nwords ~ timebonding*ses + caregiving*ses + timebonding*babysitter + caregiving*babysitter, family = poisson(link = "log"), data = dat)
summary(fit2)
```

```{r fig.width=10, fig.height=10}
plot(allEffects(fit2))
```

# 4. Model fitting with `MASS::glm.nb()`

There is still evidence for overdispersion, even after including all predictors and a series of interactions. Let's assume that this is our most complex model, we need to take into account the overdispersion:

```{r}
performance::check_overdispersion(fit2)

fit3 <- MASS::glm.nb(nwords ~ timebonding*ses + caregiving*ses + timebonding*babysitter + caregiving*babysitter, data = dat)

summary(fit3)
```

Now overdispersion is taken into account and standard errors are larger:

```{r}
car::compareCoefs(fit2, fit3)

# test statistics for the poisson model

data.frame(
    poisson = fit2$coefficients/sqrt(diag(vcov(fit2))),
    negative_binomial = fit3$coefficients/sqrt(diag(vcov(fit3)))
)
```










