Ordered Logit Models with Average Marginal Effects

Read Me

This is the data analysis RMarkdown file for the QCQ data after ISIRC 2022.

Mainly changes: incorporating the contol of before and after COVID

Please read the subtitles and notes added as normal text in this document. Blocks with darker backgrounds are code chunks, mostly with their outputs. I tried to add #comments within the chunck just before the code line to explain the purpose of the code line. Meanwhile I explain the purpose of the section in the texts, so that is where to find information to give you fuller pictures.

# 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

Data Preparation

First, import the sampled and coded data set

#import the raw data file
socialcare_raw <- read_csv("data/sample_new_cleaned.csv")

Assign orders to the ordinal level variables and name the organizational form in a reader-friendly way.

#select relevant columns, rename and relabel 
socialcare <- socialcare_raw %>% 
  # recode legal form types to be more readable / easier to present
  mutate(location_id = factor(location_id),
         form = case_when(form_num == 1 ~ "FPO",
                          form_num == 2 ~ "NPO",
                          form_num == 3 ~ "GOV",
                          form_num == 4 ~ "CIC",
                          form_num == 5 ~ "IND",
                          TRUE ~ NA_character_),
         inherited = ifelse(inherited == "Y", TRUE, FALSE),
         rating = recode(rating, 
                         "Insufficient evidence to rate" = "NA",
                         "Requires improvement" = "Req improv"),
         publication_date = dmy(publication_date)) %>% 
  # set the order of the values in the factors 
  mutate(form = ordered(form, levels = c("FPO", "NPO", "GOV", "CIC", "IND")),
         
         # assume the order of the ratings as follows but need to double check with the source 
         rating = ordered(rating, levels = c("Inadequate","Req improv", "Good", "Outstanding"))) %>% 
  
  # creating a new dummy variable for facility category
  mutate(category = case_when(primary_cat == "Community based adult social care services" ~ "community",
                              primary_cat == "Residential social care" ~ "residential",
                              TRUE ~ as.character(NA)),
         # deriving year column and dummy variable for before_covid
         year = year(publication_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)) %>% 
  # derive the rating dummy
  mutate(rating_higher = ifelse(rating_num > 2, 1, 0))

# show first several rows of the data set derived 
head(socialcare)

## # A tibble: 6 × 40
##    ...1 locati…¹ Locat…² Locat…³ Care …⁴ prima…⁵ Locat…⁶ Locat…⁷ Locat…⁸ Locat…⁹
##   <dbl> <fct>    <chr>   <chr>   <chr>   <chr>   <chr>   <chr>   <chr>   <chr>  
## 1     1 1-10000… VNH4P   69 Ten… N       Commun… 69 Ten… Warmsw… Doncas… DN4 9PE
## 2     2 1-10000… VNH4P   69 Ten… N       Commun… 69 Ten… Warmsw… Doncas… DN4 9PE
## 3     3 1-10000… VNH4P   69 Ten… N       Commun… 69 Ten… Warmsw… Doncas… DN4 9PE
## 4     4 1-10000… VNH4P   69 Ten… N       Commun… 69 Ten… Warmsw… Doncas… DN4 9PE
## 5     5 1-10000… VNH4P   69 Ten… N       Commun… 69 Ten… Warmsw… Doncas… DN4 9PE
## 6     6 1-10000… VNH4P   69 Ten… N       Commun… 69 Ten… Warmsw… Doncas… DN4 9PE
## # … with 30 more variables: `Location Local Authority` <chr>, region <chr>,
## #   `Location NHS Region` <chr>, `Location ONSPD CCG Code` <chr>,
## #   `Location ONSPD CCG` <chr>, `Location Commissioning CCG Code` <lgl>,
## #   `Location Commissioning CCG Name` <lgl>,
## #   `Service / Population Group` <chr>, domain <chr>, rating <ord>,
## #   publication_date <date>, `Report Type` <chr>, inherited <lgl>, URL <chr>,
## #   `Provider ID` <chr>, location_type <chr>, form_num <dbl>, …

whole model without control

model_order_overall <- clm(rating ~ form + category + region + inherited,
                data = filter(socialcare, domain == "Overall"),
                link = "logit")
model_order_safe <- clm(rating ~ form + category + region + inherited,
                data = filter(socialcare, domain == "Safe"),
                link = "logit")
model_order_effective <- clm(rating ~ form + category + region + inherited,
                data = filter(socialcare, domain == "Effective"),
                link = "logit")
model_order_caring <- clm(rating ~ form + category + region + inherited,
                data = filter(socialcare, domain == "Caring"),
                link = "logit")
model_order_well_led <- clm(rating ~ form + category + region + inherited,
                data = filter(socialcare, domain == "Well-led"),
                link = "logit")
model_order_responsive <- clm(rating ~ form + category + region + inherited,
                data = filter(socialcare, 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}")

## Warning in !is.null(rmarkdown::metadata$output) && rmarkdown::metadata$output
## %in% : 'length(x) = 3 > 1' in coercion to 'logical(1)'

ordinal_models

	overall	safe	effective	caring	well-led	responsive
Inadequate\|Req improv	−4.449***	−4.253***	−5.932***	−7.094***	−4.102***	−6.856***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
Req improv\|Good	−1.702***	−1.444***	−2.079***	−3.431***	−1.334***	−2.365***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
Good\|Outstanding	3.144***	5.856***	4.395***	3.024***	3.178***	2.980***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
formNPO	0.457***	0.580***	0.413***	0.316***	0.475***	0.507***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
formGOV	0.429**	0.685***	0.397*	0.328+	0.445***	0.415**
	(0.001) **	(0.000) ***	(0.029) *	(0.066) +	(0.000) ***	(0.006) **
formCIC	0.586*	0.797*	0.188	0.819*	0.478+	0.306
	(0.045) *	(0.040) *	(0.622)	(0.013) *	(0.085) +	(0.340)
formIND	−0.183*	−0.117	−0.275*	−0.050	−0.162*	−0.057
	(0.036) *	(0.197)	(0.011) *	(0.710)	(0.045) *	(0.592)
categoryresidential	−0.387***	−0.423***	−0.374***	−0.574***	−0.302***	−0.237***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
inheritedTRUE	−0.108	−0.140	−0.288**	−0.397**	−0.076	−0.251*
	(0.238)	(0.146)	(0.010) **	(0.003) **	(0.372)	(0.017) *
Num.Obs.	12942	12966	12931	12927	12964	12953
AIC	16608.8	13440.5	10365.2	8969.0	18898.5	12726.4
BIC	16735.8	13567.5	10492.1	9096.0	19025.5	12853.4
RMSE	2.26	2.12	2.13	2.21	2.26	2.27

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}", output = "data/ordinal_models_exp.docx")
ordinal_models_exp

## [1] "data/ordinal_models_exp.docx"

whole model with control

model_order_overall_covid <- clm(rating ~ form + after_covid + 
                                   category + region + inherited,
                data = filter(socialcare, domain == "Overall"),
                link = "logit")
model_order_safe_covid <- clm(rating ~ form + after_covid +
                                category + region + inherited,
                data = filter(socialcare, domain == "Safe"),
                link = "logit")
model_order_effective_covid <- clm(rating ~ form + after_covid +
                                     category + region + inherited,
                data = filter(socialcare, domain == "Effective"),
                link = "logit")
model_order_caring_covid <- clm(rating ~ form + after_covid + 
                                  category + region + inherited,
                data = filter(socialcare, domain == "Caring"),
                link = "logit")
model_order_well_led_covid <- clm(rating ~ form + after_covid +
                                    category + region + inherited,
                data = filter(socialcare, domain == "Well-led"),
                link = "logit")
model_order_responsive_covid <- clm(rating ~ form + after_covid +
                                      category + region + inherited,
                data = filter(socialcare, 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	−5.398***	−5.000***	−6.333***	−7.446***	−4.796***	−7.122***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
Req improv\|Good	−2.498***	−2.079***	−2.442***	−3.765***	−1.905***	−2.611***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
Good\|Outstanding	2.794***	5.545***	4.214***	2.847***	2.866***	2.833***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
formNPO	0.307***	0.450***	0.350***	0.274**	0.357***	0.467***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.002) **	(0.000) ***	(0.000) ***
formGOV	0.290*	0.561**	0.337+	0.290	0.316*	0.378*
	(0.036) *	(0.001) **	(0.065) +	(0.107)	(0.014) *	(0.012) *
formCIC	0.494+	0.756+	0.183	0.802*	0.408	0.296
	(0.097) +	(0.056) +	(0.633)	(0.015) *	(0.146)	(0.354)
formIND	−0.339***	−0.248**	−0.344**	−0.097	−0.282***	−0.097
	(0.000) ***	(0.009) **	(0.002) **	(0.469)	(0.001) ***	(0.364)
after_covidTRUE	−1.776***	−1.496***	−1.132***	−1.042***	−1.412***	−0.843***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
categoryresidential	−0.377***	−0.420***	−0.386***	−0.587***	−0.294***	−0.244***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
inheritedTRUE	−0.322***	−0.335***	−0.398***	−0.460***	−0.238**	−0.306**
	(0.001) ***	(0.001) ***	(0.000) ***	(0.001) ***	(0.007) **	(0.004) **
Num.Obs.	12942	12966	12931	12927	12964	12953
AIC	15296.7	12507.3	10045.6	8808.0	17890.6	12548.2
BIC	15431.1	12641.8	10180.0	8942.4	18025.1	12682.7
RMSE	2.24	2.10	2.13	2.21	2.25	2.27

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.005***	0.007***	0.002***	0.001***	0.008***	0.001***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
Req improv\|Good	0.082***	0.125***	0.087***	0.023***	0.149***	0.073***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
Good\|Outstanding	16.352***	256.051***	67.643***	17.239***	17.568***	16.988***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
formNPO	1.359***	1.568***	1.419***	1.316**	1.428***	1.595***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.002) **	(0.000) ***	(0.000) ***
formGOV	1.336*	1.753**	1.401+	1.337	1.371*	1.459*
	(0.036) *	(0.001) **	(0.065) +	(0.107)	(0.014) *	(0.012) *
formCIC	1.638+	2.130+	1.200	2.231*	1.504	1.345
	(0.097) +	(0.056) +	(0.633)	(0.015) *	(0.146)	(0.354)
formIND	0.713***	0.780**	0.709**	0.908	0.754***	0.908
	(0.000) ***	(0.009) **	(0.002) **	(0.469)	(0.001) ***	(0.364)
after_covidTRUE	0.169***	0.224***	0.323***	0.353***	0.244***	0.430***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
categoryresidential	0.686***	0.657***	0.680***	0.556***	0.746***	0.784***
	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***	(0.000) ***
inheritedTRUE	0.725***	0.715***	0.672***	0.632***	0.788**	0.736**
	(0.001) ***	(0.001) ***	(0.000) ***	(0.001) ***	(0.007) **	(0.004) **
Num.Obs.	12942	12966	12931	12927	12964	12953
AIC	15296.7	12507.3	10045.6	8808.0	17890.6	12548.2
BIC	15431.1	12641.8	10180.0	8942.4	18025.1	12682.7
RMSE	2.24	2.10	2.13	2.21	2.25	2.27

Ordered Logit Models with Average Marginal Effects

Meng Ye

10/28/2022

Read Me

Data Preparation

whole model without control

whole model with control