## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----eval = FALSE-------------------------------------------------------------
# devtools::install()

## ----setup, message = FALSE, warning = FALSE----------------------------------
library(tripaccess)
library(tidyverse)

## ----message = FALSE, warning = FALSE, out.width = "70%"----------------------
#> Filtered to households with at least one driver
house_with_drivers <- house |>
  filter(number_drivers > 0)

#> Summary statistics of vehicles by number of drivers
house_with_drivers |>
  group_by(number_drivers) |>
  summarize(
    households = n(),
    mean_vehicles = mean(number_vehicles),
    median_vehicles = median(number_vehicles),
    sd_vehicles = sd(number_vehicles)
  )

#> Filtered to households with at least one vehicle
house_with_vehicles <- house_with_drivers |>
  filter(number_vehicles > 0)

#> Plot household vehicles by number of drivers
ggplot(data = house_with_vehicles,
       aes(x = number_drivers,
           y = number_vehicles)) +
  geom_jitter(alpha = 0.08, width = 0.15, height = 0.15) +
  geom_smooth(method = lm, se = FALSE, formula = y ~ x, color = "blue") +
  labs(title = "Household Vehicles \nby Number of Drivers",
       x = "Number of Drivers in Household",
       y = "Number of Household Vehicles") +
  theme_bw()

#> Fit a simple linear regression model
vehicles_drivers_model <- lm(number_vehicles ~ number_drivers,
                             data = house_with_vehicles)
summary(vehicles_drivers_model)

#> Correlation between number of drivers and number of vehicles
cor(house_with_vehicles$number_drivers, house_with_vehicles$number_vehicles)

#> The slope estimates the average change in household vehicles for one
#> additional driver in the household.
coef(vehicles_drivers_model)["number_drivers"]

