Data Preparation

Libraries

library(tidyverse)
library(ggplot2)
library(gridExtra)
library(MASS)
library(ordinal)
library(pscl)
library(survey)
library(srvyr)
library(Rcpp)
library(tidycensus)
library(acs)
library(GGally)
library(regclass)
library(brant)
library(lemon)
library(caret)
library(yardstick)

census_api_key("a68070de348db9ce6811fd50529b751fdfae0f05")

Data

data <- read.csv("DataRegressionAnalysis.csv")
data <- data %>% 
  mutate(CityCategory = as.factor(Code)) %>% 
  mutate(Population = as.numeric(Population)) %>% 
  mutate(WhitePop = as.numeric(WhitePop)) %>% 
  mutate(MHI = as.numeric(MHI)) %>% 
  mutate(DemVote = as.numeric(DemVote)) %>% 
  mutate(MedianAge = as.numeric(MedianAge)) %>% 
  mutate(ClimateNetwork = as.numeric(ClimateNetwork))

Selecting Data for Ordinal Model

model3_data <- data %>% 
  select(CityCategory, YearDummy, ClimatePlan, FirstPlan, 
         PopDummy, MHIDummy, Poverty, Gini, NonWhitePop,
         MedianAge, Bachelor, GovernmentType, DemVote, 
         LegacyCity, CoastalCity, Engagement, ClimateNetwork)

Ordinal Model

Model 3 (3 levels)

POLR

m3 <- polr(CityCategory ~ ., 
            data = model3_data,
            Hess = TRUE,
            method = "logistic")
summary(m3)
## Call:
## polr(formula = CityCategory ~ ., data = model3_data, Hess = TRUE, 
##     method = "logistic")
## 
## Coefficients:
##                   Value Std. Error t value
## YearDummy       2.31936    0.80815  2.8700
## ClimatePlan    -0.49267    0.74517 -0.6611
## FirstPlan      -0.37156    0.77339 -0.4804
## PopDummy        1.20912    0.76374  1.5832
## MHIDummy        2.83900    1.24445  2.2813
## Poverty         0.45045    0.16275  2.7677
## Gini           -5.47227   16.42225 -0.3332
## NonWhitePop    -0.04810    0.03327 -1.4457
## MedianAge      -0.09308    0.15137 -0.6149
## Bachelor        0.05699    0.07359  0.7745
## GovernmentType -0.16554    0.85534 -0.1935
## DemVote        -0.01101    0.04649 -0.2367
## LegacyCity     -1.71420    1.09447 -1.5662
## CoastalCity     2.89155    1.00580  2.8749
## Engagement      3.83973    1.77586  2.1622
## ClimateNetwork  0.96184    0.87779  1.0957
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2  6.3597  6.8685     0.9259
## 2|3  9.2161  6.9292     1.3301
## 
## Residual Deviance: 76.27723 
## AIC: 112.2772
odds <- exp(coef(m3))
odds
##      YearDummy    ClimatePlan      FirstPlan       PopDummy       MHIDummy 
##   10.169194868    0.610993420    0.689660908    3.350542044   17.098645573 
##        Poverty           Gini    NonWhitePop      MedianAge       Bachelor 
##    1.569016582    0.004201666    0.953043128    0.911117073    1.058648582 
## GovernmentType        DemVote     LegacyCity    CoastalCity     Engagement 
##    0.847433832    0.989055315    0.180106941   18.021271787   46.512736334 
## ClimateNetwork 
##    2.616497691

CLM

m3b <- clm(CityCategory ~ ., 
            data = model3_data,
            na.action = na.omit,
            link = "logit")
summary(m3b)
## formula: 
## CityCategory ~ YearDummy + ClimatePlan + FirstPlan + PopDummy + MHIDummy + Poverty + Gini + NonWhitePop + MedianAge + Bachelor + GovernmentType + DemVote + LegacyCity + CoastalCity + Engagement + ClimateNetwork
## data:    model3_data
## 
##  link  threshold nobs logLik AIC    niter max.grad cond.H 
##  logit flexible  58   -38.14 112.28 6(0)  3.14e-13 3.2e+07
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)   
## YearDummy       2.31936    0.80814   2.870  0.00410 **
## ClimatePlan    -0.49267    0.74516  -0.661  0.50851   
## FirstPlan      -0.37154    0.77340  -0.480  0.63094   
## PopDummy        1.20911    0.76374   1.583  0.11339   
## MHIDummy        2.83900    1.24452   2.281  0.02254 * 
## Poverty         0.45045    0.16276   2.768  0.00565 **
## Gini           -5.47226   16.42430  -0.333  0.73900   
## NonWhitePop    -0.04809    0.03327  -1.446  0.14831   
## MedianAge      -0.09308    0.15133  -0.615  0.53851   
## Bachelor        0.05699    0.07359   0.774  0.43867   
## GovernmentType -0.16556    0.85535  -0.194  0.84653   
## DemVote        -0.01100    0.04648  -0.237  0.81285   
## LegacyCity     -1.71420    1.09454  -1.566  0.11732   
## CoastalCity     2.89154    1.00583   2.875  0.00404 **
## Engagement      3.83972    1.77580   2.162  0.03060 * 
## ClimateNetwork  0.96185    0.87781   1.096  0.27319   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Threshold coefficients:
##     Estimate Std. Error z value
## 1|2    6.360      6.867   0.926
## 2|3    9.216      6.928   1.330
pR2(m3b)
## fitting null model for pseudo-r2
##         llh     llhNull          G2    McFadden        r2ML        r2CU 
## -38.1386153 -63.6497120  51.0221934   0.4008046   0.5850897   0.6584243

Model 3 Selection (Both)

m3_both <- stepAIC(m3b, direction = "both", trace = FALSE)
summary(m3_both)
## formula: 
## CityCategory ~ YearDummy + PopDummy + MHIDummy + Poverty + NonWhitePop + LegacyCity + CoastalCity + Engagement
## data:    model3_data
## 
##  link  threshold nobs logLik AIC    niter max.grad cond.H 
##  logit flexible  58   -40.28 100.55 5(0)  3.94e-07 7.9e+05
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## YearDummy    2.39837    0.73169   3.278 0.001046 ** 
## PopDummy     1.45276    0.64462   2.254 0.024217 *  
## MHIDummy     3.22871    1.05846   3.050 0.002285 ** 
## Poverty      0.39969    0.11758   3.399 0.000675 ***
## NonWhitePop -0.06602    0.02453  -2.691 0.007119 ** 
## LegacyCity  -1.52192    0.92220  -1.650 0.098878 .  
## CoastalCity  2.74062    0.80630   3.399 0.000676 ***
## Engagement   3.69670    1.44539   2.558 0.010541 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Threshold coefficients:
##     Estimate Std. Error z value
## 1|2    9.143      2.755   3.319
## 2|3   11.816      2.969   3.980
pR2(m3_both)
## fitting null model for pseudo-r2
##         llh     llhNull          G2    McFadden        r2ML        r2CU 
## -40.2773517 -63.6497120  46.7447207   0.3672029   0.5533337   0.6226880
m3_both2 <- polr(CityCategory ~ YearDummy + PopDummy + MHIDummy + Poverty + NonWhitePop + LegacyCity + CoastalCity + Engagement, 
            data = model3_data,
            Hess = TRUE,
            method = "logistic")
summary(m3_both2)
## Call:
## polr(formula = CityCategory ~ YearDummy + PopDummy + MHIDummy + 
##     Poverty + NonWhitePop + LegacyCity + CoastalCity + Engagement, 
##     data = model3_data, Hess = TRUE, method = "logistic")
## 
## Coefficients:
##                Value Std. Error t value
## YearDummy    2.39837    0.73169   3.278
## PopDummy     1.45276    0.64462   2.254
## MHIDummy     3.22872    1.05846   3.050
## Poverty      0.39969    0.11758   3.399
## NonWhitePop -0.06602    0.02453  -2.691
## LegacyCity  -1.52194    0.92220  -1.650
## CoastalCity  2.74062    0.80627   3.399
## Engagement   3.69670    1.44538   2.558
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2  9.1426  2.7548     3.3188
## 2|3 11.8165  2.9689     3.9801
## 
## Residual Deviance: 80.5547 
## AIC: 100.5547
pR2(m3_both2)
## fitting null model for pseudo-r2
##         llh     llhNull          G2    McFadden        r2ML        r2CU 
## -40.2773517 -63.6497120  46.7447207   0.3672029   0.5533337   0.6226880
#Confirming with other formulas

m3_both2b <- polr(CityCategory ~ YearDummy + PopDummy + MHIDummy + Poverty + NonWhitePop + LegacyCity + CoastalCity + Engagement + DemVote, 
            data = model3_data,
            Hess = TRUE,
            method = "logistic")
summary(m3_both2b)
## Call:
## polr(formula = CityCategory ~ YearDummy + PopDummy + MHIDummy + 
##     Poverty + NonWhitePop + LegacyCity + CoastalCity + Engagement + 
##     DemVote, data = model3_data, Hess = TRUE, method = "logistic")
## 
## Coefficients:
##                Value Std. Error t value
## YearDummy    2.39197    0.73236  3.2661
## PopDummy     1.38578    0.66926  2.0706
## MHIDummy     3.15372    1.07672  2.9290
## Poverty      0.39739    0.11678  3.4030
## NonWhitePop -0.06789    0.02527 -2.6864
## LegacyCity  -1.67402    1.01665 -1.6466
## CoastalCity  2.72159    0.81106  3.3556
## Engagement   3.87043    1.55412  2.4904
## DemVote      0.01177    0.03313  0.3554
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2  9.7959  3.3507     2.9236
## 2|3 12.4752  3.5430     3.5211
## 
## Residual Deviance: 80.42835 
## AIC: 102.4284
pR2(m3_both2b)
## fitting null model for pseudo-r2
##         llh     llhNull          G2    McFadden        r2ML        r2CU 
## -40.2141768 -63.6497120  46.8710704   0.3681955   0.5543057   0.6237818

Analysis of Ordinal Model

exp(m3_both$coefficients)
##          1|2          2|3    YearDummy     PopDummy     MHIDummy      Poverty 
## 9.345413e+03 1.354649e+05 1.100526e+01 4.274907e+00 2.524718e+01 1.491362e+00 
##  NonWhitePop   LegacyCity  CoastalCity   Engagement 
## 9.361135e-01 2.182922e-01 1.549658e+01 4.031386e+01
round(confint(m3_both), 4)
##               2.5 %  97.5 %
## YearDummy    1.0547  3.9620
## PopDummy     0.2290  2.7848
## MHIDummy     1.2821  5.4796
## Poverty      0.1881  0.6550
## NonWhitePop -0.1180 -0.0207
## LegacyCity  -3.3971  0.2662
## CoastalCity  1.2570  4.4563
## Engagement   1.0965  6.7902
round(exp(coef(m3_both)), 4)
##         1|2         2|3   YearDummy    PopDummy    MHIDummy     Poverty 
##   9345.4134 135464.8602     11.0053      4.2749     25.2472      1.4914 
## NonWhitePop  LegacyCity CoastalCity  Engagement 
##      0.9361      0.2183     15.4966     40.3139
fit_m3_both <- lm(Code ~ YearDummy + PopDummy + MHIDummy + Poverty + NonWhitePop + CoastalCity + LegacyCity + Engagement, data = data)
VIF(fit_m3_both)
##   YearDummy    PopDummy    MHIDummy     Poverty NonWhitePop CoastalCity 
##    1.110251    1.154398    1.964765    3.083966    1.484721    1.252787 
##  LegacyCity  Engagement 
##    1.769131    1.210297
view_cov(m3, corr = T)
view_cov(m3_both, corr = T)
anova(m3b, m3_both)
brant(m3_both2)
## -------------------------------------------- 
## Test for X2  df  probability 
## -------------------------------------------- 
## Omnibus      3.57    8   0.89
## YearDummy    1.4 1   0.24
## PopDummy 0.19    1   0.66
## MHIDummy 0.41    1   0.52
## Poverty      0.67    1   0.41
## NonWhitePop  0.27    1   0.61
## LegacyCity   0.37    1   0.54
## CoastalCity  0.21    1   0.65
## Engagement   0   1   0.99
## -------------------------------------------- 
## 
## H0: Parallel Regression Assumption holds