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>, …
socialcare %>% 
  group_by(year, form) %>% 
  drop_na(rating_num) %>% 
  summarise(count = n(),
            perform = mean(rating_num))
## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.
## # A tibble: 34 × 4
## # Groups:   year [6]
##     year form  count perform
##    <dbl> <ord> <int>   <dbl>
##  1  2016 FPO     168    2.96
##  2  2016 NPO      36    2.89
##  3  2016 GOV      12    3   
##  4  2016 CIC       6    2.83
##  5  2016 IND      30    3   
##  6  2017 FPO    4251    3.02
##  7  2017 NPO    1171    3.04
##  8  2017 GOV     210    3.05
##  9  2017 CIC      30    3.07
## 10  2017 IND     376    3.00
## # … with 24 more rows

Descriptive analysis

table1 <- datasummary((category + inherited + before_covid + after_covid +  region) ~ form,
                      data = socialcare, fmt = 0)
## Warning in !is.null(rmarkdown::metadata$output) && rmarkdown::metadata$output
## %in% : 'length(x) = 3 > 1' in coercion to 'logical(1)'
table1
FPO NPO GOV CIC IND
category community 22454 3365 1194 390 798
residential 41721 9022 1231 114 4355
inherited FALSE 60318 11772 2185 492 5069
TRUE 3869 627 246 12 84
before_covid FALSE 18976 2588 510 115 1197
TRUE 45211 9811 1921 389 3956
after_covid FALSE 45211 9811 1921 389 3956
TRUE 18976 2588 510 115 1197
region East Midlands 6510 896 306 24 396
East of England 7451 1128 181 108 509
London 7055 1180 246 60 600
North East 2646 654 186 18 144
North West 7911 1448 348 132 774
South East 11835 2400 408 36 930
South West 7304 2070 180 36 858
West Midlands 7547 1303 258 42 498
Yorkshire and The Humber 5928 1320 318 48 444

Running OLS models

The OLS models and the ordered logit models can be written as follows

\(rating_{numerical} = \beta_0 + \beta_1form + \beta_2category+ \beta_3region + \beta_4inherited + u\)

\(log-odds(rating_{ordinal} \leq j) = \beta_{j0} + \beta_1form + \beta_2category+ \beta_3region + \beta_4inherited + u\)

In this section, we first run the OLS models. In the OLS models, we kind of “cheat” R by treating the four rating levels with orders as if they are numbers 1-4. There are the flowing reasons that we report the results from OLS models, even though the more suitable methods should be ordered logit models, about which we will discuss in a while.

  1. The purpose of fitting OLS models is to use them as benchmarks.
  2. Since there are issues like heteroscedasticity, the standard errors calculated are not reliable. But the correlation relationships between the independent variables and dependent variables are still true. So the results are still informative
  3. Plus, compared with the ordered logit models we run later, the results are more straightforward, and more easily give us intuition about how different legal forms of social care providers impact the service quality ratings.
  4. The OLS models are intended to be compared with the ordered logit models. As shown later, the results are generally consistent between the two model families, confirming that our model specification is robust between different models.

OLS with the sub-domain ratings one by one

# run the model loops with nest() method
models_ols <- socialcare %>% 
  group_by(domain) %>% 
  nest()%>% 
  mutate(ols_models = map(data, 
                          ~lm(rating_num ~ form + 
                              category + region + inherited , 
                              data = .x))) %>% 
  mutate(results = map(ols_models, ~tidy(.x, conf.int = TRUE))) %>% 
  pull(ols_models, name = domain)
# run the model loops with nest() method
# add before_covid as control
# add also interaction term
table_ols <- modelsummary(models_ols, statistic = "({p.value}) {stars}")
table_ols
Safe Effective Caring Responsive Well-led Overall
(Intercept) 2.789*** 2.892*** 3.016*** 2.961*** 2.814*** 2.873***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
formNPO 0.082*** 0.040*** 0.027*** 0.064*** 0.100*** 0.082***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
formGOV 0.097*** 0.039* 0.029+ 0.053** 0.097*** 0.081***
(0.000) *** (0.023) * (0.061) + (0.006) ** (0.000) *** (0.001) ***
formCIC 0.100* 0.019 0.079* 0.039 0.109+ 0.104*
(0.030) * (0.604) (0.014) * (0.339) (0.055) + (0.043) *
formIND −0.019 −0.033** −0.004 −0.007 −0.040* −0.036*
(0.214) (0.007) ** (0.712) (0.621) (0.033) * (0.035) *
categoryresidential −0.064*** −0.037*** −0.047*** −0.030*** −0.067*** −0.070***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionEast of England 0.047** 0.041** 0.015 0.035* 0.030 0.043*
(0.004) ** (0.002) ** (0.189) (0.016) * (0.132) (0.017) *
regionLondon 0.048** 0.054*** 0.003 −0.003 0.025 0.019
(0.004) ** (0.000) *** (0.807) (0.820) (0.206) (0.308)
regionNorth East 0.123*** 0.095*** 0.069*** 0.082*** 0.155*** 0.133***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionNorth West 0.067*** 0.064*** 0.022* 0.040** 0.052** 0.052**
(0.000) *** (0.000) *** (0.050) * (0.004) ** (0.007) ** (0.003) **
regionSouth East 0.098*** 0.049*** 0.040*** 0.051*** 0.057** 0.062***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.002) ** (0.000) ***
regionSouth West 0.124*** 0.094*** 0.085*** 0.091*** 0.144*** 0.142***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionWest Midlands 0.019 0.020 −0.015 0.000 −0.075*** −0.030+
(0.233) (0.114) (0.173) (0.982) (0.000) *** (0.092) +
regionYorkshire and The Humber −0.012 0.031* 0.036** 0.008 −0.033 −0.025
(0.467) (0.020) * (0.002) ** (0.576) (0.115) (0.186)
inheritedTRUE −0.018 −0.032** −0.032** −0.032* −0.018 −0.020
(0.248) (0.008) ** (0.003) ** (0.020) * (0.346) (0.260)
Num.Obs. 12966 12931 12927 12953 12964 12942
R2 0.023 0.013 0.017 0.013 0.025 0.023
R2 Adj. 0.021 0.012 0.016 0.012 0.024 0.022
AIC 14402.5 8181.6 5099.3 11037.1 19530.8 16998.2
BIC 14522.0 8301.1 5218.8 11156.6 19650.3 17117.7
Log.Lik. −7185.242 −4074.798 −2533.650 −5502.563 −9749.405 −8483.091
RMSE 0.42 0.33 0.29 0.37 0.51 0.47 
# run the model loops with nest() method
# add after_covid as control
models_ols_covid <- socialcare %>% 
  group_by(domain) %>% 
  nest()%>% 
  mutate(ols_models = map(data, 
                          ~lm(rating_num ~ form + after_covid + 
                              category + region + inherited, 
                              data = .x))) %>% 
  mutate(results = map(ols_models, ~tidy(.x, conf.int = TRUE))) %>% 
  pull(ols_models, name = domain) 
table_ols_covid <- modelsummary(models_ols_covid, statistic = "({p.value}) {stars}")
table_ols_covid
Safe Effective Caring Responsive Well-led Overall
(Intercept) 2.879*** 2.925*** 3.036*** 2.986*** 2.923*** 2.982***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
formNPO 0.056*** 0.032*** 0.023** 0.058*** 0.069*** 0.050***
(0.000) *** (0.000) *** (0.002) ** (0.000) *** (0.000) *** (0.000) ***
formGOV 0.069** 0.032+ 0.025 0.048* 0.065* 0.050*
(0.001) ** (0.058) + (0.102) (0.013) * (0.011) * (0.028) *
formCIC 0.084+ 0.016 0.079* 0.037 0.085 0.080+
(0.058) + (0.661) (0.014) * (0.357) (0.118) (0.100) +
formIND −0.038* −0.039** −0.008 −0.012 −0.062*** −0.059***
(0.011) * (0.001) ** (0.454) (0.387) (0.001) *** (0.000) ***
after_covidTRUE −0.271*** −0.135*** −0.083*** −0.107*** −0.329*** −0.332***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
categoryresidential −0.060*** −0.037*** −0.048*** −0.031*** −0.061*** −0.063***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionEast of England 0.038* 0.040** 0.014 0.034* 0.019 0.032+
(0.016) * (0.002) ** (0.212) (0.016) * (0.325) (0.058) +
regionLondon 0.045** 0.055*** 0.002 −0.001 0.024 0.018
(0.004) ** (0.000) *** (0.834) (0.919) (0.209) (0.303)
regionNorth East 0.133*** 0.104*** 0.074*** 0.089*** 0.169*** 0.150***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionNorth West 0.062*** 0.063*** 0.020+ 0.039** 0.047* 0.046**
(0.000) *** (0.000) *** (0.068) + (0.006) ** (0.012) * (0.006) **
regionSouth East 0.085*** 0.047*** 0.038*** 0.050*** 0.042* 0.047**
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.014) * (0.002) **
regionSouth West 0.110*** 0.092*** 0.082*** 0.088*** 0.128*** 0.126***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionWest Midlands 0.007 0.018 −0.017 −0.002 −0.090*** −0.045**
(0.662) (0.142) (0.125) (0.892) (0.000) *** (0.008) **
regionYorkshire and The Humber −0.016 0.031* 0.036** 0.008 −0.036+ −0.029
(0.334) (0.018) * (0.003) ** (0.596) (0.066) + (0.105)
inheritedTRUE −0.045** −0.042*** −0.037*** −0.038** −0.051** −0.053**
(0.003) ** (0.001) *** (0.001) *** (0.005) ** (0.006) ** (0.001) **
Num.Obs. 12966 12931 12927 12953 12964 12942
R2 0.104 0.041 0.030 0.027 0.107 0.124
R2 Adj. 0.103 0.039 0.029 0.026 0.106 0.123
AIC 13271.6 7815.4 4933.2 10857.2 18399.2 15587.3
BIC 13398.6 7942.4 5060.2 10984.2 18526.2 15714.3
Log.Lik. −6618.794 −3890.704 −2449.623 −5411.619 −9182.592 −7776.674
RMSE 0.40 0.33 0.29 0.37 0.49 0.44 
# run the model loops with nest() method
# add after_covid as control
models_ols_covid_b <- socialcare %>% 
  group_by(domain) %>% 
  nest()%>% 
  mutate(ols_models = map(data, 
                          ~lm(rating_num ~ form + before_covid + 
                              category + region + inherited, 
                              data = .x))) %>% 
  mutate(results = map(ols_models, ~tidy(.x, conf.int = TRUE))) %>% 
  pull(ols_models, name = domain) 
table_ols_covid_b <- modelsummary(models_ols_covid_b, statistic = "({p.value}) {stars}")
table_ols_covid_b
Safe Effective Caring Responsive Well-led Overall
(Intercept) 2.609*** 2.790*** 2.953*** 2.879*** 2.594*** 2.650***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
formNPO 0.056*** 0.032*** 0.023** 0.058*** 0.069*** 0.050***
(0.000) *** (0.000) *** (0.002) ** (0.000) *** (0.000) *** (0.000) ***
formGOV 0.069** 0.032+ 0.025 0.048* 0.065* 0.050*
(0.001) ** (0.058) + (0.102) (0.013) * (0.011) * (0.028) *
formCIC 0.084+ 0.016 0.079* 0.037 0.085 0.080+
(0.058) + (0.661) (0.014) * (0.357) (0.118) (0.100) +
formIND −0.038* −0.039** −0.008 −0.012 −0.062*** −0.059***
(0.011) * (0.001) ** (0.454) (0.387) (0.001) *** (0.000) ***
before_covidTRUE 0.271*** 0.135*** 0.083*** 0.107*** 0.329*** 0.332***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
categoryresidential −0.060*** −0.037*** −0.048*** −0.031*** −0.061*** −0.063***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionEast of England 0.038* 0.040** 0.014 0.034* 0.019 0.032+
(0.016) * (0.002) ** (0.212) (0.016) * (0.325) (0.058) +
regionLondon 0.045** 0.055*** 0.002 −0.001 0.024 0.018
(0.004) ** (0.000) *** (0.834) (0.919) (0.209) (0.303)
regionNorth East 0.133*** 0.104*** 0.074*** 0.089*** 0.169*** 0.150***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionNorth West 0.062*** 0.063*** 0.020+ 0.039** 0.047* 0.046**
(0.000) *** (0.000) *** (0.068) + (0.006) ** (0.012) * (0.006) **
regionSouth East 0.085*** 0.047*** 0.038*** 0.050*** 0.042* 0.047**
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.014) * (0.002) **
regionSouth West 0.110*** 0.092*** 0.082*** 0.088*** 0.128*** 0.126***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionWest Midlands 0.007 0.018 −0.017 −0.002 −0.090*** −0.045**
(0.662) (0.142) (0.125) (0.892) (0.000) *** (0.008) **
regionYorkshire and The Humber −0.016 0.031* 0.036** 0.008 −0.036+ −0.029
(0.334) (0.018) * (0.003) ** (0.596) (0.066) + (0.105)
inheritedTRUE −0.045** −0.042*** −0.037*** −0.038** −0.051** −0.053**
(0.003) ** (0.001) *** (0.001) *** (0.005) ** (0.006) ** (0.001) **
Num.Obs. 12966 12931 12927 12953 12964 12942
R2 0.104 0.041 0.030 0.027 0.107 0.124
R2 Adj. 0.103 0.039 0.029 0.026 0.106 0.123
AIC 13271.6 7815.4 4933.2 10857.2 18399.2 15587.3
BIC 13398.6 7942.4 5060.2 10984.2 18526.2 15714.3
Log.Lik. −6618.794 −3890.704 −2449.623 −5411.619 −9182.592 −7776.674
RMSE 0.40 0.33 0.29 0.37 0.49 0.44 
# run the model loops with nest() method
# add after_covid as control
# and interaction term
models_ols_inter <- socialcare %>% 
  group_by(domain) %>% 
  nest()%>% 
  mutate(ols_models = map(data, 
                          ~lm(rating_num ~ form * after_covid + 
                              category + region + inherited, 
                              data = .x))) %>% 
  mutate(results = map(ols_models, ~tidy(.x, conf.int = TRUE))) 

models_ols_inter_named <- models_ols_inter %>% 
  pull(ols_models, name = domain)
table_ols_inter <- modelsummary(models_ols_inter_named, statistic = "({p.value}) {stars}")
table_ols_inter
Safe Effective Caring Responsive Well-led Overall
(Intercept) 2.881*** 2.927*** 3.038*** 2.989*** 2.925*** 2.986***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
formNPO 0.044*** 0.019* 0.014+ 0.043*** 0.056*** 0.030*
(0.000) *** (0.030) * (0.073) + (0.000) *** (0.000) *** (0.018) *
formGOV 0.035 0.009 0.012 0.035+ 0.034 0.011
(0.134) (0.621) (0.464) (0.096) + (0.244) (0.664)
formCIC 0.036 0.001 0.026 0.032 0.047 0.028
(0.483) (0.986) (0.481) (0.481) (0.447) (0.614)
formIND −0.011 −0.026* −0.006 −0.008 −0.039+ −0.038*
(0.518) (0.050) * (0.632) (0.614) (0.060) + (0.040) *
after_covidTRUE −0.276*** −0.143*** −0.092*** −0.119*** −0.336*** −0.344***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
categoryresidential −0.059*** −0.037*** −0.048*** −0.031*** −0.061*** −0.063***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionEast of England 0.037* 0.039** 0.014 0.034* 0.018 0.032+
(0.017) * (0.002) ** (0.230) (0.017) * (0.342) (0.064) +
regionLondon 0.045** 0.055*** 0.002 −0.002 0.024 0.018
(0.004) ** (0.000) *** (0.843) (0.893) (0.207) (0.306)
regionNorth East 0.132*** 0.103*** 0.074*** 0.089*** 0.168*** 0.149***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionNorth West 0.062*** 0.063*** 0.021+ 0.039** 0.047* 0.046**
(0.000) *** (0.000) *** (0.063) + (0.005) ** (0.012) * (0.006) **
regionSouth East 0.085*** 0.046*** 0.038*** 0.049*** 0.042* 0.047**
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.015) * (0.002) **
regionSouth West 0.109*** 0.091*** 0.082*** 0.088*** 0.128*** 0.125***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionWest Midlands 0.006 0.018 −0.018 −0.002 −0.091*** −0.046**
(0.691) (0.157) (0.119) (0.867) (0.000) *** (0.007) **
regionYorkshire and The Humber −0.016 0.031* 0.035** 0.007 −0.037+ −0.029
(0.314) (0.021) * (0.003) ** (0.623) (0.063) + (0.100)
inheritedTRUE −0.045** −0.042*** −0.037*** −0.038** −0.051** −0.054**
(0.003) ** (0.001) *** (0.001) *** (0.005) ** (0.006) ** (0.001) **
formNPO × after_covidTRUE 0.055* 0.070*** 0.049* 0.086*** 0.056+ 0.088***
(0.020) * (0.001) *** (0.010) * (0.000) *** (0.050) + (0.001) ***
formGOV × after_covidTRUE 0.160** 0.124** 0.071+ 0.067 0.145* 0.176**
(0.002) ** (0.005) ** (0.076) + (0.181) (0.019) * (0.001) **
formCIC × after_covidTRUE 0.192+ 0.070 0.248** 0.023 0.156 0.215+
(0.059) + (0.424) (0.002) ** (0.818) (0.217) (0.058) +
formIND × after_covidTRUE −0.109** −0.075* −0.017 −0.027 −0.093* −0.084*
(0.001) ** (0.015) * (0.542) (0.446) (0.024) * (0.023) *
Num.Obs. 12966 12931 12927 12953 12964 12942
R2 0.106 0.043 0.032 0.028 0.108 0.126
R2 Adj. 0.105 0.041 0.030 0.027 0.107 0.125
AIC 13249.0 7796.5 4921.7 10849.4 18390.5 15564.4
BIC 13405.9 7953.3 5078.5 11006.3 18547.3 15721.2
Log.Lik. −6603.505 −3877.241 −2439.846 −5403.713 −9174.233 −7761.182
RMSE 0.40 0.33 0.29 0.37 0.49 0.44 
# run the model loops with nest() method
# add after_covid as control
# and interaction term
models_ols_inter_b <- socialcare %>% 
  group_by(domain) %>% 
  nest()%>% 
  mutate(ols_models = map(data, 
                          ~lm(rating_num ~ form * before_covid + 
                              category + region + inherited, 
                              data = .x))) %>% 
  mutate(results = map(ols_models, ~tidy(.x, conf.int = TRUE))) 

models_ols_inter_named_b <- models_ols_inter_b %>% 
  pull(ols_models, name = domain)
table_ols_inter_b <- modelsummary(models_ols_inter_named_b, statistic = "({p.value}) {stars}")
table_ols_inter_b
Safe Effective Caring Responsive Well-led Overall
(Intercept) 2.605*** 2.784*** 2.946*** 2.871*** 2.590*** 2.642***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
formNPO 0.099*** 0.089*** 0.064*** 0.129*** 0.112*** 0.117***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
formGOV 0.196*** 0.134*** 0.084* 0.102* 0.179** 0.187***
(0.000) *** (0.001) *** (0.022) * (0.024) * (0.001) ** (0.000) ***
formCIC 0.228** 0.070 0.274*** 0.055 0.203+ 0.243*
(0.010) ** (0.363) (0.000) *** (0.530) (0.066) + (0.014) *
formIND −0.120*** −0.101*** −0.023 −0.034 −0.132*** −0.122***
(0.000) *** (0.000) *** (0.376) (0.281) (0.000) *** (0.000) ***
before_covidTRUE 0.276*** 0.143*** 0.092*** 0.119*** 0.336*** 0.344***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
categoryresidential −0.059*** −0.037*** −0.048*** −0.031*** −0.061*** −0.063***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionEast of England 0.037* 0.039** 0.014 0.034* 0.018 0.032+
(0.017) * (0.002) ** (0.230) (0.017) * (0.342) (0.064) +
regionLondon 0.045** 0.055*** 0.002 −0.002 0.024 0.018
(0.004) ** (0.000) *** (0.843) (0.893) (0.207) (0.306)
regionNorth East 0.132*** 0.103*** 0.074*** 0.089*** 0.168*** 0.149***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionNorth West 0.062*** 0.063*** 0.021+ 0.039** 0.047* 0.046**
(0.000) *** (0.000) *** (0.063) + (0.005) ** (0.012) * (0.006) **
regionSouth East 0.085*** 0.046*** 0.038*** 0.049*** 0.042* 0.047**
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.015) * (0.002) **
regionSouth West 0.109*** 0.091*** 0.082*** 0.088*** 0.128*** 0.125***
(0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) *** (0.000) ***
regionWest Midlands 0.006 0.018 −0.018 −0.002 −0.091*** −0.046**
(0.691) (0.157) (0.119) (0.867) (0.000) *** (0.007) **
regionYorkshire and The Humber −0.016 0.031* 0.035** 0.007 −0.037+ −0.029
(0.314) (0.021) * (0.003) ** (0.623) (0.063) + (0.100)
inheritedTRUE −0.045** −0.042*** −0.037*** −0.038** −0.051** −0.054**
(0.003) ** (0.001) *** (0.001) *** (0.005) ** (0.006) ** (0.001) **
formNPO × before_covidTRUE −0.055* −0.070*** −0.049* −0.086*** −0.056+ −0.088***
(0.020) * (0.001) *** (0.010) * (0.000) *** (0.050) + (0.001) ***
formGOV × before_covidTRUE −0.160** −0.124** −0.071+ −0.067 −0.145* −0.176**
(0.002) ** (0.005) ** (0.076) + (0.181) (0.019) * (0.001) **
formCIC × before_covidTRUE −0.192+ −0.070 −0.248** −0.023 −0.156 −0.215+
(0.059) + (0.424) (0.002) ** (0.818) (0.217) (0.058) +
formIND × before_covidTRUE 0.109** 0.075* 0.017 0.027 0.093* 0.084*
(0.001) ** (0.015) * (0.542) (0.446) (0.024) * (0.023) *
Num.Obs. 12966 12931 12927 12953 12964 12942
R2 0.106 0.043 0.032 0.028 0.108 0.126
R2 Adj. 0.105 0.041 0.030 0.027 0.107 0.125
AIC 13249.0 7796.5 4921.7 10849.4 18390.5 15564.4
BIC 13405.9 7953.3 5078.5 11006.3 18547.3 15721.2
Log.Lik. −6603.505 −3877.241 −2439.846 −5403.713 −9174.233 −7761.182
RMSE 0.40 0.33 0.29 0.37 0.49 0.44

Average Marginal Effects

For curvy fitted model line, the fitted coefficients is less meaningful than a average marginal effect concept.

Dr. Heiss has a very detailed blog post on this: https://www.andrewheiss.com/blog/2022/05/20/marginalia/

marginal_ols_inter <- models_ols_inter %>% 
  mutate(ame = map(ols_models,
                   ~ summary(marginaleffects(.)))) %>%
  mutate(nice_name = paste(domain, "ame"))
marginal_ols_inter_b <- models_ols_inter_b %>% 
  mutate(ame = map(ols_models,
                   ~ summary(marginaleffects(.)))) %>%
  mutate(nice_name = paste(domain, "ame"))
ame_ols_inter_named <- marginal_ols_inter %>% 
  pull(ame, name = nice_name)
ame_ols_inter_named_b <- marginal_ols_inter_b %>% 
  pull(ame, name = nice_name)
table_ame_ols_inter <- tribble(
  ~domain, ~ame_df,
  "safe",  ame_ols_inter_named_b[["Safe ame"]],
  "effective", ame_ols_inter_named_b[["Effective ame"]],
  "caring",  ame_ols_inter_named_b[["Caring ame"]],
  "responsive", ame_ols_inter_named_b[["Responsive ame"]],
  "well-led",  ame_ols_inter_named_b[["Well-led ame"]],
  "overall", ame_ols_inter_named_b[["Overall ame"]]
) %>% 
  unnest(ame_df) %>% 
  filter(term == "form")

table_ame_ols_inter
## # A tibble: 24 × 10
##    domain   type  term  contr…¹ estim…² std.e…³ stati…⁴ p.value conf.low conf.…⁵
##    <chr>    <chr> <chr> <chr>     <dbl>   <dbl>   <dbl>   <dbl>    <dbl>   <dbl>
##  1 safe     resp… form  NPO - …  0.0595 0.0101    5.89  3.78e-9  0.0397   0.0793
##  2 safe     resp… form  GOV - …  0.0816 0.0214    3.82  1.33e-4  0.0398   0.123 
##  3 safe     resp… form  CIC - …  0.0914 0.0445    2.05  3.99e-2  0.00421  0.179 
##  4 safe     resp… form  IND - … -0.0426 0.0149   -2.86  4.20e-3 -0.0718  -0.0134
##  5 effecti… resp… form  NPO - …  0.0347 0.00815   4.26  2.01e-5  0.0188   0.0507
##  6 effecti… resp… form  GOV - …  0.0365 0.0171    2.14  3.26e-2  0.00303  0.0700
##  7 effecti… resp… form  CIC - …  0.0160 0.0360    0.445 6.57e-1 -0.0545   0.0865
##  8 effecti… resp… form  IND - … -0.0426 0.0121   -3.52  4.30e-4 -0.0663  -0.0189
##  9 caring   resp… form  NPO - …  0.0243 0.00729   3.34  8.50e-4  0.0100   0.0386
## 10 caring   resp… form  GOV - …  0.0268 0.0153    1.75  7.99e-2 -0.00319  0.0568
## # … with 14 more rows, and abbreviated variable names ¹​contrast, ²​estimate,
## #   ³​std.error, ⁴​statistic, ⁵​conf.high
table_ame_ols_inter_b <- tribble(
  ~domain, ~ame_df,
  "safe",  ame_ols_inter_named[["Safe ame"]],
  "effective", ame_ols_inter_named[["Effective ame"]],
  "caring",  ame_ols_inter_named[["Caring ame"]],
  "responsive", ame_ols_inter_named[["Responsive ame"]],
  "well-led",  ame_ols_inter_named[["Well-led ame"]],
  "overall", ame_ols_inter_named[["Overall ame"]]
) %>% 
  unnest(ame_df) %>% 
  filter(term == "form")

table_ame_ols_inter
## # A tibble: 24 × 10
##    domain   type  term  contr…¹ estim…² std.e…³ stati…⁴ p.value conf.low conf.…⁵
##    <chr>    <chr> <chr> <chr>     <dbl>   <dbl>   <dbl>   <dbl>    <dbl>   <dbl>
##  1 safe     resp… form  NPO - …  0.0595 0.0101    5.89  3.78e-9  0.0397   0.0793
##  2 safe     resp… form  GOV - …  0.0816 0.0214    3.82  1.33e-4  0.0398   0.123 
##  3 safe     resp… form  CIC - …  0.0914 0.0445    2.05  3.99e-2  0.00421  0.179 
##  4 safe     resp… form  IND - … -0.0426 0.0149   -2.86  4.20e-3 -0.0718  -0.0134
##  5 effecti… resp… form  NPO - …  0.0347 0.00815   4.26  2.01e-5  0.0188   0.0507
##  6 effecti… resp… form  GOV - …  0.0365 0.0171    2.14  3.26e-2  0.00303  0.0700
##  7 effecti… resp… form  CIC - …  0.0160 0.0360    0.445 6.57e-1 -0.0545   0.0865
##  8 effecti… resp… form  IND - … -0.0426 0.0121   -3.52  4.30e-4 -0.0663  -0.0189
##  9 caring   resp… form  NPO - …  0.0243 0.00729   3.34  8.50e-4  0.0100   0.0386
## 10 caring   resp… form  GOV - …  0.0268 0.0153    1.75  7.99e-2 -0.00319  0.0568
## # … with 14 more rows, and abbreviated variable names ¹​contrast, ²​estimate,
## #   ³​std.error, ⁴​statistic, ⁵​conf.high
# commented out for the moment
# table_ols_inter <- modelsummary(ame_ols_inter_named, statistic = "({p.value}) {stars}")
#table_ols_inter

The gap of before vs. after covid rating (before is higher) is larger for for-profit smaller for NPO, GOV, which is consistent with the NYT articles

update: more complex when interpreting non-linear marginal effects http://datacolada.org/57

Running logit models

No Covid Control

# run the model loops with nest() method
models_logit <- socialcare %>% 
  group_by(domain) %>% 
  nest()%>% 
  mutate(logit_models = map(data, 
                          ~glm(rating_higher ~ form + 
                              category + region + inherited , 
                              data = .x,  family = binomial(link = "logit")))) %>% 
  mutate(results = map(logit_models, ~tidy(.x, conf.int = TRUE))) 

models_logit_named <- models_logit %>% 
  pull(logit_models, name = domain)
table_logit <- modelsummary(models_logit_named, statistic = "({p.value}) {stars}")
table_logit

Add Covid control

# run the model loops with nest() method
models_logit <- socialcare %>% 
  group_by(domain) %>% 
  nest()%>% 
  mutate(logit_models = map(data, 
                          ~glm(rating_higher ~ form + 
                              category + region + inherited , 
                              data = .x,  family = binomial(link = "logit")))) %>% 
  mutate(results = map(logit_models, ~tidy(.x, conf.int = TRUE))) 

models_logit_named <- models_logit %>% 
  pull(logit_models, name = domain)
table_logit <- modelsummary(models_logit_named, statistic = "({p.value}) {stars}")
table_logit
# run the model loops with nest() method
models_logit_covid <- socialcare %>% 
  group_by(domain) %>% 
  nest()%>% 
  mutate(logit_models = map(data, 
                          ~glm(rating_higher ~ form + after_covid +
                              category + region + inherited , 
                              data = .x,  family = binomial(link = "logit")))) %>% 
  mutate(results = map(logit_models, ~tidy(.x, conf.int = TRUE))) 

models_logit_covid_named <- models_logit_covid %>% 
  pull(logit_models, name = domain)
table_logit_covid <- modelsummary(models_logit_covid_named, statistic = "({p.value}) {stars}")
table_logit_covid
# run the model loops with nest() method
models_logit_inter <- socialcare %>% 
  group_by(domain) %>% 
  nest()%>% 
  mutate(logit_models = map(data, 
                          ~glm(rating_higher ~ form * after_covid +
                              category + region + inherited , 
                              data = .x,  family = binomial(link = "logit")))) %>% 
  mutate(results = map(logit_models, ~tidy(.x, conf.int = TRUE))) 

models_logit_inter_named <- models_logit_inter %>% 
  pull(logit_models, name = domain)
table_logit_inter <- modelsummary(models_logit_inter_named, statistic = "({p.value}) {stars}")
table_logit_inter