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>, …
whole model without control
model_order_overall <- clm(rating ~ form_spinout + socialcare + region + inherited,
data = filter(locations_cleaned, domain == "Overall"),
link = "logit")
model_order_safe <- clm(rating ~ form_spinout + socialcare + region + inherited,
data = filter(locations_cleaned, domain == "Safe"),
link = "logit")
model_order_effective <- clm(rating ~ form_spinout + socialcare + region + inherited,
data = filter(locations_cleaned, domain == "Effective"),
link = "logit")
model_order_caring <- clm(rating ~ form_spinout + socialcare + region + inherited,
data = filter(locations_cleaned, domain == "Caring"),
link = "logit")
model_order_well_led <- clm(rating ~ form_spinout + socialcare + region + inherited,
data = filter(locations_cleaned, domain == "Well-led"),
link = "logit")
model_order_responsive <- clm(rating ~ form_spinout + socialcare + region + inherited,
data = filter(locations_cleaned, 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
| Inadequate |
Req improv |
-3.811 |
-3.457 |
-4.162 |
-6.936 |
-3.620 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| Req improv |
Good |
-0.957 |
-0.478 |
-1.092 |
-3.579 |
-1.015 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| Good |
Outstanding |
3.011 |
5.611 |
3.966 |
2.792 |
2.841 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| form_spinoutNSP_CIC |
0.781 |
1.011 |
0.236 |
0.185 |
0.675 |
0.771 |
|
(<0.001) *** |
(<0.001) *** |
(0.366) |
(0.521) |
(0.002) ** |
(<0.001) *** |
| form_spinoutSP_CIC |
1.063 |
1.488 |
0.806 |
0.537 |
0.469 |
0.897 |
|
(<0.001) *** |
(<0.001) *** |
(0.025) * |
(0.142) |
(0.107) |
(0.003) ** |
| socialcareTRUE |
0.459 |
1.363 |
0.501 |
-0.714 |
0.232 |
0.845 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(0.009) ** |
(<0.001) *** |
| inheritedTRUE |
0.117 |
0.552 |
-0.089 |
-0.201 |
0.224 |
0.096 |
|
(0.386) |
(0.004) ** |
(0.663) |
(0.356) |
(0.173) |
(0.568) |
| :———————- |
————-: |
————-: |
————-: |
————-: |
————-: |
————-: |
| Num.Obs. |
3766 |
3262 |
2920 |
3242 |
3261 |
3249 |
| AIC |
6015.1 |
4456.2 |
3485.3 |
2826.8 |
5402.8 |
4841.4 |
| BIC |
6108.6 |
4547.5 |
3575.0 |
2918.0 |
5494.2 |
4932.7 |
| RMSE |
2.38 |
2.14 |
2.25 |
2.40 |
2.37 |
2.35 |
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
| Inadequate |
Req improv |
0.022 |
0.032 |
0.016 |
0.001 |
0.027 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| Req improv |
Good |
0.384 |
0.620 |
0.336 |
0.028 |
0.362 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| Good |
Outstanding |
20.312 |
273.430 |
52.786 |
16.315 |
17.136 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| form_spinoutNSP_CIC |
2.184 |
2.749 |
1.266 |
1.204 |
1.964 |
2.163 |
|
(<0.001) *** |
(<0.001) *** |
(0.366) |
(0.521) |
(0.002) ** |
(<0.001) *** |
| form_spinoutSP_CIC |
2.896 |
4.429 |
2.238 |
1.710 |
1.598 |
2.451 |
|
(<0.001) *** |
(<0.001) *** |
(0.025) * |
(0.142) |
(0.107) |
(0.003) ** |
| socialcareTRUE |
1.582 |
3.907 |
1.650 |
0.490 |
1.261 |
2.327 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(0.009) ** |
(<0.001) *** |
| inheritedTRUE |
1.124 |
1.736 |
0.915 |
0.818 |
1.252 |
1.101 |
|
(0.386) |
(0.004) ** |
(0.663) |
(0.356) |
(0.173) |
(0.568) |
| :———————- |
————-: |
————-: |
————-: |
————-: |
————-: |
————-: |
| Num.Obs. |
3766 |
3262 |
2920 |
3242 |
3261 |
3249 |
| AIC |
6015.1 |
4456.2 |
3485.3 |
2826.8 |
5402.8 |
4841.4 |
| BIC |
6108.6 |
4547.5 |
3575.0 |
2918.0 |
5494.2 |
4932.7 |
| RMSE |
2.38 |
2.14 |
2.25 |
2.40 |
2.37 |
2.35 |
whole model with control
model_order_overall_covid <- clm(rating ~ form_spinout + after_covid +
socialcare + region + inherited,
data = filter(locations_cleaned, domain == "Overall"),
link = "logit")
model_order_safe_covid <- clm(rating ~ form_spinout + after_covid + socialcare + region + inherited,
data = filter(locations_cleaned, domain == "Safe"),
link = "logit")
model_order_effective_covid <- clm(rating ~ form_spinout + after_covid + socialcare + region + inherited,
data = filter(locations_cleaned, domain == "Effective"),
link = "logit")
model_order_caring_covid <- clm(rating ~ form_spinout + after_covid + socialcare + region + inherited,
data = filter(locations_cleaned, domain == "Caring"),
link = "logit")
model_order_well_led_covid <- clm(rating ~ form_spinout + after_covid + socialcare + region + inherited,
data = filter(locations_cleaned, domain == "Well-led"),
link = "logit")
model_order_responsive_covid <- clm(rating ~ form_spinout + after_covid + socialcare + region + inherited,
data = filter(locations_cleaned, 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
| Inadequate |
Req improv |
-4.355 |
-3.940 |
-4.422 |
-7.111 |
-4.028 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| Req improv |
Good |
-1.383 |
-0.830 |
-1.324 |
-3.750 |
-1.349 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| Good |
Outstanding |
2.768 |
5.434 |
3.813 |
2.680 |
2.635 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| form_spinoutNSP_CIC |
0.810 |
1.251 |
0.337 |
0.246 |
0.812 |
0.864 |
|
(<0.001) *** |
(<0.001) *** |
(0.193) |
(0.398) |
(<0.001) *** |
(<0.001) *** |
| form_spinoutSP_CIC |
0.979 |
1.372 |
0.735 |
0.476 |
0.405 |
0.840 |
|
(<0.001) *** |
(0.002) ** |
(0.042) * |
(0.194) |
(0.165) |
(0.005) ** |
| after_covidTRUE |
-1.400 |
-1.286 |
-0.835 |
-0.676 |
-1.153 |
-0.761 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| socialcareTRUE |
0.517 |
1.440 |
0.489 |
-0.714 |
0.242 |
0.845 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(0.007) ** |
(<0.001) *** |
| inheritedTRUE |
-0.033 |
0.447 |
-0.185 |
-0.255 |
0.132 |
0.027 |
|
(0.813) |
(0.023) * |
(0.372) |
(0.245) |
(0.431) |
(0.874) |
| :———————- |
————-: |
————-: |
————-: |
————-: |
————-: |
————-: |
| Num.Obs. |
3766 |
3262 |
2920 |
3242 |
3261 |
3249 |
| AIC |
5781.0 |
4275.6 |
3436.5 |
2807.9 |
5255.2 |
4786.0 |
| BIC |
5880.8 |
4373.0 |
3532.1 |
2905.2 |
5352.7 |
4883.3 |
| RMSE |
2.36 |
2.12 |
2.25 |
2.40 |
2.36 |
2.35 |
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
| Inadequate |
Req improv |
0.013 |
0.019 |
0.012 |
0.001 |
0.018 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| Req improv |
Good |
0.251 |
0.436 |
0.266 |
0.024 |
0.259 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| Good |
Outstanding |
15.924 |
229.112 |
45.272 |
14.582 |
13.949 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| form_spinoutNSP_CIC |
2.248 |
3.495 |
1.401 |
1.279 |
2.251 |
2.373 |
|
(<0.001) *** |
(<0.001) *** |
(0.193) |
(0.398) |
(<0.001) *** |
(<0.001) *** |
| form_spinoutSP_CIC |
2.663 |
3.943 |
2.085 |
1.610 |
1.500 |
2.317 |
|
(<0.001) *** |
(0.002) ** |
(0.042) * |
(0.194) |
(0.165) |
(0.005) ** |
| after_covidTRUE |
0.247 |
0.276 |
0.434 |
0.509 |
0.316 |
0.467 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| socialcareTRUE |
1.676 |
4.222 |
1.630 |
0.490 |
1.274 |
2.328 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(0.007) ** |
(<0.001) *** |
| inheritedTRUE |
0.968 |
1.563 |
0.831 |
0.775 |
1.141 |
1.027 |
|
(0.813) |
(0.023) * |
(0.372) |
(0.245) |
(0.431) |
(0.874) |
| :———————- |
————-: |
————-: |
————-: |
————-: |
————-: |
————-: |
| Num.Obs. |
3766 |
3262 |
2920 |
3242 |
3261 |
3249 |
| AIC |
5781.0 |
4275.6 |
3436.5 |
2807.9 |
5255.2 |
4786.0 |
| BIC |
5880.8 |
4373.0 |
3532.1 |
2905.2 |
5352.7 |
4883.3 |
| RMSE |
2.36 |
2.12 |
2.25 |
2.40 |
2.36 |
2.35 |