---
title: "SOCI 620: Lab 11"
author: "Peter McMahan"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(rethinking)
# this will tell rethinking to use as many cores as your computer has
options(mc.cores = parallel::detectCores())
set.seed(789)
```

# Predicting video game use
We will use another set of variables from the first wave of Add Health (collected in the United States from 1990 to 1991) to look at the time students spent playing video games.

Start by loading the data and looking at the first few rows to get a sense of it.


```{r load_data}
d <- read.csv("https://soci620.netlify.app/data/media.csv")
head(d)
```




In additionn to demographic data on gender, grade, race, ethnicity, and family income, this data lists the number of hours per week the students report spending consuming various media. We will be looking at `hours_games`, which codes students' responses to the question "How many hours a week do you play video or computer games?"

Start by looking at an overall summary variable and seeing how many values are missing.

```{r inspect_data}
summary(d$hours_games)

# how many missing (NA)?
sum( is.na(d$hours_games) )

# proportion NA?
mean( is.na(d$hours_games) )
# less than 0.2% missing
```


```{r inspect_data_detail}
# for discrete data a table can give a quick summary of the values
tg <- table(d$hours_games)
# if a column-based format is easier to read, you can turn
# the table into a martrix
as.matrix(tg)
# What proportion reported playing no games?
mean(d$hours_games==0,na.rm=TRUE)
```


```{r inspect_data_plot}
# plot the table for a nice visualization
plot(tg)
```


So to summarize, about half of the students say they played no video games. A significant number also reported spending a suspiciously large amount of time playing, at least for the 1990s (8 hours per weekday and 16 hours per day on the weekend would 72 hours, but 10 students report spending at least 70 hours per week playing games).


## Intercept-only model

To build a Poisson model, we will start by building a model with just an intercept to get a feel for the approach:

$$
\begin{aligned}
H_i &\sim Pois(\lambda)\\
\log(\lambda) &= \alpha\\
\alpha &\sim Norm(\mbox{??},\mbox{??})
\end{aligned}
$$

Before we estimate this model, we need to figure out a good prior for α by looking at different priors' implications for λ. To get the implied prior for λ, you simply exponentiate the values of the prior for α. Performing this transformation on a normal distribution gives you a new distribution that is already built in to R: the log-normal distribution. This makes it easy to experiment with.


```{r intercept_prior_search1}
# what is the range of lambda that is reasonable in our situation?
# Since there are only 168 hours in a week, we can look at the 
# range [0,170]
x <- seq(0,170,length.out=5000)

# Try alpha ~ Norm(3,10)
# the function `dlnorm` is the density of the log-normal distribution
plot(x,dlnorm(x,3,10), type='l')
```



That's pretty bad. It assumes that the average time spent playing video games is very low. Let's try some more values for the hyperparameters on the same plot


```{r intercept_prior_search2}
# plotting NA gives you a blank plot. You can specify
# xlim and ylim to set the range of the axes
plot(NA,xlim=range(x),ylim=c(0,0.045))

# lines() adds a line plot to the currently active
# plotting device
lines(x,dlnorm(x,3,0.5), lty=1, lwd=3)
lines(x,dlnorm(x,3,1),   lty=2, lwd=3)
lines(x,dlnorm(x,4,0.5), lty=3, lwd=3)
lines(x,dlnorm(x,4,1),   lty=4, lwd=3)

# legend() adds a legend to the plot. It's an awkward
# function but a very useful one.
legend('topright',lty=1:4,lwd=3,legend=c(
    'N(3,0.5)',
    'N(3,1)',
    'N(4,0.5)',
    'N(4,1)'
))
```



Remember, our prior over $\alpha$ should allow for a range of possible values of $\lambda$, and $\lambda$ is the *average* number of hours per week across all of the students. 

If we pick the second line from the above plot ($\mathrm{Norm}(3,1)$), that will give the highest prior likelihood to an average of just over 0 to 50 hours per week. But the 'fat' tail on the right side means that our prior isn't too opinionated, and it will not insist on a low value if the data disagrees with it.

Using that prior, we can estimate our model:


```{r estimate_intercept_only_model, cache=TRUE, output=FALSE}
# Step 1: specify the model
# (NOTE this does not make any reference to our data `d`.)
m1 <- alist(
    hours_games ~ dpois(lambda),
    log(lambda) <- alpha,
    alpha ~ dnorm(3,1)
)

# make a 'slim' copy of the data (just the columns we care about)
d_slim <- d[c('male','hours_games')]
# get complete cases just on the 'slim' data
d_slim <- d_slim[complete.cases(d_slim),]

# estimate the model using ulam on only complete cases
fit1 <- ulam(m1,data=d_slim,chains=4,cores=4,iter=4000, refresh=0)
```


```{r intercept_only_results}
# examine results
precis(fit1,prob=0.95)
# exponentiate the coefficient(s)
exp(coef(fit1))
```


$\exp(\hat\alpha)$ is the mean of the posterior distribution of the average rate of video game consumption by our students. This means that we might predict that an "average" student would play about 2.9 hours of video games per week.

We can see more by extracting samples.

Note, when using `ulam()`, it extracts samples during estimation. If you try to get more samples than it originally gathered (default 1000 ⨉ number of chains), it will give return `NA` for any beyond that limit. To get more you have to give it a value of `iter` that is at least twice the number of samples you want. (We will talk about why this is when we discuss MCMC and other Monte-Carlo methods.)



```{r intercept_only_pplot1}
samp1 <- extract.samples(fit1)
plot(density(samp1$alpha))
```



To get this in terms of actual rates, just exponentiate the posterior sample.


```{r intercept_only_pplot2}
lambda_post <- exp(samp1$alpha)
plot(density(lambda_post))
```



```{r intercept_only_CI}
quantile(lambda_post,c(0.05,0.95))
```

This tells us a little bit more: we are 90% confident that the average rate of game use among these students is between about 2.80 and 2.88 hours per week.

## Adding a covariate: do boys play more?

Let's test the cultural stereotype that adolescent boys spend more time playing games than adolescent girls. We can us a Poisson regression to estimate the difference in weekly rates of video game playing by boys and girls (keeping in mind the problematic way that Add Health collects data on gender).

$$
\begin{aligned}
H_i &\sim Pois(\lambda_i)\\
\log(\lambda_i) &= \alpha + \beta M_i\\
\alpha &\sim Norm(4,1)\\
\beta &\sim Norm(0,0.5)\\
\end{aligned}
$$

It is usally a good idea for your non-intercept coefficients to have a mean of zero. This is because $\exp(0)=1$, so a value of zero for $\beta$ would mean gender has no association.

For the standard deviation on a coefficient prior, you almost always want a pretty narrow distribution when using a Poisson regression. We'll use 0.5 here, but you can experiment with others. Prior predictive plots can be useful here.


```{r gender_model, cache=TRUE}
# build the model
m2 <- alist(
    hours_games ~ dpois(lambda),
    log(lambda) <- alpha + beta*male,
    alpha ~ dnorm(4,1),
    beta ~ dnorm(0,0.5)
)

# fit
fit2 <- ulam(m2,data=d_slim,chains=4,cores=4,iter=4000)
```


```{r gender_model_precis}
# estimates on the scale of the linear model
precis(fit2, prob=0.95)
# point estimates in terms of lambda:
exp(coef(fit2))
```



This tells a very different story. Girls in the sample play an average of about 1.50 hours per week, while boys play about three times more than girls:
$$
\begin{aligned}
\exp(\alpha + \beta) &= \exp(\alpha)\times\exp(\beta)\\
\exp(0.38 + 1.10) &= \exp(0.38)\times\exp(1.10) \\
1.46\times2.99&=4.37
\end{aligned}
$$

Let's compare the posterior distributions of $\lambda$ for boys and girls. (Remeber, these are the posterior distributions for the _average_ number of hours for boys and girls, not the predicted distribution across a sample of boys and girls.)


```{r gender_model_pplots}
samp2 <- extract.samples(fit2)

# the plot will look nicer if we estimate both densities
# on the same range, so get the limits for both distributions
xlim <- range(0, exp(samp2$alpha), exp(samp2$alpha + samp2$beta))

# get 'density' objects for the two distributions
# (increasing n for better resolution)
pdens_girls <- density(exp(samp2$alpha), from=xlim[1], to=xlim[2], n = 2000)
pdens_boys <- density(exp(samp2$alpha + samp2$beta), from=xlim[1], to=xlim[2], n = 2000)

# figure out proper height to fit both densities on the same plot
ylim <- range(pdens_girls$y,pdens_boys$y,0)

# plot a blank canvas with those axes
plot(NA,xlim=xlim,ylim=ylim,xlab='rate',ylab='density')
# add each density
lines(pdens_girls,lty=1,lwd=4)
lines(pdens_boys, lty=3,lwd=4)
# and a legend
legend('topright',legend=c('boys','girls'),lty=c(3,1))
```


```{r gender_model_pplot2}
# plot them separately to see more detail
plot(pdens_girls)
plot(pdens_boys)
```

All told, it looks like the early 1990s stereotype of video games being a boys' passtime may have some support.

## Adding grade

In the early 1990s, video games were still widely considered a kids' passtime. Let's investigate that idea by seeing if students 'age out' of video games as they progress from 7th to 12th grade. We will use students' grade level, centered at 10th grade, as a further predictor.


```{r create_centered_grade}
# center
d$grade_c10 <- d$grade - 10
# plot
plot(table(d$grade_c10))
```


```{r gender_grade_model, cache=TRUE}
# build the model
m3 <- alist(
    hours_games ~ dpois(lambda),
    log(lambda) <- alpha + beta_m*male + beta_g*grade_c10,
    alpha ~ dnorm(4,1),
    beta_m ~ dnorm(0,0.5),
    beta_g ~ dnorm(0,0.3) # grade has a bigger range
)
# slim data
d_slim <- d[,c('hours_games','male','grade_c10')]
d_slim <- d_slim[complete.cases(d_slim),]

# fit
fit3 <- ulam(m3,data=d_slim,chains=4,cores=4,iter=4000)
```


```{r gender_grade_model_precis}
precis(fit3,prob=0.95)
exp(coef(fit3))
```

The values for $\alpha$ and $\beta_M$ did not change much, but the estimate for $\beta_G$ (the coefficient on grade level) is negative. Older students tend to report playing games for less time than younger students.

Specifically, our model predicts that for each grade level a student achieves, their rate of playing video games will be multiplied by about 0.86. This corresponds to an annual decrease of about 14%.



# Estimating with `glm` and `brm`

Soon, we will start shifting away from using the `rethinking` package for estimation in this class. While `rethinking` is great for learning about the models and forces us to think carefully through the assumptions of our models, it is neither practical nor efficient for everyday research. Below, we will estimate the same Poisson model of time spent playing video games (using gender and grade as predictors) using the `glm` function from base R and the `brm` function from the powerful `brms` package. Leaving off the prior, the model we are estimating is:

$$
\begin{aligned}
H_i &\sim Pois(\lambda_i)\\
\log(\lambda_i) &= \alpha + \beta_M M_i + \beta_G G_i
\end{aligned}
$$

## Keeping it basic: Poisson with `glm`

The `glm` function is built in to R (it's part of the `stats` package that is automatically loaded whenever R starts). It is used for maximum-likelihood (ML) estimation of generalized linear models such as logistic and poisson (see `?family` for a list of the distributions that it can handle). Estimating models with `glm` is fast and easy, and, because ML estimation considers regularizing parameters an estimation issue rather than a modeling issue, there are no priors to specify.

Notably, `glm` (as well as `brms` and `lme4`), expects your model to be specified in the standard "formula" format in R: `<outcome> ~ <predictor 1> + <precictor 2> + ...`

```{r glm_gender_grade}
# give `glm()` the model formula, the "family" of
# the regression model, and the data
fit3_glm <- glm(
    hours_games ~ male + grade_c10,
    family = poisson,
    data = d
)

# summary table
summary(fit3_glm)

# to get a confidence interval, use `confint()`
confint(fit3_glm, level= 0.95)

# exponentiate the coefficients
exp(coef(fit3_glm))
# you can even exponentiate the confidence intervals,
# BUT be aware that this is not formally meaningful
# for frequentist estimates so use only for rough information
exp(confint(fit3_glm, level= 0.95))
```


## More powerful and much slower: Poisson with `brm`

The `brm` function is the workhorse Bayesian estimator in the extremely powerful and flexible `brms` package. It can do a lot (including simple Poisson regressions). It uses the same `stan` backend as the `ulam` function in `rethinking`, but makes model specification more streamlined and does, in general, a much more efficient job with estimation. But this power and flexibility come at a price -- `brm` is _slow_. For each estimation, the default estimator in `brm` writes an entire program in C++, compiles it into a standalone app in the background, runs that app several times in the background, then collects and transforms those results for you, the hard-working researcher. Despite that, it is quite easy to use! For this simple model, we can use it exactly like `glm`.

```{r brm_gender_grade}
# load brms (normally we would do this at the top of the script)
library(brms)

# give `glm()` the model formula, the "family" of
# the regression model, and the data
# (the `cores` argument can be omitted, but it speeds things up)
fit3_brm <- brm(
    hours_games ~ male + grade_c10,
    family = poisson,
    data = d,
    cores=4 # run the estimators in parallel
)

# brms model summaries include a 95% CI by default
summary(fit3_brm)

# use `fixef` ("fixed effects") instead of 'coef'
# NOTE: the exponentiated "Est.Error" values should
# not be interpreted to mean anything
exp(fixef(fit3_brm))
```

Note that because we didn't specify priors for this model, `brm` picked priors for us. The output of `get_prior(fit3_brm)` will tell us what priors were used -- notably, it chose improper 'flat' priors for the coefficients on grade and gender! Letting `brm` choose priors is sometimes fine, especially for simple models like this or when just getting a feel for a model. But for more complex models this is a bad idea both conceptually and practically.

## Comparing estimators

We can compare the estimates produced by the three estimators we used here (`ulam`, `glm`, `brm`).

```{r compare_estimators}
print('ulam')
precis(fit3, prob = 0.95)[,c(1,3,4)]

print('glm')
cbind(coef(fit3_glm),confint(fit3_glm)) 

print('brm')
fixef(fit3_brm)[,c(1,3,4)]
```