---
title: "BMI510 Lab 6"
author: "McKay/Bozkurt"
date: "2025-04-02"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = F, warning = F, message = F)
library(tidyverse)
```

```{r}
# read data, do a little reorganization
d = read_csv("https://jlucasmckay.bmi.emory.edu/global/bmi510/Labs/bai-fog.csv") |> 
  mutate(sex = if_else(male == 1, "Male", "Female")) |> 
  select(-male) |> 
  select(-contains("trails")) |> 
  select(-bditotal) |> 
  filter(!is.na(baitotal))
```

## Demographics

```{r,results='asis'}
arsenal::tableby(foggroup~.,data = d |> select(-fogid)) |> summary(digits.p = 3, digits.pct = 0, digits = 0)
```


## Visualization

```{r}
d |>
  filter(foggroup != "NOFOG") |> 
  droplevels() |> 
  ggplot(aes(nfogqtotal,baitotal,color = pdduration,label = paste(" ",pdduration))) +
  geom_point() +
  geom_smooth(method = 'lm', se = F) + theme_minimal() + geom_text(hjust = 0, vjust = 1)
```

## Linear Model

```{r}
# Fit the linear model
m1 = d |> 
  filter(foggroup != "NOFOG") |> 
  droplevels() |> 
  lm(baitotal ~ nfogqtotal, data = _)

# View model summary
summary(m1)

```

```{r}
# Cuter
t1 = m1 |> broom::tidy()
t1 |> knitr::kable()
```


## Linear Model, Scaled Inputs and Outputs

```{r}
# Fit the linear model
m2 = d |> 
  filter(foggroup != "NOFOG") |> 
  droplevels() |> 
  lm(scale(baitotal) ~ scale(nfogqtotal), data = _)

# View model summary
summary(m2)

```

```{r}
# Cuter
t2 = m2 |> broom::tidy()
t2 |> knitr::kable()
```





## Manually adjust `baitotal` for `pdduration`

Notice what happens with the intercept here

```{r}
# adjust the baitotal variable for linear effects of pdduration
d$baitotal.adj = d |> lm(baitotal ~ pdduration, data = _) |> residuals()

d |>
  filter(foggroup != "NOFOG") |> 
  droplevels() |> 
  ggplot(aes(nfogqtotal,baitotal.adj,color = pdduration,label = paste(" ",pdduration))) +
  geom_point() +
  geom_smooth(method = 'lm', se = F) + theme_minimal() + geom_text(hjust = 0, vjust = 1)

```


## Linear Model, Manually Adjusted

```{r}
# Fit the linear model
m3 = d |> 
  filter(foggroup != "NOFOG") |> 
  droplevels() |> 
  lm(baitotal.adj ~ nfogqtotal, data = _)

# View model summary
summary(m3)
```

```{r}
# Cuter
t3 = m3 |> broom::tidy()
t3 |> knitr::kable()
```

## Use another term to do the adjustment in `lm`


```{r}
# Fit the linear model
m4 = d |> 
  filter(foggroup != "NOFOG") |> 
  droplevels() |> 
  lm(baitotal ~ nfogqtotal + pdduration, data = _)

# View model summary
summary(m4)
```

```{r}
# Cuter
t4 = m4 |> broom::tidy()
t4 |> knitr::kable()
```


```{r}
# Can even use "partial residuals" to illustrate this:

# Get partial residuals for nfogqtotal
library(car)
d.resid <- d |> 
  filter(foggroup != "NOFOG") |> 
  droplevels() |> 
  mutate(baitotal.resid = residuals(m4) + coef(m4)["nfogqtotal"] * nfogqtotal)


d.resid |>
  ggplot(aes(nfogqtotal,baitotal.resid,color = pdduration,label = paste(" ",pdduration))) +
  geom_point() +
  geom_smooth(method = 'lm', se = F) + theme_minimal() + geom_text(hjust = 0, vjust = 1)

```



## Compare models

```{r}

bind_rows(
  t1 |> mutate(model = "Simple") |> select(model,everything()),
  t2 |> mutate(model = "Scaled") |> select(model,everything()),
  t3 |> mutate(model = "Manually-Adjusted") |> select(model,everything()),
  t4 |> mutate(model = "Two Predictors") |> select(model,everything())) |> 
  filter(str_detect(term,"nfogq")) |> knitr::kable()

```