Setting up the R session and Loading packages

# environment setup to run ordered logit properly
options(contrasts = rep("contr.treatment", 2))

This chunk loads all the packages to use

#packages for ordered logit
library(ordinal) # package for ordinal logit regression
library(brant) # brant test for the parallel assumption for ordered logit
library(MASS) # models that work with the brant test

library(tidyverse) # package for data cleaning and plotting
library(readxl) # package for reading excel file
library(broom) # extracting model summary as data frame
library(modelsummary) # deriving model tables
library(scales) # label percent
library(lubridate) # working with dates
library(marginaleffects) #to calculate marginal effects

Import the cleaned data

locations_coded <- read_csv("D:/Archive/Socialcare UK/data/locations_coded.csv")

Prepare the data for ordinal regression

#select relevant columns, rename and relabel 
locations_cleaned <- locations_coded %>% 
  # recode legal form types to be more readable / easier to present
  mutate(inherited = ifelse(inherited == "Y", TRUE, FALSE),
         rating = recode(rating, 
                         "Insufficient evidence to rate" = "NA",
                         "Requires improvement" = "Req improv"),
         date = ymd(publication_date)) %>% 
  
  # assign order in the rating levels
  mutate(rating = ordered(rating, levels = c("Inadequate","Req improv", "Good", "Outstanding")),
         form_spinout = case_when(form == "GOV" ~ "GOV",
                           form == "CIC" & spin_out == "TRUE" ~ "SP_CIC",
                           form == "CIC" & spin_out == "FALSE" ~ "NSP_CIC"),
         socialcare = ifelse(type == "Social Care Org", TRUE, FALSE)) %>% 
  
  # creating a new dummy variable for facility category
  mutate(year = year(date),
         after_covid = ifelse(year >= 2020, TRUE, FALSE),
         before_covid = ifelse(year <= 2019, TRUE, FALSE)) %>%

  # converting the ordinal variable to numerical 
  mutate(rating_num = case_when(rating == "Inadequate" ~ 1,
                                rating == "Req improv" ~ 2,
                                rating == "Good" ~ 3,
                                rating == "Outstanding" ~ 4)) 
# show first several rows of the data set derived 
head(locations_cleaned)
## # A tibble: 6 × 39
##   `Location ID` `Location ODS Code` `Location Name` `Care Home?` type           
##   <chr>         <chr>               <chr>           <chr>        <chr>          
## 1 1-1017570228  VN638               Bluebell Centre N            Social Care Org
## 2 1-1017570228  VN638               Bluebell Centre N            Social Care Org
## 3 1-1017570228  VN638               Bluebell Centre N            Social Care Org
## 4 1-1017570228  VN638               Bluebell Centre N            Social Care Org
## 5 1-1017570228  VN638               Bluebell Centre N            Social Care Org
## 6 1-1017570228  VN638               Bluebell Centre N            Social Care Org
## # ℹ 34 more variables: primary_cat <chr>, `Location Street Address` <chr>,
## #   `Location Address Line 2` <chr>, `Location City` <chr>,
## #   `Location Post Code` <chr>, `Location Local Authority` <chr>, region <chr>,
## #   `Location NHS Region` <chr>, `Location ONSPD CCG Code` <chr>,
## #   `Location ONSPD CCG` <chr>, `Location Commissioning CCG Code` <chr>,
## #   `Location Commissioning CCG Name` <chr>,
## #   `Service / Population Group` <chr>, domain <chr>, rating <ord>, …

drop GOV providers in extra primary_cad

cic_cats <- locations_cleaned %>%
  filter(form == "CIC") %>%
  distinct(primary_cat) %>% 
  pull(primary_cat)
locations_short <- locations_cleaned %>% 
  filter(primary_cat %in% cic_cats)
nrow(locations_cleaned)
## [1] 19931
nrow(locations_short)
## [1] 6807

whole model without control

model_order_overall <- clm(rating ~ form_spinout  + socialcare + region + inherited,
                data = filter(locations_short, domain == "Overall"),
                link = "logit")


model_order_safe <- clm(rating ~ form_spinout  + socialcare + region + inherited,
                data = filter(locations_short, domain == "Safe"),
                link = "logit")
model_order_effective <- clm(rating ~ form_spinout  + socialcare + region + inherited,
                data = filter(locations_short, domain == "Effective"),
                link = "logit")
model_order_caring <- clm(rating ~ form_spinout  + socialcare + region + inherited,
                data = filter(locations_short, domain == "Caring"),
                link = "logit")
model_order_well_led <- clm(rating ~ form_spinout + socialcare + region + inherited,
                data = filter(locations_short, domain == "Well-led"),
                link = "logit")
model_order_responsive <- clm(rating ~ form_spinout + socialcare + region + inherited,
                data = filter(locations_short, domain == "Responsive"),
                link = "logit")
ordinal_models <-
  modelsummary(
    list(
      "overall" = model_order_overall,
      "safe" = model_order_safe,
      "effective" = model_order_effective,
      "caring" = model_order_caring,
      "well-led" = model_order_well_led,
      "responsive" = model_order_responsive
    ),
    coef_omit = "region",
    exponentiate = F,
    statistic = "({p.value}) {stars}")
ordinal_models
overall safe effective caring well-led responsive
Inadequate Req improv -6.647 -5.776 -5.351
(<0.001) *** (<0.001) *** (<0.001) ***
Req improv Good -2.306 -1.577 -2.410 -3.626 -1.896
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Good Outstanding 2.794 5.100 4.029 3.427 2.793
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
form_spinoutNSP_CIC 0.296 0.309 -0.243 0.642 0.486 0.362
(0.155) (0.365) (0.504) (0.064) + (0.082) + (0.264)
form_spinoutSP_CIC 0.583 0.997 0.744 0.999 0.161 0.479
(0.033) * (0.053) + (0.126) (0.018) * (0.653) (0.252)
socialcareTRUE -0.322 0.264 -0.261 -0.174 -0.480 -0.317
(0.042) * (0.285) (0.366) (0.541) (0.023) * (0.212)
inheritedTRUE -0.767 -0.412 -0.695 -0.377 -0.193 -0.633
(<0.001) *** (0.179) (0.046) * (0.374) (0.478) (0.080) +
:———————- ————-: ————-: ————-: ————-: ————-: ————-:
Num.Obs. 1527 1055 1051 1052 1057 1054
AIC 1714.0 898.6 781.1 701.7 1362.0 897.6
BIC 1794.0 973.0 850.5 771.1 1436.4 967.1
RMSE 2.33 2.13 1.18 1.29 2.32 1.33 
ordinal_models_exp <-
  modelsummary(
    list(
      "overall" = model_order_overall,
      "safe" = model_order_safe,
      "effective" = model_order_effective,
      "caring" = model_order_caring,
      "well-led" = model_order_well_led,
      "responsive" = model_order_responsive
    ),
    coef_omit = "region",
    exponentiate = T,
    statistic = "({p.value}) {stars}")
ordinal_models_exp
overall safe effective caring well-led responsive
Inadequate Req improv 0.001 0.003 0.005
(<0.001) *** (<0.001) *** (<0.001) ***
Req improv Good 0.100 0.207 0.090 0.027 0.150
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Good Outstanding 16.354 164.070 56.220 30.770 16.323
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
form_spinoutNSP_CIC 1.345 1.363 0.784 1.900 1.626 1.436
(0.155) (0.365) (0.504) (0.064) + (0.082) + (0.264)
form_spinoutSP_CIC 1.791 2.709 2.105 2.715 1.174 1.615
(0.033) * (0.053) + (0.126) (0.018) * (0.653) (0.252)
socialcareTRUE 0.725 1.302 0.770 0.841 0.619 0.729
(0.042) * (0.285) (0.366) (0.541) (0.023) * (0.212)
inheritedTRUE 0.464 0.662 0.499 0.686 0.825 0.531
(<0.001) *** (0.179) (0.046) * (0.374) (0.478) (0.080) +
:———————- ————-: ————-: ————-: ————-: ————-: ————-:
Num.Obs. 1527 1055 1051 1052 1057 1054
AIC 1714.0 898.6 781.1 701.7 1362.0 897.6
BIC 1794.0 973.0 850.5 771.1 1436.4 967.1
RMSE 2.33 2.13 1.18 1.29 2.32 1.33

whole model with control

model_order_overall_covid <- clm(rating ~ form_spinout  + after_covid +
                                  socialcare + region + inherited,
                data = filter(locations_short, domain == "Overall"),
                link = "logit")

model_order_safe_covid <- clm(rating ~ form_spinout  + after_covid + socialcare + region + inherited,
                data = filter(locations_short, domain == "Safe"),
                link = "logit")
model_order_effective_covid <- clm(rating ~ form_spinout  + after_covid + socialcare + region + inherited,
                data = filter(locations_short, domain == "Effective"),
                link = "logit")
model_order_caring_covid <- clm(rating ~ form_spinout  + after_covid + socialcare + region + inherited,
                data = filter(locations_short, domain == "Caring"),
                link = "logit")
model_order_well_led_covid <- clm(rating ~ form_spinout + after_covid + socialcare + region + inherited,
                data = filter(locations_short, domain == "Well-led"),
                link = "logit")
model_order_responsive_covid <- clm(rating ~ form_spinout + after_covid + socialcare + region + inherited,
                data = filter(locations_short, domain == "Responsive"),
                link = "logit")
ordinal_models_covid <-
  modelsummary(
    list(
      "overall" = model_order_overall_covid,
      "safe" = model_order_safe_covid,
      "effective" = model_order_effective_covid,
      "caring" = model_order_caring_covid,
      "well-led" = model_order_well_led_covid,
      "responsive" = model_order_responsive_covid
    ),
    coef_omit = "region",
    exponentiate = F,
    statistic = "({p.value}) {stars}"
  )
ordinal_models_covid
overall safe effective caring well-led responsive
Inadequate Req improv -7.314 -6.037 -5.644
(<0.001) *** (<0.001) *** (<0.001) ***
Req improv Good -2.861 -1.803 -2.546 -3.660 -2.132
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Good Outstanding 2.572 5.024 3.967 3.399 2.696
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
form_spinoutNSP_CIC 0.373 0.489 -0.126 0.667 0.632 0.417
(0.072) + (0.161) (0.731) (0.055) + (0.023) * (0.198)
form_spinoutSP_CIC 0.581 0.931 0.714 0.983 0.122 0.451
(0.034) * (0.075) + (0.139) (0.020) * (0.733) (0.280)
after_covidTRUE -1.715 -1.062 -0.764 -0.222 -1.109 -0.447
(<0.001) *** (<0.001) *** (0.004) ** (0.464) (<0.001) *** (0.093) +
socialcareTRUE -0.217 0.339 -0.231 -0.165 -0.423 -0.301
(0.179) (0.176) (0.424) (0.563) (0.045) * (0.236)
inheritedTRUE -1.032 -0.589 -0.804 -0.406 -0.317 -0.689
(<0.001) *** (0.061) + (0.023) * (0.340) (0.251) (0.058) +
:———————- ————-: ————-: ————-: ————-: ————-: ————-:
Num.Obs. 1527 1055 1051 1052 1057 1054
AIC 1630.7 876.9 775.2 703.1 1329.9 896.7
BIC 1716.0 956.3 849.5 777.5 1409.3 971.1
RMSE 2.33 2.12 1.18 1.29 2.31 1.33 
ordinal_models_covid_exp <-
  modelsummary(
    list(
      "overall" = model_order_overall_covid,
      "safe" = model_order_safe_covid,
      "effective" = model_order_effective_covid,
      "caring" = model_order_caring_covid,
      "well-led" = model_order_well_led_covid,
      "responsive" = model_order_responsive_covid
    ),
    coef_omit = "region",
    exponentiate = T,
    statistic = "({p.value}) {stars}"
  )
ordinal_models_covid_exp
overall safe effective caring well-led responsive
Inadequate Req improv 0.001 0.002 0.004
(<0.001) *** (<0.001) *** (<0.001) ***
Req improv Good 0.057 0.165 0.078 0.026 0.119
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Good Outstanding 13.099 151.977 52.848 29.930 14.826
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
form_spinoutNSP_CIC 1.452 1.631 0.882 1.948 1.881 1.517
(0.072) + (0.161) (0.731) (0.055) + (0.023) * (0.198)
form_spinoutSP_CIC 1.789 2.536 2.042 2.672 1.130 1.570
(0.034) * (0.075) + (0.139) (0.020) * (0.733) (0.280)
after_covidTRUE 0.180 0.346 0.466 0.801 0.330 0.640
(<0.001) *** (<0.001) *** (0.004) ** (0.464) (<0.001) *** (0.093) +
socialcareTRUE 0.805 1.404 0.794 0.848 0.655 0.740
(0.179) (0.176) (0.424) (0.563) (0.045) * (0.236)
inheritedTRUE 0.356 0.555 0.448 0.666 0.728 0.502
(<0.001) *** (0.061) + (0.023) * (0.340) (0.251) (0.058) +
:———————- ————-: ————-: ————-: ————-: ————-: ————-:
Num.Obs. 1527 1055 1051 1052 1057 1054
AIC 1630.7 876.9 775.2 703.1 1329.9 896.7
BIC 1716.0 956.3 849.5 777.5 1409.3 971.1
RMSE 2.33 2.12 1.18 1.29 2.31 1.33