# 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) %>% 
  mutate(rating = ordered(rating, levels = c("Inadequate","Req improv", "Good", "Outstanding"))) %>% 
     filter(form!="NA" & rating!="NA")

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} \]

Modeling on rating with staff control

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

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

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


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

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

staff_direct_order <- clm(rating ~ form + staff_level + category + region + during_covid,
                data = filter(cqc_skills, domain == "Overall"),
                link = "logit")
staff_direct_order_weighted <- clm(rating ~ form + staff_level + category + region + during_covid,
                data = filter(cqc_skills, domain == "Overall"),
                link = "logit",
                weights = ipw)
staff_models <- modelsummary(list("full" = model_order, "full weighted" = model_order_weighted, 
                                  "full weighted 2" = model_order_weighted2, 
                                    "staffing" = model_staff, "staffing weighted"= model_staff_weighted, 
                                    "direct" = staff_direct_order, "direct weighted" = staff_direct_order_weighted),
                               coef_omit = "region", exponentiate = F,
                               statistic = "({p.value}) {stars}")
staff_models
full full weighted  full weighted 2 staffing staffing weighted direct direct weighted
Inadequate|Req improv −5.830 −5.584 −5.561 −5.782 −5.544
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Req improv|Good −2.862 −2.732 −2.719 −2.820 −2.703
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Good|Outstanding 2.383 2.521 2.522 2.419 2.555
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
formCIC 0.368 0.348 0.376 0.584 0.547 0.284 0.241
(0.336) (0.244) (0.205) (0.043) * (0.075) + (0.480) (0.445)
formGOV 0.125 0.264 0.276 0.783 0.785 0.006 0.134
(0.391) (0.054) + (0.044) * (<0.001) *** (<0.001) *** (0.966) (0.344)
formIND −0.184 −0.249 −0.221 −0.179 −0.180 −0.128 −0.189
(0.259) (0.008) ** (0.018) * (0.168) (0.084) + (0.464) (0.061) +
formNPO 0.273 0.336 0.344 0.572 0.571 0.215 0.269
(0.001) ** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (0.018) * (<0.001) ***
categoryresidential −0.491 −0.472 −0.480 0.462 0.473 −0.543 −0.538
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
during_covidTRUE −1.687 −1.789 −1.785 −0.169 −0.166 −1.665 −1.771
(<0.001) *** (<0.001) *** (<0.001) *** (0.004) ** (0.003) ** (<0.001) *** (<0.001) ***
(Intercept) 1.300 1.283
(<0.001) *** (<0.001) ***
staff_level 0.048 0.058
(0.005) ** (<0.001) ***
Num.Obs. 6504 12933 12941 5767 5767 5767 11430
R2 0.037 0.036
R2 Adj. 0.035 0.033
AIC 7518.8 15461.8 15516.9 24025.6 23772.0 6650.4 13575.8
BIC 7634.1 15588.7 15643.9 24132.2 23878.5 6770.3 13708.0
Log.Lik. −11996.820 −11869.988
RMSE 2.27 2.27 2.27 1.94 1.94 2.27 2.28

Modeling on rating with turnover

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

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

turnover_direct_order <- clm(rating ~ form + turnover + category + region + during_covid,
                data = filter(cqc_skills, domain == "Overall"),
                link = "logit")
turnover_direct_order_weighted <- clm(rating ~ form + turnover + category + region + during_covid,
                data = filter(cqc_skills, domain == "Overall"),
                link = "logit",
                weights = ipw)
staff_models <- modelsummary(list("full" = model_order, "full weighted" = model_order_weighted, 
                                    "turnover" = model_turnover, "turnover weighted"= model_turnover_weighted, 
                                    "direct" = turnover_direct_order, "direct weighted" = turnover_direct_order_weighted),
                               coef_omit = "region", exponentiate = F,
                               statistic = "({p.value}) {stars}")
staff_models
full full weighted turnover turnover weighted direct direct weighted
Inadequate|Req improv −5.830 −5.584 −5.964 −5.733
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Req improv|Good −2.862 −2.732 −2.916 −2.809
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
Good|Outstanding 2.383 2.521 2.366 2.489
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
formCIC 0.368 0.348 −5.100 −5.935 0.219 0.130
(0.336) (0.244) (0.501) (0.486) (0.683) (0.756)
formGOV 0.125 0.264 −19.803 −20.050 0.065 0.199
(0.391) (0.054) + (<0.001) *** (<0.001) *** (0.670) (0.157)
formIND −0.184 −0.249 −11.525 −11.785 −0.165 −0.248
(0.259) (0.008) ** (<0.001) *** (<0.001) *** (0.376) (0.021) *
formNPO 0.273 0.336 −7.255 −7.369 0.260 0.315
(0.001) ** (<0.001) *** (<0.001) *** (<0.001) *** (0.010) ** (<0.001) ***
categoryresidential −0.491 −0.472 −5.656 −6.003 −0.561 −0.554
(<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) *** (<0.001) ***
during_covidTRUE −1.687 −1.789 3.175 3.249 −1.642 −1.754
(<0.001) *** (<0.001) *** (0.012) * (0.011) * (<0.001) *** (<0.001) ***
(Intercept) 35.652 35.499
(<0.001) *** (<0.001) ***
turnover −0.001 −0.001
(0.248) (0.127)
Num.Obs. 6504 12933 4939 4939 4939 9774
R2 0.031 0.026
R2 Adj. 0.029 0.023
AIC 7518.8 15461.8 50234.9 50613.3 5599.5 11414.2
BIC 7634.1 15588.7 50339.0 50717.3 5716.6 11543.5
Log.Lik. −25101.448 −25290.627
RMSE 2.27 2.27 38.99 38.99 2.27 2.27