# 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(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
library(ordinal) # package for ordinal logit regression

Data Preparation

First, import the sampled and coded data set

#import the raw data file
cqc_skills <- read_csv("data/cleaned_sfc_ipw.csv") %>% 
  mutate(form = fct_relevel(form, "FPO"),
         rating = factor(rating),
         Year = factor(year),
         during_covid = after_covid,
         staff_level = staff_level/100)

Detailes about new variables

Employee turnover rate is the percent of employees who leave a company within a specific time period. Turnover rate is commonly calculated by month, quarter, or year and includes both voluntary and involuntary losses.

\[ \text{Employee turnover rate} = \frac{\text{Number of employees who left}}{\text{Average number of employees}} \times 100 \]

\[ \text{Average number of employees} = \frac{\text{(Headcount at the begining of the timeframe + headcount at the end of the timeframe)}}{2} \]

ordered logit models

Modeling on rating withough staff control

model_order_overall <- clm(rating ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Overall"),
                link = "logit",
                weights = ipw)

model_order_safe <- clm(rating ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Safe"),
                link = "logit",
                weights = ipw)

model_order_effective <- clm(rating ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Effective"),
                link = "logit",
                weights = ipw)

model_order_caring <- clm(rating ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Caring"),
                link = "logit",
                weights = ipw)

model_order_well_led <- clm(rating ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Well-led"),
                link = "logit",
                weights = ipw)

model_order_responsive <- clm(rating ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Responsive"),
                link = "logit",
                weights = ipw)
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
Good|Inadequate 2.144 2.299 2.715 2.404 1.713 2.338
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Inadequate|Outstanding 2.229 2.410 2.738 2.409 1.814 2.345
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Outstanding|Req improv 2.610 2.440 2.938 3.507 2.120 3.100
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
formCIC 0.250 0.045 −0.418 0.837 0.529 0.299
(0.369) (0.888) (0.363) (0.003) ** (0.027) * (0.328)
formGOV −0.231 −0.414 −0.018 −0.251 −0.153 −0.374
(0.120) (0.014) * (0.922) (0.230) (0.242) (0.040) *
formIND 0.172 0.205 0.221 −0.405 0.216 −0.137
(0.066) + (0.035) * (0.060) + (0.011) * (0.012) * (0.232)
formNPO −0.341 −0.534 −0.303 −0.101 −0.364 −0.151
(<0.001) *** (<0.001) *** (<0.001) *** (0.270) (<0.001) *** (0.046) *
categoryresidential 0.421 0.563 0.507 −0.020 0.329 0.456
(<0.001) *** (<0.001) *** (<0.001) *** (0.785) (<0.001) *** (<0.001) ***
during_covidTRUE 1.461 1.373 1.149 0.670 1.257 0.790
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
utility 0.001 0.001 0.001 0.001 0.001 0.001
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Num.Obs. 12933 12903 12452 12361 12911 12429
AIC 15923.7 12692.7 9345.3 8559.1 18116.2 12042.2
BIC 16058.1 12827.1 9479.0 8692.7 18250.6 12175.9
RMSE 1.52 1.45 1.15 0.97 1.71 1.24

Modeling on staffing

model_staff_overall <- lm(staff_level ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Overall"),
                weights = ipw)

model_staff_safe <- lm(staff_level ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Safe"),
                weights = ipw)

model_staff_effective <- lm(staff_level ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Effective"),
                weights = ipw)

model_staff_caring <- lm(staff_level ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Caring"),
                weights = ipw)

model_staff_well_led <- lm(staff_level ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Well-led"),
                weights = ipw)

model_staff_responsive <- lm(staff_level ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Responsive"),
                weights = ipw)
staff_models <- modelsummary(list("overall" = model_staff_overall, "safe" = model_staff_safe, 
                                    "effective" = model_staff_effective, "caring"= model_staff_caring, 
                                    "well-led" = model_staff_well_led, "responsive" = model_staff_responsive),
                               coef_omit = "region", exponentiate = F,
                               statistic = "({p.value}) {stars}")
staff_models
overall safe effective caring well-led responsive
(Intercept) 1.641 1.641 1.641 1.641 1.641 1.641
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
formCIC 0.487 0.487 0.487 0.487 0.487 0.487
(0.104) (0.104) (0.104) (0.104) (0.104) (0.104)
formGOV 0.700 0.700 0.700 0.700 0.700 0.700
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
formIND −0.235 −0.235 −0.235 −0.235 −0.235 −0.235
(0.017) * (0.017) * (0.017) * (0.017) * (0.017) * (0.017) *
formNPO 0.531 0.531 0.531 0.531 0.531 0.531
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
categoryresidential 0.210 0.210 0.210 0.210 0.210 0.210
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
during_covidTRUE −0.103 −0.103 −0.103 −0.103 −0.103 −0.103
(0.048) * (0.048) * (0.048) * (0.048) * (0.048) * (0.048) *
utility −0.004 −0.004 −0.004 −0.004 −0.004 −0.004
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Num.Obs. 6162 6162 6162 6162 6162 6162
R2 0.064 0.064 0.064 0.064 0.064 0.064
R2 Adj. 0.061 0.061 0.061 0.061 0.061 0.061
AIC 25042.2 25042.2 25042.2 25042.2 25042.2 25042.2
BIC 25156.5 25156.5 25156.5 25156.5 25156.5 25156.5
Log.Lik. −12504.081 −12504.081 −12504.081 −12504.081 −12504.081 −12504.081
RMSE 1.88 1.88 1.88 1.88 1.88 1.88

Modeling on staffing level direct effects

model_staff_direct_overall <- clm(rating ~ form + category + region + during_covid + utility + staff_level,
                data = filter(cqc_skills, domain == "Overall"),
                link = "logit",
                weights = ipw)

model_staff_direct_safe <- clm(rating ~ form + category + region + during_covid + utility + staff_level,
                data = filter(cqc_skills, domain == "Safe"),
                link = "logit",
                weights = ipw)

model_staff_direct_effective <- clm(rating ~ form + category + region + during_covid + utility + staff_level,
                data = filter(cqc_skills, domain == "Effective"),
                link = "logit",
                weights = ipw)

model_staff_direct_caring <- clm(rating ~ form + category + region + during_covid + utility + staff_level,
                data = filter(cqc_skills, domain == "Caring"),
                link = "logit",
                weights = ipw)

model_staff_direct_well_led <- clm(rating ~ form + category + region + during_covid + utility + staff_level,
                data = filter(cqc_skills, domain == "Well-led"),
                link = "logit",
                weights = ipw)

model_staff_direct_responsive <- clm(rating ~ form + category + region + during_covid + utility + staff_level,
                data = filter(cqc_skills, domain == "Responsive"),
                link = "logit",
                weights = ipw)
staff_direct_models <- modelsummary(list("overall" = model_staff_direct_overall, "safe" = model_staff_direct_safe, 
                                    "effective" = model_staff_direct_effective, "caring"= model_staff_direct_caring, 
                                    "well-led" = model_staff_direct_well_led, "responsive" = model_staff_direct_responsive),
                               coef_omit = "region", exponentiate = F,
                               statistic = "({p.value}) {stars}")
staff_direct_models
overall safe effective caring well-led responsive
Good|Inadequate 2.001 2.074 2.703 2.317 1.543 2.287
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Inadequate|Outstanding 2.084 2.181 2.722 2.323 1.641 2.294
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Outstanding|Req improv 2.488 2.212 2.938 3.486 1.969 3.089
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
formCIC 0.484 0.339 −0.250 0.946 0.771 0.452
(0.088) + (0.296) (0.590) (0.001) ** (0.002) ** (0.145)
formGOV −0.132 −0.261 0.068 −0.323 −0.028 −0.371
(0.393) (0.133) (0.720) (0.148) (0.834) (0.050) +
formIND 0.101 0.073 0.015 −0.455 0.161 −0.109
(0.321) (0.494) (0.913) (0.008) ** (0.082) + (0.372)
formNPO −0.251 −0.429 −0.192 −0.066 −0.312 −0.130
(<0.001) *** (<0.001) *** (0.047) * (0.502) (<0.001) *** (0.108)
categoryresidential 0.453 0.660 0.599 −0.031 0.348 0.464
(<0.001) *** (<0.001) *** (<0.001) *** (0.693) (<0.001) *** (<0.001) ***
during_covidTRUE 1.452 1.358 1.184 0.640 1.239 0.787
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
utility 0.001 0.001 0.001 0.001 0.001 0.002
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
staff_level −0.099 −0.178 −0.044 −0.017 −0.102 0.007
(<0.001) *** (<0.001) *** (0.089) + (0.454) (<0.001) *** (0.660)
Num.Obs. 11430 11402 11015 10944 11411 10997
AIC 13952.4 11020.1 8161.7 7611.2 15971.3 10666.3
BIC 14091.9 11159.6 8300.5 7749.9 16110.8 10805.1
RMSE 1.51 1.44 1.15 0.98 1.70 1.24

Modeling on turnover

model_turnover_overall <- lm(turnover ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Overall"),
                weights = ipw)

model_turnover_safe <- lm(turnover ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Safe"),
                weights = ipw)

model_turnover_effective <- lm(turnover ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Effective"),
                weights = ipw)

model_turnover_caring <- lm(turnover ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Caring"),
                weights = ipw)

model_turnover_well_led <- lm(turnover ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Well-led"),
                weights = ipw)

model_turnover_responsive <- lm(turnover ~ form + category + region + during_covid + utility,
                data = filter(cqc_skills, domain == "Responsive"),
                weights = ipw)
turnover_models <- modelsummary(list("overall" = model_turnover_overall, "safe" = model_turnover_safe, 
                                    "effective" = model_turnover_effective, "caring"= model_turnover_caring, 
                                    "well-led" = model_turnover_well_led, "responsive" = model_turnover_responsive),
                               coef_omit = "region", exponentiate = F,
                               statistic = "({p.value}) {stars}")
turnover_models
overall safe effective caring well-led responsive
(Intercept) 35.122 35.122 35.122 35.122 35.122 35.122
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
formCIC −6.468 −6.468 −6.468 −6.468 −6.468 −6.468
(0.453) (0.453) (0.453) (0.453) (0.453) (0.453)
formGOV −20.268 −20.268 −20.268 −20.268 −20.268 −20.268
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
formIND −12.288 −12.288 −12.288 −12.288 −12.288 −12.288
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
formNPO −7.427 −7.427 −7.427 −7.427 −7.427 −7.427
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
categoryresidential −5.667 −5.667 −5.667 −5.667 −5.667 −5.667
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
during_covidTRUE 2.641 2.641 2.641 2.641 2.641 2.641
(0.035) * (0.035) * (0.035) * (0.035) * (0.035) * (0.035) *
utility 0.009 0.009 0.009 0.009 0.009 0.009
(0.171) (0.171) (0.171) (0.171) (0.171) (0.171)
Num.Obs. 5265 5265 5265 5265 5265 5265
R2 0.026 0.026 0.026 0.026 0.026 0.026
R2 Adj. 0.023 0.023 0.023 0.023 0.023 0.023
AIC 54039.4 54039.4 54039.4 54039.4 54039.4 54039.4
BIC 54151.1 54151.1 54151.1 54151.1 54151.1 54151.1
Log.Lik. −27002.691 −27002.691 −27002.691 −27002.691 −27002.691 −27002.691
RMSE 39.27 39.27 39.27 39.27 39.27 39.27

Modeling on turnover rate direct effects

model_turnover_direct_overall <- clm(rating ~ form + category + region + during_covid + utility + turnover,
                data = filter(cqc_skills, domain == "Overall"),
                link = "logit",
                weights = ipw)

model_turnover_direct_safe <- clm(rating ~ form + category + region + during_covid + utility + turnover,
                data = filter(cqc_skills, domain == "Safe"),
                link = "logit",
                weights = ipw)

model_turnover_direct_effective <- clm(rating ~ form + category + region + during_covid + utility + turnover,
                data = filter(cqc_skills, domain == "Effective"),
                link = "logit",
                weights = ipw)

model_turnover_direct_caring <- clm(rating ~ form + category + region + during_covid + utility + turnover,
                data = filter(cqc_skills, domain == "Caring"),
                link = "logit",
                weights = ipw)

model_turnover_direct_well_led <- clm(rating ~ form + category + region + during_covid + utility + turnover,
                data = filter(cqc_skills, domain == "Well-led"),
                link = "logit",
                weights = ipw)

model_turnover_direct_responsive <- clm(rating ~ form + category + region + during_covid + utility + turnover,
                data = filter(cqc_skills, domain == "Responsive"),
                link = "logit",
                weights = ipw)
turnover_direct_models <- modelsummary(list("overall" = model_turnover_direct_overall, "safe" = model_turnover_direct_safe, 
                                    "effective" = model_turnover_direct_effective, "caring"= model_turnover_direct_caring, 
                                    "well-led" = model_turnover_direct_well_led, "responsive" = model_turnover_direct_responsive),
                               coef_omit = "region", exponentiate = F,
                               statistic = "({p.value}) {stars}")
turnover_direct_models
overall safe effective caring well-led responsive
Good|Inadequate 2.164 2.333 2.757 2.339 1.770 2.359
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Inadequate|Outstanding 2.238 2.428 2.776 2.347 1.867 2.369
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Outstanding|Req improv 2.638 2.463 3.000 3.560 2.190 3.168
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
formCIC 0.666 0.576 −0.703 0.479 0.885 0.475
(0.059) + (0.138) (0.356) (0.278) (0.005) ** (0.254)
formGOV −0.150 −0.339 0.082 −0.197 −0.064 −0.249
(0.321) (0.045) * (0.656) (0.352) (0.630) (0.176)
formIND 0.259 0.253 0.062 −0.399 0.254 0.039
(0.014) * (0.021) * (0.664) (0.027) * (0.009) ** (0.758)
formNPO −0.248 −0.439 −0.183 0.049 −0.331 −0.038
(0.003) ** (<0.001) *** (0.084) + (0.638) (<0.001) *** (0.667)
categoryresidential 0.394 0.562 0.603 −0.026 0.319 0.469
(<0.001) *** (<0.001) *** (<0.001) *** (0.756) (<0.001) *** (<0.001) ***
during_covidTRUE 1.425 1.353 1.144 0.590 1.237 0.697
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
utility 0.001 0.002 0.001 0.001 0.001 0.001
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
turnover 0.002 0.002 0.002 0.000 0.002 0.002
(<0.001) *** (0.008) ** (0.014) * (0.626) (<0.001) *** (0.004) **
Num.Obs. 9774 9748 9403 9348 9757 9388
AIC 11757.9 9416.8 6935.9 6325.2 13419.8 8838.4
BIC 11894.5 9553.3 7071.8 6460.9 13556.3 8974.2
RMSE 1.49 1.44 1.14 0.95 1.67 1.21