Assignment 6

# Set up environment
rm(list = ls())
gc()

# Load necessary libraries
library(clarify)
library(tidyverse)
library(ggplot2)
library(dplyr)
library(lubridate)
library(AER)
library(janitor)
library(pscl)

# Import the dataset
DHS_report <- read.csv("C:/Users/Shamp/OneDrive/Desktop/Data 712/Dataset for HW6/DHS_Daily_Report_20250320.csv")

## Data Cleaning

# Clean column names
DHS_report <- DHS_report %>%
  clean_names()

# Check for missing values
sum(is.na(DHS_report))

# Convert 'date_of_census' to Date format
DHS_report$date_of_census <- as.Date(DHS_report$date_of_census, format = "%m/%d/%Y")

# Verify the column exists and contains data
str(DHS_report)

head(DHS_report$Adults_in_Families_with_Children_in_Shelter)

# Derive predictors
DHS_report <- DHS_report %>%
  mutate(
    day_of_week = weekdays(date_of_census),  # Extract day of the week
    season = case_when(
      month(date_of_census) %in% c(12, 1, 2) ~ "Winter",
      month(date_of_census) %in% c(3, 4, 5) ~ "Spring",
      month(date_of_census) %in% c(6, 7, 8) ~ "Summer",
      month(date_of_census) %in% c(9, 10, 11) ~ "Fall"
    ),
    proportion_adults_in_families = adults_in_families_with_children_in_shelter / total_adults_in_shelter,
    family_size = total_individuals_in_families_with_children_in_shelter / families_with_children_in_shelter
  )

# View the updated dataset
head(DHS_report)

# Summary Statistics
# Calculate summary statistics for the dependent variable
summary(DHS_report$Adults_in_Families_with_Children_in_Shelter)

# Calculate the mean and variance
mean_count <- mean(DHS_report$adults_in_families_with_children_in_shelter, na.rm = TRUE)
var_count <- var(DHS_report$adults_in_families_with_children_in_shelter, na.rm = TRUE)

# Print the results
cat("Mean:", mean_count, "\nVariance:", var_count, "\n")

# Interpretation: The mean and variance of the dependent variable are calculated to check for overdispersion. The result shows that the dependent variable, adults_in_families_with_children_in_shelter, has a mean of 20,811 and a variance of 62,574,034. Since the variance is significantly greater than the mean, the data exhibits overdispersion. This means that the variability in the number of adults in families with children in shelters is much higher than what would be expected under a Poisson distribution. For this case a Negative Binomial regression model is more appropriate.

# Relationship with Season
ggplot(DHS_report, aes(x = season, y = adults_in_families_with_children_in_shelter)) +
  geom_boxplot(fill = "blue") +
  labs(title = "Shelter Occupancy by Season",
       x = "Season",
       y = "Number of Adults")

# This boxplot compares the number of adults in families with children in shelters across different seasons. It helps us identify whether shelter usage varies significantly by season, which could indicate seasonal trends such as higher usage in winter due to colder weather.

# Relationship Between Family Size and Shelter Occupancy
ggplot(DHS_report, aes(x = family_size, y = adults_in_families_with_children_in_shelter)) +
  geom_point(color = "blue", alpha = 0.6) +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  labs(title = "Relationship Between Family Size and Shelter Occupancy",
       x = "Family Size",
       y = "Number of Adults in Families with Children")

# Interpretation: The scatterplot shows the relationship between family size and the number of adults in families with children in shelters. Each point represents a data observation, and the linear regression line (Red) indicates the general trend. As family size increases, shelter occupancy also tends to increase, suggesting that larger families are more likely to seek shelter.

# Scatterplot for Proportion of Adults in Families
ggplot(DHS_report, aes(x = proportion_adults_in_families, y = adults_in_families_with_children_in_shelter)) +
  geom_point(color = "green", alpha = 0.6) +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  labs(title = "Relationship Between Proportion of Adults and Shelter Occupancy",
       x = "Proportion of Adults in Families",
       y = "Number of Adults in Families with Children")

# Interpretation: The scatterplot illustrates the relationship between the proportion of adults in families and the number of adults in families with children in shelters. Each point represents a different observation, showing how the proportion of adults in a family affects the number of adults in shelters. The linear regression line (Red) indicates a strong positive relationship, as the proportion of adults increases, shelter occupancy rises sharply.

# Regression Modeling

# Load necessary libraries
library(MASS)  # For Negative Binomial regression
library(car)   # For VIF calculation

# Fit Negative Binomial Regression Model
nb_model <- glm.nb(adults_in_families_with_children_in_shelter ~ family_size + proportion_adults_in_families,
                   data = DHS_report)

# Summarize the model
summary(nb_model)

# Check for Multicollinearity using Variance Inflation Factor (VIF)
vif(nb_model)

# Compare with Poisson Model (for reference)
poisson_model <- glm(adults_in_families_with_children_in_shelter ~ family_size + proportion_adults_in_families,
                     family = poisson, data = DHS_report)
summary(poisson_model)

# Overdispersion Test for Poisson Model
library(AER)
dispersiontest(poisson_model)

# Interpretation: The model coefficients reveal significant insights into shelter occupancy. The intercept (1.98219 on the log scale) indicates that when family size and the proportion of adults are zero, the expected number of adults in shelters is exp(1.98219) ≈ 7.26, serving as a baseline reference. Each additional family member multiplies the expected number of adults in shelters by exp(2.06972) ≈ 7.92, highlighting the greater housing challenges faced by larger families. Similarly, a one-unit increase in the proportion of adults multiplies the expected number of adults in shelters by exp(2.88417) ≈ 17.88, suggesting that families with fewer children and more adults experience higher housing instability. The model fit statistics indicate a good model fit. The residual deviance (1476.6) is significantly lower than the null deviance (131150.7), suggesting that the predictors improve the model. Additionally, the AIC value (24063) is much lower than that of the Poisson model (77036), confirming that the Negative Binomial model is a better fit for the data. The theta value (559.3174) confirms that the data exhibits overdispersion, justifying the use of Negative Binomial regression instead of Poisson regression. While the theta value is large, the significant reduction in deviance and AIC compared to the Poisson model supports the use of the Negative Binomial model. Key findings from the model show that larger families (p < 2e-16) and a higher proportion of adults in families (p < 2e-16) are significant risk factors for homelessness among families with children. These results suggest that housing policies should prioritize support for larger families and those with a higher proportion of adults.

# Simulate coefficients for the Negative Binomial model
sim_results <- sim(nb_model, n = 1000)

# Define new data for prediction
new_data <- data.frame(
  family_size = c(2, 3, 4, 5),  # Example family sizes
  proportion_adults_in_families = mean(DHS_report$proportion_adults_in_families)  # Hold proportion constant
)

# Calculate predicted counts
predicted_counts <- sim_apply(sim_results, FUN = predict, newdata = new_data, type = "response")

# Summarize the predicted counts
summary(predicted_counts)

# Interpretation of Predicted Counts: Using the clarify package, we simulated predictions from the Negative Binomial model to estimate shelter occupancy for different family sizes while keeping the proportion of adults in families constant. The results confirm that as family size increases, shelter occupancy grows exponentially, highlighting the heightened risk for larger families.

# Check the structure of predicted_counts to identify the correct column names
str(predicted_counts)

# Define new data for prediction
new_data <- data.frame(
  family_size = c(2, 3, 4, 5),  # Example family sizes
  proportion_adults_in_families = mean(DHS_report$proportion_adults_in_families)  # Hold proportion constant
)

# Calculate predicted counts
predicted_counts <- sim_apply(sim_results, FUN = predict, newdata = new_data, type = "response")

# Calculate mean and confidence intervals manually
predicted_counts_df <- data.frame(
  family_size = c(2, 3, 4, 5),
  mean = mean(predicted_counts),
  lower = apply(predicted_counts, 2, quantile, 0.025),
  upper = apply(predicted_counts, 2, quantile, 0.975)
)

# Create a bar plot with the extracted values
ggplot(predicted_counts_df, aes(x = factor(family_size), y = mean)) +
  geom_bar(stat = "identity", fill = "blue") +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
  labs(title = "Predicted Shelter Occupancy by Family Size",
       x = "Family Size",
       y = "Predicted Number of Adults") +
  theme_minimal()

# Summarize the predicted counts
summary(predicted_counts)

# The graph illustrates the predicted number of adults in shelters for families of sizes 2, 3, 4, and 5. The predicted number of adults generally increases with family size. However, the prediction for families of size 5 has a significantly wider confidence interval, indicating a higher degree of uncertainty in the model's prediction for this family size. This suggests that the model's predictions for larger families are less precise. Possible reasons for this uncertainty could be limitations in the model's fit or higher variability in the data for larger families.

# Visualize the distribution of the dependent variable
ggplot(DHS_report, aes(x = adults_in_families_with_children_in_shelter)) +
  geom_histogram(binwidth = 10, fill = "blue", color = "black") +
  labs(title = "Distribution of Adults in Families with Children in Shelter",
       x = "Number of Adults",
       y = "Frequency")

# This histogram shows the distribution of the number of adults in families with children in shelters. Most families have close to 10,000 adults, as indicated by the high peak. The data is skewed right, with a long tail extending to 30,000 adults. This suggests that while most families have a lower number of adults, there's significant variability, and some families have considerably higher numbers. The right skew indicates the mean number of adults is likely higher than the median. Further investigation is needed to understand the factors causing these higher counts.

# Plot 1: Adults in Families by Day of the Week
ggplot(DHS_report, aes(x = day_of_week, y = adults_in_families_with_children_in_shelter)) +
  geom_boxplot(fill = "orange") +
  labs(title = "Adults in Families by Day of the Week",
       x = "Day of the Week",
       y = "Number of Adults")

# The graph shows that the number of adults in families residing in shelters remains stable throughout the week. The boxplots for each day are similar, with median values hovering around the same level. The interquartile ranges (IQR) and whiskers also show comparable spread and overlap, indicating minimal impact of the day of the week on shelter occupancy. There are no significant outliers or extreme variations in the data.

# Plot 2: Adults in Families by Season
ggplot(DHS_report, aes(x = season, y = adults_in_families_with_children_in_shelter)) +
  geom_boxplot(fill = "darkcyan") +
  labs(title = "Adults in Families by Season",
       x = "Season",
       y = "Number of Adults")

# The boxplot illustrates that the number of adults in families with children in shelters varies significantly by season. Winter shows the highest median and interquartile range, indicating increased shelter occupancy during colder months. In contrast, summer has the lowest median and IQR, suggesting fewer adults in shelters. This seasonal trend aligns with expectations, as families may seek shelter more during winter due to harsh weather conditions.

# Fit Negative Binomial Regression Model
nb_model <- glm.nb(adults_in_families_with_children_in_shelter ~ day_of_week + season + proportion_adults_in_families + family_size,
                    data = DHS_report)
summary(nb_model)

# Use clarify to simulate and interpret results
sim_results <- sim(nb_model, n = 1000)
clarify_summary <- summary(sim_results)
clarify_summary

# Extract quantities of interest
new_data <- data.frame(
  day_of_week = "Monday",
  season = c("Winter", "Spring", "Summer", "Fall"),
  proportion_adults_in_families = mean(DHS_report$proportion_adults_in_families),
  family_size = mean(DHS_report$family_size)
)

# Corrected line: Use sim_apply() instead of predict()
predicted_counts <- sim_apply(sim_results, FUN = predict, newdata = new_data, type = "response")
predicted_counts

# Calculate mean and CI for plotting or further analysis
predicted_counts_df <- data.frame(
  season = c("Winter", "Spring", "Summer", "Fall"),
  mean = mean(predicted_counts),
  lower = apply(predicted_counts, 2, quantile, 0.025),
  upper = apply(predicted_counts, 2, quantile, 0.975)
)

predicted_counts_df

# Plot the predicted counts by season
library(ggplot2)
ggplot(predicted_counts_df, aes(x = season, y = mean)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2) +
  labs(title = "Predicted Shelter Occupancy by Season (Holding Monday Constant)",
       x = "Season",
       y = "Predicted Number of Adults") +
  theme_minimal()

# This bar graph displays the predicted number of adults in families with children in shelters across different seasons, specifically for Mondays. The error bars represent the confidence intervals of these predictions. Winter shows the highest predicted adult occupancy, while summer shows the lowest, with fall and spring being relatively close. This aligns with our earlier findings that shelter usage tends to be higher in winter, likely due to colder weather and increased need for shelter. The model's predictions support the hypothesis that seasonal factors significantly influence shelter occupancy.