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
merged <- read_csv("D:/Archive/Socialcare UK/data/spinout_2017.csv")
Prepare the data for ordinal regression
#select relevant columns, rename and relabel
merged_cleaned <- merged %>%
# 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),
year2 = year-2013,
Year = factor(year)) %>%
# 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(merged_cleaned)
## # A tibble: 6 × 21
## project_id provider_name location_name type primary_cat form region domain
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 location1 Solihull Metro… Bluebell Cen… Soci… Community … GOV West … Safe
## 2 location2 Solihull Metro… Bluebell Cen… Soci… Community … GOV West … Effec…
## 3 location3 Solihull Metro… Bluebell Cen… Soci… Community … GOV West … Caring
## 4 location4 Solihull Metro… Bluebell Cen… Soci… Community … GOV West … Respo…
## 5 location5 Solihull Metro… Bluebell Cen… Soci… Community … GOV West … Well-…
## 6 location6 Solihull Metro… Bluebell Cen… Soci… Community … GOV West … Overa…
## # ℹ 13 more variables: rating <ord>, publication_date <dttm>,
## # report_type <chr>, std_name <chr>, level <chr>, spin_out <lgl>,
## # date <date>, form_spinout <chr>, socialcare <lgl>, year <dbl>, year2 <dbl>,
## # Year <fct>, rating_num <dbl>
Overall distribution
datasummary(form_spinout + socialcare + region + Year ~ 1, data = merged_short, fmt = 0)
| |
|
All |
| form_spinout |
GOV |
5963 |
|
NSP_CIC |
813 |
|
SP_CIC |
843 |
| socialcare |
FALSE |
2561 |
|
TRUE |
5058 |
| region |
East Midlands |
730 |
|
East of England |
851 |
|
London |
756 |
|
North East |
381 |
|
North West |
1304 |
|
South East |
1204 |
|
South West |
677 |
|
West Midlands |
690 |
|
Yorkshire and The Humber |
1026 |
| Year |
2013 |
0 |
|
2014 |
276 |
|
2015 |
786 |
|
2016 |
3459 |
|
2017 |
3098 |
Comparison table on Primary Inspection Category
datasummary(primary_cat ~ form, data = merged_cleaned, fmt = 0)
| primary_cat |
CIC |
GOV |
| Acute hospital - Independent specialist |
24 |
0 |
| Acute hospital - NHS non-specialist |
0 |
13135 |
| Acute hospital - NHS specialist |
0 |
1092 |
| Ambulance service |
0 |
361 |
| Community based adult social care services |
360 |
1944 |
| Community health - NHS & Independent |
300 |
954 |
| Dentists |
0 |
24 |
| GP Practices |
546 |
203 |
| Hospice services |
12 |
6 |
| Independent consulting doctors |
24 |
0 |
| Mental health - community & hospital - independent |
156 |
0 |
| Mental health - community & residential - NHS |
0 |
4242 |
| Out of hours |
6 |
6 |
| Prison Healthcare |
30 |
0 |
| Residential social care |
186 |
2832 |
| Urgent care services & mobile doctors |
12 |
18 |
datasummary(primary_cat ~ form_spinout, data = merged_cleaned, fmt = 0)
| primary_cat |
GOV |
NSP_CIC |
SP_CIC |
| Acute hospital - Independent specialist |
0 |
0 |
24 |
| Acute hospital - NHS non-specialist |
13135 |
0 |
0 |
| Acute hospital - NHS specialist |
1092 |
0 |
0 |
| Ambulance service |
361 |
0 |
0 |
| Community based adult social care services |
1944 |
252 |
108 |
| Community health - NHS & Independent |
954 |
102 |
198 |
| Dentists |
24 |
0 |
0 |
| GP Practices |
203 |
285 |
261 |
| Hospice services |
6 |
6 |
6 |
| Independent consulting doctors |
0 |
0 |
24 |
| Mental health - community & hospital - independent |
0 |
24 |
132 |
| Mental health - community & residential - NHS |
4242 |
0 |
0 |
| Out of hours |
6 |
6 |
0 |
| Prison Healthcare |
0 |
30 |
0 |
| Residential social care |
2832 |
102 |
84 |
| Urgent care services & mobile doctors |
18 |
6 |
6 |
whole models with continous year trend
model_order_overall <- clm(rating ~ form_spinout + socialcare + region + year2 + level,
data = filter(merged_short, domain == "Overall"),
link = "logit")
model_order_safe <- clm(rating ~ form_spinout + socialcare + region + year2 + level,
data = filter(merged_short, domain == "Safe"),
link = "logit")
model_order_effective <- clm(rating ~ form_spinout + socialcare + region + year2 + level,
data = filter(merged_short, domain == "Effective"),
link = "logit")
model_order_caring <- clm(rating ~ form_spinout + socialcare + region + year2 + level,
data = filter(merged_short, domain == "Caring"),
link = "logit")
model_order_well_led <- clm(rating ~ form_spinout + socialcare + region + year2 + level,
data = filter(merged_short, domain == "Well-led"),
link = "logit")
model_order_responsive <- clm(rating ~ form_spinout + socialcare + region + year2 + level,
data = filter(merged_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 |
-7.056 |
-5.504 |
-6.800 |
|
-6.312 |
-6.037 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
|
(<0.001) *** |
(<0.001) *** |
| Req improv|Good |
-4.270 |
-2.467 |
-3.189 |
-4.662 |
-3.416 |
-2.539 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| Good|Outstanding |
0.940 |
5.019 |
3.171 |
2.970 |
1.746 |
3.142 |
|
(0.012) * |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| form_spinoutNSP_CIC |
0.809 |
0.470 |
0.258 |
0.794 |
0.194 |
0.705 |
|
(<0.001) *** |
(0.116) |
(0.434) |
(0.044) * |
(0.465) |
(0.034) * |
| form_spinoutSP_CIC |
0.788 |
0.790 |
0.645 |
0.241 |
1.195 |
0.921 |
|
(<0.001) *** |
(0.010) * |
(0.058) + |
(0.559) |
(<0.001) *** |
(0.007) ** |
| socialcareTRUE |
-0.092 |
0.138 |
-0.768 |
-0.810 |
-0.504 |
0.170 |
|
(0.587) |
(0.575) |
(0.011) * |
(0.023) * |
(0.030) * |
(0.531) |
| year2 |
-0.778 |
-0.362 |
-0.215 |
0.061 |
-0.446 |
-0.096 |
|
(<0.001) *** |
(<0.001) *** |
(0.059) + |
(0.692) |
(<0.001) *** |
(0.396) |
| levelprovider |
-1.106 |
-1.359 |
-1.644 |
0.987 |
-1.339 |
-0.707 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(0.011) * |
(<0.001) *** |
(0.026) * |
| Num.Obs. |
1580 |
1209 |
1206 |
1206 |
1209 |
1209 |
| AIC |
1927.3 |
1295.8 |
1143.7 |
668.7 |
1571.4 |
1183.4 |
| BIC |
2013.1 |
1377.4 |
1225.3 |
745.2 |
1653.0 |
1265.0 |
| RMSE |
2.26 |
2.11 |
2.14 |
1.23 |
2.21 |
2.21 |
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.004 |
0.001 |
|
0.002 |
0.002 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
|
(<0.001) *** |
(<0.001) *** |
| Req improv|Good |
0.014 |
0.085 |
0.041 |
0.009 |
0.033 |
0.079 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| Good|Outstanding |
2.559 |
151.271 |
23.821 |
19.486 |
5.730 |
23.147 |
|
(0.012) * |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| form_spinoutNSP_CIC |
2.245 |
1.600 |
1.294 |
2.213 |
1.214 |
2.025 |
|
(<0.001) *** |
(0.116) |
(0.434) |
(0.044) * |
(0.465) |
(0.034) * |
| form_spinoutSP_CIC |
2.199 |
2.204 |
1.905 |
1.273 |
3.302 |
2.513 |
|
(<0.001) *** |
(0.010) * |
(0.058) + |
(0.559) |
(<0.001) *** |
(0.007) ** |
| socialcareTRUE |
0.912 |
1.148 |
0.464 |
0.445 |
0.604 |
1.185 |
|
(0.587) |
(0.575) |
(0.011) * |
(0.023) * |
(0.030) * |
(0.531) |
| year2 |
0.459 |
0.696 |
0.807 |
1.062 |
0.640 |
0.908 |
|
(<0.001) *** |
(<0.001) *** |
(0.059) + |
(0.692) |
(<0.001) *** |
(0.396) |
| levelprovider |
0.331 |
0.257 |
0.193 |
2.682 |
0.262 |
0.493 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(0.011) * |
(<0.001) *** |
(0.026) * |
| Num.Obs. |
1580 |
1209 |
1206 |
1206 |
1209 |
1209 |
| AIC |
1927.3 |
1295.8 |
1143.7 |
668.7 |
1571.4 |
1183.4 |
| BIC |
2013.1 |
1377.4 |
1225.3 |
745.2 |
1653.0 |
1265.0 |
| RMSE |
2.26 |
2.11 |
2.14 |
1.23 |
2.21 |
2.21 |
Models with Year fixed effects
model_order_overall_fe <- clm(rating ~ form_spinout + socialcare + region + Year + level,
data = filter(merged_short, domain == "Overall"),
link = "logit")
model_order_safe_fe <- clm(rating ~ form_spinout + socialcare + region + Year + level,
data = filter(merged_short, domain == "Safe"),
link = "logit")
model_order_effective_fe <- clm(rating ~ form_spinout + socialcare + region + Year + level,
data = filter(merged_short, domain == "Effective"),
link = "logit")
model_order_caring_fe <- clm(rating ~ form_spinout + socialcare + region + Year + level,
data = filter(merged_short, domain == "Caring"),
link = "logit")
model_order_well_led_fe <- clm(rating ~ form_spinout + socialcare + region + Year + level,
data = filter(merged_short, domain == "Well-led"),
link = "logit")
model_order_responsive_fe <- clm(rating ~ form_spinout + socialcare + region + Year + level,
data = filter(merged_short, domain == "Responsive"),
link = "logit")
ordinal_models_fe <-
modelsummary(
list(
"overall" = model_order_overall_fe,
"safe" = model_order_safe_fe,
"effective" = model_order_effective_fe,
"caring" = model_order_caring_fe,
"well-led" = model_order_well_led_fe,
"responsive" = model_order_responsive_fe
),
coef_omit = "region",
exponentiate = F,
statistic = "({p.value}) {stars}")
ordinal_models_fe
| |
overall |
safe |
effective |
caring |
well-led |
responsive |
| Inadequate|Req improv |
-5.746 |
-4.520 |
-6.258 |
|
-5.176 |
-5.669 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
|
(<0.001) *** |
(<0.001) *** |
| Req improv|Good |
-2.957 |
-1.484 |
-2.647 |
-4.115 |
-2.282 |
-2.172 |
|
(<0.001) *** |
(0.002) ** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| Good|Outstanding |
2.295 |
6.044 |
3.724 |
3.552 |
2.920 |
3.513 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| form_spinoutNSP_CIC |
0.839 |
0.506 |
0.267 |
0.831 |
0.202 |
0.713 |
|
(<0.001) *** |
(0.092) + |
(0.419) |
(0.035) * |
(0.447) |
(0.032) * |
| form_spinoutSP_CIC |
0.735 |
0.721 |
0.613 |
0.164 |
1.125 |
0.891 |
|
(<0.001) *** |
(0.020) * |
(0.074) + |
(0.694) |
(<0.001) *** |
(0.010) ** |
| socialcareTRUE |
-0.140 |
0.105 |
-0.783 |
-0.855 |
-0.525 |
0.154 |
|
(0.412) |
(0.672) |
(0.010) ** |
(0.017) * |
(0.025) * |
(0.570) |
| Year2015 |
-0.499 |
0.303 |
0.224 |
0.445 |
0.600 |
0.250 |
|
(0.269) |
(0.492) |
(0.653) |
(0.550) |
(0.184) |
(0.613) |
| Year2016 |
-0.786 |
0.054 |
-0.056 |
0.937 |
-0.168 |
0.120 |
|
(0.055) + |
(0.887) |
(0.899) |
(0.172) |
(0.673) |
(0.788) |
| Year2017 |
-1.918 |
-0.570 |
-0.348 |
0.615 |
-0.663 |
-0.046 |
|
(<0.001) *** |
(0.143) |
(0.437) |
(0.373) |
(0.101) |
(0.919) |
| levelprovider |
-1.061 |
-1.257 |
-1.588 |
1.104 |
-1.215 |
-0.656 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(0.005) ** |
(<0.001) *** |
(0.043) * |
| Num.Obs. |
1580 |
1209 |
1206 |
1206 |
1209 |
1209 |
| AIC |
1921.6 |
1294.2 |
1146.8 |
669.2 |
1570.0 |
1186.8 |
| BIC |
2018.1 |
1385.9 |
1238.5 |
755.9 |
1661.8 |
1278.5 |
| RMSE |
2.27 |
2.11 |
2.14 |
1.23 |
2.21 |
2.21 |
ordinal_models_fe_exp <-
modelsummary(
list(
"overall" = model_order_overall_fe,
"safe" = model_order_safe_fe,
"effective" = model_order_effective_fe,
"caring" = model_order_caring_fe,
"well-led" = model_order_well_led_fe,
"responsive" = model_order_responsive_fe
),
coef_omit = "region",
exponentiate = T,
statistic = "({p.value}) {stars}")
ordinal_models_fe_exp
| |
overall |
safe |
effective |
caring |
well-led |
responsive |
| Inadequate|Req improv |
0.003 |
0.011 |
0.002 |
|
0.006 |
0.003 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
|
(<0.001) *** |
(<0.001) *** |
| Req improv|Good |
0.052 |
0.227 |
0.071 |
0.016 |
0.102 |
0.114 |
|
(<0.001) *** |
(0.002) ** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| Good|Outstanding |
9.928 |
421.397 |
41.410 |
34.868 |
18.543 |
33.553 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
| form_spinoutNSP_CIC |
2.314 |
1.658 |
1.306 |
2.297 |
1.224 |
2.040 |
|
(<0.001) *** |
(0.092) + |
(0.419) |
(0.035) * |
(0.447) |
(0.032) * |
| form_spinoutSP_CIC |
2.085 |
2.057 |
1.846 |
1.178 |
3.079 |
2.439 |
|
(<0.001) *** |
(0.020) * |
(0.074) + |
(0.694) |
(<0.001) *** |
(0.010) ** |
| socialcareTRUE |
0.869 |
1.111 |
0.457 |
0.425 |
0.592 |
1.167 |
|
(0.412) |
(0.672) |
(0.010) ** |
(0.017) * |
(0.025) * |
(0.570) |
| Year2015 |
0.607 |
1.354 |
1.251 |
1.560 |
1.823 |
1.284 |
|
(0.269) |
(0.492) |
(0.653) |
(0.550) |
(0.184) |
(0.613) |
| Year2016 |
0.456 |
1.056 |
0.946 |
2.553 |
0.845 |
1.128 |
|
(0.055) + |
(0.887) |
(0.899) |
(0.172) |
(0.673) |
(0.788) |
| Year2017 |
0.147 |
0.566 |
0.706 |
1.850 |
0.515 |
0.955 |
|
(<0.001) *** |
(0.143) |
(0.437) |
(0.373) |
(0.101) |
(0.919) |
| levelprovider |
0.346 |
0.284 |
0.204 |
3.016 |
0.297 |
0.519 |
|
(<0.001) *** |
(<0.001) *** |
(<0.001) *** |
(0.005) ** |
(<0.001) *** |
(0.043) * |
| Num.Obs. |
1580 |
1209 |
1206 |
1206 |
1209 |
1209 |
| AIC |
1921.6 |
1294.2 |
1146.8 |
669.2 |
1570.0 |
1186.8 |
| BIC |
2018.1 |
1385.9 |
1238.5 |
755.9 |
1661.8 |
1278.5 |
| RMSE |
2.27 |
2.11 |
2.14 |
1.23 |
2.21 |
2.21 |