---
title: 'Multilevel Logistic Regression: Respiritory Data'
author: "Carolyn J. Anderson"
date: "4/19/2022"
output: 
    html_document:
      toc: yes
---

The R in this document reproduces the results in the lecture on multilevel logistic regression.  The data are from Skrondal & Rabe-Hesketh (2004), which were also analyzed by Zeger \& Karim (1991), Diggle et al. (2002), but originally from Sommer et al. (1983).

# Setup

The package that we'll use

```{r setup, message=FALSE, warning=FALSE,include=TRUE}
library(lme4)
library(lmerTest)
library(gee)    # for generalized estimating equations
library(MASS)   # For confidence limits for parameters
```

And read in the data

```{r data}
resp <- read.table(file="D:/Dropbox/edps587/lectures/10 logreg/respitory.txt", header=TRUE)

head(resp)
```

# Logistic regression

This is the example shown on page 26, except that in notes I give Wald 

```{r logreg}

# glm = Generalized Linear model
model1 <- glm(resp ~ 1 + age, data=resp, family=binomial)
summary(model1)
```

Note that z^2 above equals Wald statistics in the lecture notes.  These are compared to N(0,1) to get the p-value.

The following puts 95% confidence levels on parameters (i.e., the coefficients or beta's)

```{r}
round(confint(model1, level = 0.95),digits=2)
```
For interpetation of the results

```{r} 
odds <- exp(coefficients(model1))
round(odds,digits=2)

odds2 <- 1/odds
round(odds2,digits=2) 
```

Out of curiosity and as a check on whether variance function of mean (i.e., binomial distribution OK).

```{r quasibinom}
model1b <- glm(resp ~ 1 + age, data=resp, family=quasibinomial)
summary(model1b)
```

Since the dispersion paramters is  is close to 1, binomial is probably just fine for these data.

# Generalized Estimating Equations

This yields a population averge model

```{r gee}
model.gee <- gee(resp ~1 + age, id, data=resp,family=binomial,corstr="exchangeable")
summary(model.gee)
```

# Random effects logistic regression via glmer

The command glmer is an option in the lme4 package.

### LaPlace

Now for the multilevel model or random effects model using the default estimation method (i.e., LaPlace)

```{r laplace}
mod1.laplace <- glmer(resp ~ 1 + age  + (1 | id), data=resp, family=binomial)
summary(mod1.laplace)
```

### MLE: quass quadrature
If you want MLE, glmer uses Guass quadrature --- MLE gold standard.  glmer will use adaptive guass quadrature and the number of quadrature points is indicated by nAGQ.  These results are reported ~ page 50,

```{r mle}
mod1.quad <- glmer(resp ~ 1 + age +  (1 | id),data=resp, family=binomial, nAGQ=10 )
summary(mod1.quad)
```

Computing the confidence intervals take a bit of time. The following gives 95% profile confidence intervals for all parameters

```{r}
round(confint(mod1.quad),digits=4)
```

For interpretations, we use the odds

```{r mleodds}
odds <- exp(fixef(mod1.quad))
round(odds,digits=2)
```

# Extract R script from this document

This is not run here, but you can enter it into your console with correct path and file name.
```{r extractR, eval=FALSE}
setwd("D:/Dropbox/edps587/lectures/10 logreg")
knitr::purl("Rmarkdown_respitory_example.Rmd","R_script_respitory_example.txt")
```