model-01.qmd

---
title: "Thesis - Model"
toc: true
number-sections: true
code-fold: true
execute: 
  eval: true
  output: true
  error: true
format: gfm
---

# Load libraries

```{r}
library(tidyverse)
library(tidymodels)
library(modelr)
library(corrr)
library(ggfortify)
```

# Select Variables to work with
```{r}
mdl_df <- df %>% 
  select(
    c(
      name,
      muni_id,
      type,
      age_65_plus_pct,
      age_0_17_pct,
      min_wage_pct,
      high_educ_35_55_pct,
      ses_i_2015,
      peri_i_2015,
      likud_pct,
      coal_pct,
      budget_approved,
      sector, pop_log10, immig_1990_pct_log10, income_wage_log10, budget_approved_capita_log10
    )
  ) %>% 
  relocate(budget_approved, .after = last_col())
```


# First varaible to start with

Let's check which single variable is best to start with.
```{r}
df_mdl1a <- mdl_df %>% 
  select(-c(name, muni_id, budget_approved_capita_log10)) %>% 
  map(~ lm(mdl_df$budget_approved_capita_log10 ~ .x, data = mdl_df)) %>% 
  map(glance) %>% 
  map("r.squared") %>% 
  as_tibble() %>% 
  pivot_longer(everything() ,names_to = "var", values_to = "r.squared") %>% 
  arrange(desc(r.squared))

df_mdl1b <- mdl_df %>% 
  select(-c(name, muni_id, budget_approved)) %>% 
  map(~ lm(mdl_df$budget_approved ~ .x, data = mdl_df)) %>% 
  map(glance) %>% 
  map("r.squared") %>% 
  as_tibble() %>% 
  pivot_longer(everything() ,names_to = "var", values_to = "r.squared") %>% 
  arrange(desc(r.squared))

df_mdl1a
df_mdl1b

```

Looks like the population variable has the best correlation with the transformed approved SELA budget variable.
Let's examine how it fairs with linear regression assumptions.

```{r}
mdl_df %>% 
  ggplot(aes(pop_log10, budget_approved_capita_log10)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)
```
This seems like a plain result of the transformation itself, and not a result of a correlation between variables.
from now on we will examine the untransformed SELA budget variable.
Let's examine the same plot for the untransformed variable.
```{r}
mdl_df %>% 
  ggplot(aes(pop_log10, budget_approved)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)
```
This makes more sense. Also. it uncovers two more patterns we will try to identify: First, it seems like municipalities are grouped by some way, as noted in the dense dotted line-like groups. Second, it seems there are two separate linear trends.
Because the log-transformed variable is no longer necessary, we will drop it.
```{r}
mdl_df <- mdl_df %>% 
  select(-budget_approved_capita_log10)
```


Because we are now working with the original SELA budget variable, Let's examine the most correlated variables with the outcome variable.
as it seems, the 4 most correlated variables are socio-economic variables.
Let's start by examining their multicollinearity.
```{r}
mdl_df %>% 
  select(
    income_wage_log10, ses_i_2015, high_educ_35_55_pct, min_wage_pct
  ) %>% 
  correlate()
```
As it Seems, the four of them are highly correlated with one another. Since the SES variable is already an index variable taking into account multiple different variables, and since we do not have a lot of observations, we will continue to use only the SES variable.
Let's get rid of the unnecessary varaibles.
```{r}
mdl_df <- mdl_df %>% 
  select(-c(
    min_wage_pct,
    high_educ_35_55_pct,
    income_wage_log10
  ))
```

Let's visualize its correlation with SELA budget.
```{r}
mdl_df %>% 
  ggplot(aes(ses_i_2015, budget_approved)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE)

mdl1 <- lm(budget_approved ~ ses_i_2015, data = mdl_df)
mdl1 %>% glance()
```
This seems good enough, let's continue to the second variable.

# Second variable
Let's see which second variable is best for the model.
```{r}
df_mdl2 <- mdl_df %>% 
  select(-c(name, muni_id, ses_i_2015)) %>% 
  map(~ lm(budget_approved ~ .x + ses_i_2015, data = mdl_df)) %>% 
  map(glance) %>% 
  map("r.squared") %>% 
  as_tibble() %>% 
  pivot_longer(everything() ,names_to = "var", values_to = "r.squared") %>% 
  arrange(desc(r.squared))

df_mdl2
```
Seems like the municipality type is the best second variable.
Let's check for significance, visualize and examine the model.
```{r}
mdl2 <- lm(budget_approved ~ ses_i_2015 + type, data = mdl_df)

anv_21 <- anova(mdl1, mdl2)
anv_21
anv_21_p_value <- anv_21 %>% tidy() %>%  slice_tail() %>% pull(p.value)

rsq_1 <- mdl1 %>% rsquare(mdl_df)
rsq_2 <- mdl2 %>% rsquare(mdl_df)

df_predict <- data_grid(
  mdl_df, ses_i_2015, type, .model = mdl2
) %>% 
  add_predictions(mdl2)

mdl_df %>% 
  ggplot(aes(ses_i_2015, budget_approved)) +
  geom_point() +
  geom_line(data = df_predict, aes(y = pred), size = 1) +
  facet_wrap(~ type)

mdl2 %>% glance()
mdl2 %>% tidy()
```
The second model has an R-squared value of `r rsq_2` that represents an R-squared change of `r rsq_2 - rsq_1` with a confidence level p-value of `r anv_21_p_value`.

# Third variable

```{r}
df_mdl3 <- mdl_df %>% 
  select(-c(name, muni_id, ses_i_2015, type)) %>% 
  map(~ lm(budget_approved ~ .x + ses_i_2015 + type, data = mdl_df)) %>% 
  map(glance) %>% 
  map("r.squared") %>% 
  as_tibble() %>% 
  pivot_longer(everything() ,names_to = "var", values_to = "r.squared") %>% 
  arrange(desc(r.squared))

df_mdl3
```
The next best variable is peripheral index.
Let's check if it works.
```{r}
mdl3 <- lm(budget_approved ~ ses_i_2015 + type + peri_i_2015, data = mdl_df)

anv_32 <- anova(mdl2, mdl3)
anv_32
anv_32_p_value <- anv_32 %>% tidy() %>%  slice_tail() %>% pull(p.value)

rsq_3 <- mdl3 %>% rsquare(mdl_df)

df_predict <- data_grid(
  mdl_df, ses_i_2015, type, peri_i_2015, .model = mdl3
) %>% 
  add_predictions(mdl3)

mdl_df %>% 
  ggplot(aes(ses_i_2015, budget_approved)) +
  geom_point(aes(color = peri_i_2015)) +
  geom_smooth(method = "lm", se = FALSE) +
#  geom_line(data = df_predict, aes(y = pred), size = 1) +
  facet_wrap(~ type)

mdl3 %>% glance()
mdl3 %>% tidy()
```
The third model has an R-squared value of `r rsq_3` that represents an R-squared change of `r rsq_3 - rsq_2` with a confidence level p-value of `r anv_32_p_value`.

# Fourth variable
```{r}
df_mdl4 <- mdl_df %>% 
  select(-c(name, muni_id, ses_i_2015, type, peri_i_2015)) %>% 
  map(~ lm(budget_approved ~ .x + ses_i_2015 + type + peri_i_2015, data = mdl_df)) %>% 
  map(glance) %>% 
  map("r.squared") %>% 
  as_tibble() %>% 
  pivot_longer(everything() ,names_to = "var", values_to = "r.squared") %>% 
  arrange(desc(r.squared))

df_mdl4
```
The next best varaibale is population, which makes sense and obviously needed to be taken into account.
```{r}
mdl4 <- lm(budget_approved ~ ses_i_2015 + type + peri_i_2015 + pop_log10, data = mdl_df)

anv_43 <- anova(mdl3, mdl4)
anv_43
anv_43_p_value <- anv_43 %>% tidy() %>%  slice_tail() %>% pull(p.value)

rsq_4 <- mdl4 %>% rsquare(mdl_df)

mdl4 %>% glance()
mdl4 %>% tidy()
```
The fourth model has an R-squared value of `r rsq_4` that represents an R-squared change of `r rsq_4 - rsq_3` with a confidence level p-value of `r anv_43_p_value`.

# Fifth variable
```{r}
df_mdl5 <- mdl_df %>% 
  select(-c(name, muni_id, ses_i_2015, type, peri_i_2015, pop_log10)) %>% 
  map(~ lm(budget_approved ~ .x + ses_i_2015 + type + peri_i_2015 + pop_log10, data = mdl_df)) %>% 
  map(glance) %>% 
  map("r.squared") %>% 
  as_tibble() %>% 
  pivot_longer(everything() ,names_to = "var", values_to = "r.squared") %>% 
  arrange(desc(r.squared))

df_mdl5
```

It seems like this is as far as we can get.

# Trying to create a function for next variable diagnosis
We will create a function that gets a vector of variables that already went into the model, and finds the best next variable to insert to the model.

```{r}
calc_next_rsq <- function(data, vars = NULL){
  
  frmla <- as.formula(
    paste("budget_approved ~ .x", vars, sep = " + ")
  )
  
  data %>% 
    select(-c(name, muni_id, {{ vars }})) %>% 
    map(~ lm(frmla, data = data)) %>% 
    map(glance) %>% 
    map("r.squared") %>% 
    as_tibble() %>% 
    pivot_longer(everything() ,names_to = "var", values_to = "r.squared") %>% 
    arrange(desc(r.squared))
}
```