```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, 
                      fig.width = 8, fig.height = 6, 
                      fig.align = "center")

```

# Part 1: Linear regression in R

```{r}
library(tidyverse) #ggplot2, dplyr, etc.
library(reshape2) #need this for melt()
library(knitr) #need this for kable
library(MASS) #contains dataset
```

Load the birthwt data. This data contains 189 observations, 9 predictors, and an outcome, birthweight, available both as a continuous measure and a binary indicator for low birth weight.

```{r}
data(birthwt)
head(birthwt)
```

1. Plot a scatterplot of birthweight (bwt) and mother's weight (lwt).

**Solution**

```{r}
birthwt %>%
ggplot(aes(x = lwt, y = bwt)) + 
  geom_point() +
  labs(x = "Mother's weight (pounds)", y = "Birthweight (grams)")
```


2. Use OLS to fit the regression of birthweight on mother's weight.

**Solution**

```{r}
fit <- lm(bwt ~ lwt, data = birthwt)
summary(fit)
```


3. Extract the following: estimated coefficients, standard errors, variance-covariance matrix, and confidence intervals.


**Solution**

```{r}
# Estimated coefficients.
coefficients(fit)

# Standard errors. 
summary(fit)$coeff[, 2]

# Variance-covariance matrix. 
vcov(fit)

# Confidence intervals. 
confint(fit)
```


4. Plot the regression line and interpret the intercept and slope

**Solution**

```{r}
birthwt %>%
ggplot(aes(x = lwt, y = bwt)) +
  geom_point() +
  geom_smooth(method = "lm", col = "red", se = FALSE)
```

```{r}
# Or we can manually add the regression line with geom_abline()
birthwt %>%
ggplot(aes(x = lwt, y = bwt)) +
  geom_point() +
  geom_abline(slope = fit$coefficients[2], 
              intercept = fit$coefficients[1], col = "red")

```

5. Does the interpretation of the intercept make sense? How might we change this?

**Solution**

No, the intercept does not make sense because it implies that if a mother has zero pounds of weight, the expected birthweight is 2369 grams, which is not realistic.

We can center the mother's weight variable to make the intercept more interpretable. This means we will subtract the mean of mother's weight from each observation of mother's weight.

```{r}
birthwt <- birthwt %>% mutate(lwt_star = lwt - mean(lwt))
fit.new <- lm(bwt ~ lwt_star, data = birthwt)
summary(fit.new)
```


6. Now, we want to fit a model that includes race, mother's age, and smoking status in the model. Race takes on value 1 for white, 2 for black, and 3 for other. Mother's age is continuous. Smoking status is binary. Use OLS to calculate the coefficient estimates in this model.

**Solution**

```{r}
# wrrong
summary(lm(bwt ~ race + age + smoke, data = birthwt))
```

```{r}
fit2 <- lm(bwt ~ as.factor(race) + age + smoke, data = birthwt)
summary(fit2)
```

7. Interpret all the coefficient estimates.

**Solution**

- **`(Intercept)`: 3281.673 grams**
    - This is the estimated average birth weight in grams for a baseline individual. The "baseline" here refers to a mother who is White (`race = 1`, as it's the reference level for the `as.factor(race)` variable), does not smoke (`smoke = 0`), and has an age of 0.
    - **Important Note:** An age of 0 is not practically meaningful, so the intercept is often best understood as the starting point for the model's predictions, from which the effects of other variables are then added or subtracted.

- **`as.factor(race)2`: -444.069 grams**
    - This coefficient represents the estimated average difference in birth weight for infants born to **Black mothers (`race = 2`)** compared to **White mothers (`race = 1`)**, *holding all other variables (age and smoking status) constant*.
    - Specifically, on average, infants born to Black mothers are estimated to weigh approximately 444.069 grams less than infants born to White mothers, given the same age and smoking status.
    - The p-value of 0.004973 (which is less than 0.05) indicates that this difference is statistically significant.

- **`as.factor(race)3`: -447.858 grams**
    - This coefficient represents the estimated average difference in birth weight for infants born to **mothers of other races (`race = 3`)** compared to **White mothers (`race = 1`)**, *holding all other variables (age and smoking status) constant*.
    - On average, infants born to mothers of other races are estimated to weigh approximately 447.858 grams less than infants born to White mothers, given the same age and smoking status.
    - The p-value of 0.000226 (which is much less than 0.05) indicates that this difference is highly statistically significant.

- **`age`: 2.134 grams**
    - This coefficient represents the estimated average change in birth weight for every **one-year increase in the mother's age**, *holding all other variables (race and smoking status) constant*.
    - Specifically, for each additional year in the mother's age, the infant's birth weight is estimated to increase by approximately 2.134 grams.
    - The p-value of 0.827326 (which is much greater than 0.05) indicates that this effect is **not statistically significant**. This suggests that, in this model, the mother's age does not have a significant linear relationship with birth weight when considering race and smoking status.

- **`smoke`: -426.093 grams**
    - This coefficient represents the estimated average difference in birth weight for infants born to **mothers who smoke (`smoke = 1`)** compared to **mothers who do not smoke (`smoke = 0`)**, *holding all other variables (race and age) constant*.
    - On average, infants born to mothers who smoke are estimated to weigh approximately 426.093 grams less than infants born to mothers who do not smoke, given the same race and age.
    - The p-value of 0.000149 (which is much less than 0.05) indicates that this difference is highly statistically significant.


9. Print the results in Rmarkdown using kable().

```{r}
table <- data.frame(summary(fit2)$coef)
row.names(table) <- c("Intercept","White","Black","Mother's age", "Smoker")

knitr::kable(table,digits=3,align=rep('c', 2),
  col.names = c("estimate","standard error","test statistic","p-value"))
```

# Part 2: Generalized linear model 

Suppose that 2500 pregnant women are enrolled in a study and the outcome is the occurrence of preterm birth. Possible predictors of preterm birth include age of the woman, smoking, socioeconomic status, body mass index, bleeding during pregnancy, serum level of dde, and several dietary factors. 

1. Formulate the problem of selecting the important predictors of preterm birth in a generalized linear model (GLM) framework. 
2. Show the components of the GLM, including the link function and distribution (in exponential family form). 
3. Describe (briefly) how estimation and inference could proceed via a frequentist approach.

**Solution**:
$y_{i}=1$ if woman $i$ has preterm birth and $y_{i}=0$ otherwise $(i=1, \ldots, n)$
$y_{i} \sim \operatorname{Bernoulli}\left(\pi_{i}\right)$
Probability density function:
$$
\begin{aligned}
f\left(y_{i} ; \pi_{i}\right) &=\pi_{i}^{y_{i}}\left(1-\pi_{i}\right)^{1-y_{i}} \\
&=\exp \left\{y_{i} \log \pi_{i}+\left(1-y_{i}\right) \log \left(1-\pi_{i}\right)\right\} \\
&=\exp \left\{y_{i} \log \left(\frac{\pi_{i}}{1-\pi_{i}}\right)+\log \left(1-\pi_{i}\right)\right\} \\
&=\exp \left\{\frac{y_{i} \theta_{i}-b\left(\theta_{i}\right)}{a(\phi)}+c\left(y_{i}, \phi\right)\right\}
\end{aligned}
$$
where
$$
\theta_{i}=\log \left(\frac{\pi_{i}}{1-\pi_{i}}\right), \quad b\left(\theta_{i}\right)=\log \left(1+e^{\theta_{i}}\right), \quad a(\phi)=\phi=1,
$$
and $c\left(y_{i}, \phi\right)=0$.

Link function:
Any mapping from $\Re \rightarrow[0,1]$. A convenient choice is the canonical link,
$$
\eta_{i}=\theta_{i}=\log \left(\frac{\pi_{i}}{1-\pi_{i}}\right),
$$
which is the logit. The probit and complementary $\log -\log$ are alternatives.
Frequentist Estimation:
Maximum likelihood estimates can be obtained for a given model, say
$$
\log \left(\frac{\pi_{i}}{1-\pi_{i}}\right)=\mathbf{x}_{i}^{\prime} \boldsymbol{\beta}
$$

(where $\mathbf{x}_{i}$ is a $p \times 1$ vector of predictors) by iterative weighted least squares

Frequentist Inference:
One can select the important predictors to be included in the model by stepwise selection, using the AIC or BIC criterion.

Alternatively, one can just fit the model with all the predictors and then do inferences based on the MLEs and asymptotic standard errors. For example, for continuous predictors included as linear terms in the model, we can do a Wald test. Alternatively, we could do analysis of deviance (see notes for details) to test for significant differences in fit between the nested models with and without a particular predictor.


# Part 3: GLMs in R (Logistic regression)

Consider the space shuttle data in the MASS library. Consider modeling the use of the autolander as the outcome (variable name use). 

1. Fit a logistic regression model with autolander (variable auto) use (labeled as "auto" 1) versus not (0) as predicted by wind sign (variable wind). 

2. Give the estimated odds ratio for autolander use comparing head winds, labeled as "head" in the variable headwind (numerator) to tail winds (denominator).

3. Give the estimated odds ratio for autolander use comparing head winds (numerator) to tail winds (denominator) adjusting for wind strength from the variable magn.

4. If you fit a logistic regression model to a binary variable, for example use of the autolander, then fit a logistic regression model for one minus the outcome (not using the autolander) what happens to the coefficients?

```{r}
library(MASS)
?shuttle
data(shuttle)
head(shuttle)
```

**Solution**

1-2. 

```{r}
# not wrong, but needs to be careful with the interpretation of the coefficients
fit1 <- glm(use ~ wind, family = binomial(link = logit), data = shuttle)
summary(fit1)$coefficients
```

```{r}
shuttle$use.binary <- as.integer(shuttle$use == "auto")
fit2 <- glm(use.binary ~ wind, family = binomial(link = logit), data = shuttle)
summary(fit2)$coefficients
```

```{r}
# exp(intercept + 0*coef(windtail))/exp(intercept + coef(windtail))
exp(summary(fit2)$coefficients[1, 1]) / exp(sum(summary(fit2)$coefficients[, 1]))
1 / exp(summary(fit2)$coefficients[2, 1])
```

3. 
```{r}
fit3 <- glm(use.binary ~ wind + magn, family = binomial(link = logit), data = shuttle)
summary(fit3)$coefficients

exp(summary(fit3)$coefficients[1, 1] +
  summary(fit3)$coefficients[5, 1]) /
  exp(sum(summary(fit3)$coefficients[1:2, 1]) + summary(fit3)$coefficients[5, 1])
```

4. 
```{r}
fit4 <- glm(1 - use.binary ~ wind, family = binomial(link = logit), data = shuttle)
summary(fit4)$coefficients

summary(fit2)$coefficients
```

```{r}
fit5 <- lm(use.binary ~ wind, data = shuttle)
summary(fit5)$coefficients
```

# Part 4: GLMs in R (Poisson regression)

Consider the insect spray data InsectSprays. Fit a Poisson model using spray as a factor level. 

1. Report the estimated relative rate comapring spray A (numerator) to spray B (denominator).

2. Consider a Poisson glm with an offset, t. So, for example, a model of the form glm(count $\sim x+$ offset(t), family = poisson) where $x$ is a factor variable comparing a treatment (1) to a control (0) and $t$ is the natural log of a monitoring time. What is impact of the coefficient for $x$ if we fit the model glm(count $\sim x+$ offset(t2), family $=$ poisson) where $t 2<-\log (10)+t$ ? In other words, what happens to the coefficients if we change the units of the offset variable. (Note, adding log(10) on the log scale is multiplying by 10 on the original scale.)

```{r}
data("InsectSprays")
head(InsectSprays)
```

**Solution**:

```{r}
# Load data and fit Poisson model
model <- glm(count ~ spray, family = poisson, data = InsectSprays)

# Extract coefficients
coefs <- coef(model)
# sprayB coefficient = log(rate_B / rate_A)
rate_ratio_B_to_A <- exp(coefs["sprayB"])
relative_rate_A_to_B <- 1 / rate_ratio_B_to_A

# Result
relative_rate_A_to_B
```


```{r}
# Fit Poisson model with offset
t <- rep(log(10), nrow(InsectSprays)) # Example offset
model_offset <- glm(count ~ spray, offset = t, family = poisson, data = InsectSprays)

# Fit Poisson model with modified offset
t2 <- log(10) + t # Modified offset
model_offset2 <- glm(count ~ spray, offset = t2, family = poisson, data = InsectSprays)

# Compare coefficients
coefs_offset <- coef(model_offset)
coefs_offset2 <- coef(model_offset2)
# Coefficients for all
coefs_offset
coefs_offset2
```

```{r}
# The coefficients for spray A
exp(coefs_offset["(Intercept)"] + t[1])
exp(coefs_offset2["(Intercept)"] + t2[1])
```

