Data Preparation

First, import the sampled and coded data set

#import the curated social care data file
#socialcare_raw <- read_csv("data/sample_new_cleaned.csv")
socialcare_new <-read_rds(here("data", "derived_data", "rating_clean2025.rds"))

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_new %>% 
  # recode legal form types to be more readable / easier to present
  mutate(location_id = factor(location_id),
         inherited = ifelse(inherited == "Y", TRUE, FALSE),
         rating = recode(rating, 
                         "Insufficient evidence to rate" = "NA",
                         "Requires improvement" = "Req improv"),
 #        publication_date = dmy(publication_date)
         region = location_region,
         time_pass = interval(provider_hsca_start_date, publication_date) / years(1),
         age = scale(time_pass)
         ) %>% 
  # set the order of the values in the factors 
  mutate(form = fct_relevel(form, "FPO"),
         form2 = fct_relevel(form, "CIC"),
         form3 = fct_relevel(form, "GOV"),
         
         # 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),
         during_covid = ifelse(year >= 2020 & year <2023, TRUE, FALSE),
         COVID =   publication_date >= as.Date("2020-02-10") & 
                   publication_date <= as.Date("2022-12-31")) %>%

  # 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, TRUE, FALSE),
         Year = factor(year)) %>% 
  filter(!is.na(form))

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

## # A tibble: 6 × 59
##   location_id   location_ods_code location_name care_home location_type  
##   <fct>         <chr>             <chr>         <chr>     <chr>          
## 1 1-10000302982 VNH4N             Henley House  Y         Social Care Org
## 2 1-10000302982 VNH4N             Henley House  Y         Social Care Org
## 3 1-10000302982 VNH4N             Henley House  Y         Social Care Org
## 4 1-10000302982 VNH4N             Henley House  Y         Social Care Org
## 5 1-10000302982 VNH4N             Henley House  Y         Social Care Org
## 6 1-10000302982 VNH4N             Henley House  Y         Social Care Org
## # ℹ 54 more variables: primary_cat <chr>, location_street_address <chr>,
## #   location_address_line_2 <chr>, location_city <chr>,
## #   location_post_code <chr>, location_local_authority <chr>,
## #   location_region <chr>, location_nhs_region <chr>,
## #   location_onspd_ccg_code <chr>, location_onspd_ccg <chr>,
## #   location_commissioning_ccg_code <chr>,
## #   location_commissioning_ccg_name <chr>, service_population_group <chr>, …

Descriptive data

# Define service user bands to include as controls
service_bands <- c(
  "service_user_band_sensory_impairment",
  "service_user_band_mental_health",
  "service_user_band_learning_disabilities_or_autistic_spectrum_disorder",
  "service_user_band_physical_disability",
  "service_user_band_dementia",
  "service_user_band_whole_population",
  "service_user_band_younger_adults",
  "service_user_band_older_people"
)

Summerize at different Levels

socialcare %>% 
  summarise(observation = n(),
            location = n_distinct(location_id),
            provider = n_distinct(provider_name))

## # A tibble: 1 × 3
##   observation location provider
##         <int>    <int>    <int>
## 1      151199    23498    13413

form_level <- socialcare %>% 
  group_by(form) %>% 
  summarise(observation = n(),
            location = n_distinct(location_id),
            provider = n_distinct(provider_name))
form_level

## # A tibble: 5 × 4
##   form  observation location provider
##   <fct>       <int>    <int>    <int>
## 1 FPO        123027    19011    11481
## 2 CIC           543       88       63
## 3 GOV          4192      679      133
## 4 IND          7576     1201     1025
## 5 NPO         15861     2519      711

form_pct <- form_level %>%
  mutate(across(-form, ~ round(.x / sum(.x) * 100, 2)))

form_pct

## # A tibble: 5 × 4
##   form  observation location provider
##   <fct>       <dbl>    <dbl>    <dbl>
## 1 FPO         81.4     80.9     85.6 
## 2 CIC          0.36     0.37     0.47
## 3 GOV          2.77     2.89     0.99
## 4 IND          5.01     5.11     7.64
## 5 NPO         10.5     10.7      5.3

table1 <- datasummary((form + category + inherited + COVID +  region + rating + 1) ~ 1,
                      data = socialcare, fmt = 0) 
table1

tinytable_x9x28slijxp5d0azfzxk

		All
form	FPO	123027
	CIC	543
	GOV	4192
	IND	7576
	NPO	15861
category	community	61393
	residential	89788
inherited	FALSE	142775
	TRUE	8424
COVID	FALSE	110570
	TRUE	39627
region	East Midlands	14375
	East of England	17523
	London	18354
	North East	6863
	North West	18142
	South East	27265
	South West	17837
	Unspecified	139
	West Midlands	16223
	Yorkshire and The Humber	14478
rating	Inadequate	1090
	Req improv	18960
	Good	125073
	Outstanding	4975
	All	151199

table_1b <- datasummary(
  (service_user_band_sensory_impairment +
   service_user_band_mental_health +
   service_user_band_learning_disabilities_or_autistic_spectrum_disorder +
   service_user_band_physical_disability +
   service_user_band_dementia +
   service_user_band_whole_population +
   service_user_band_younger_adults +
   service_user_band_older_people) ~ 1,
  data = socialcare, fmt = 0
)

table_1b

tinytable_q8u4gj3zthoks9suufzr

		All
service_user_band_sensory_impairment	FALSE	90474
	TRUE	60719
service_user_band_mental_health	FALSE	88345
	TRUE	62848
service_user_band_learning_disabilities_or_autistic_spectrum_disorder	FALSE	78831
	TRUE	72362
service_user_band_physical_disability	FALSE	65117
	TRUE	86076
service_user_band_dementia	FALSE	56233
	TRUE	94960
service_user_band_whole_population	FALSE	148329
	TRUE	2864
service_user_band_younger_adults	FALSE	47771
	TRUE	103422
service_user_band_older_people	FALSE	29775
	TRUE	121418

# calculating ratios
N <- nrow(socialcare)
N

## [1] 151199

# Get a *data.frame* from datasummary
table1_df <- datasummary(
  (form + category + inherited + COVID + region + rating + 1) ~ 1,
  data = socialcare,
  fmt = 0,
  output = "data.frame"   # <- key: return a data.frame, not a tinytable
)

# Ensure counts are numeric (strip commas/formatting if present), then add % column
table1_df <- table1_df %>%
  mutate(
    All = as.numeric(gsub("[^0-9.-]", "", All)),
    group_ratio = round(All / N * 100, 2)
  )

table1_df

##                                          All group_ratio
## 1       form                      FPO 123027       81.37
## 2                                 CIC    543        0.36
## 3                                 GOV   4192        2.77
## 4                                 IND   7576        5.01
## 5                                 NPO  15861       10.49
## 6   category                community  61393       40.60
## 7                         residential  89788       59.38
## 8  inherited                    FALSE 142775       94.43
## 9                                TRUE   8424        5.57
## 10     COVID                    FALSE 110570       73.13
## 11                               TRUE  39627       26.21
## 12    region            East Midlands  14375        9.51
## 13                    East of England  17523       11.59
## 14                             London  18354       12.14
## 15                         North East   6863        4.54
## 16                         North West  18142       12.00
## 17                         South East  27265       18.03
## 18                         South West  17837       11.80
## 19                        Unspecified    139        0.09
## 20                      West Midlands  16223       10.73
## 21           Yorkshire and The Humber  14478        9.58
## 22    rating               Inadequate   1090        0.72
## 23                         Req improv  18960       12.54
## 24                               Good 125073       82.72
## 25                        Outstanding   4975        3.29
## 26                                All 151199      100.00

table_ap1 <- datasummary((category + inherited + COVID +  region + rating + 1) ~ form,
                      data = socialcare, fmt = 0)
table_ap1

tinytable_874cykk0992j6ieobsxb

		FPO	CIC	GOV	IND	NPO
category	community	51851	441	2217	1707	5177
	residential	71158	102	1975	5869	10684
inherited	FALSE	115798	531	3943	7468	15035
	TRUE	7229	12	249	108	826
COVID	FALSE	88397	370	3231	5918	12654
	TRUE	33814	161	925	1640	3087
region	East Midlands	12147	24	418	702	1084
	East of England	14849	114	394	741	1425
	London	15286	72	415	820	1761
	North East	5426	42	300	214	881
	North West	14117	132	825	989	2079
	South East	22115	24	624	1371	3131
	South West	13783	33	268	1420	2333
	Unspecified	121	0	6	0	12
	West Midlands	13465	42	336	748	1632
	Yorkshire and The Humber	11718	60	606	571	1523
rating	Inadequate	901	0	17	92	80
	Req improv	16093	55	304	1086	1422
	Good	101197	433	3656	6202	13585
	Outstanding	3927	43	179	178	648
	All	123027	543	4192	7576	15861

table_ap1b <- datasummary(
  (service_user_band_sensory_impairment +
   service_user_band_mental_health +
   service_user_band_learning_disabilities_or_autistic_spectrum_disorder +
   service_user_band_physical_disability +
   service_user_band_dementia +
   service_user_band_whole_population +
   service_user_band_younger_adults +
   service_user_band_older_people) ~ form,
  data = socialcare, fmt = 0
)

table_ap1b

tinytable_6yhxyl8rga80dn38si09

		FPO	CIC	GOV	IND	NPO
service_user_band_sensory_impairment	FALSE	71795	258	2072	5800	10549
	TRUE	51232	285	2120	1776	5306
service_user_band_mental_health	FALSE	68642	240	2512	5126	11825
	TRUE	54385	303	1680	2450	4030
service_user_band_learning_disabilities_or_autistic_spectrum_disorder	FALSE	66240	90	1704	4954	5843
	TRUE	56787	453	2488	2622	10012
service_user_band_physical_disability	FALSE	49905	150	1468	4855	8739
	TRUE	73122	393	2724	2721	7116
service_user_band_dementia	FALSE	40083	210	1674	3726	10540
	TRUE	82944	333	2518	3850	5315
service_user_band_whole_population	FALSE	120463	519	4162	7510	15675
	TRUE	2564	24	30	66	180
service_user_band_younger_adults	FALSE	34676	84	1226	4530	7255
	TRUE	88351	459	2966	3046	8600
service_user_band_older_people	FALSE	21193	120	808	1693	5961
	TRUE	101834	423	3384	5883	9894

Regression analysis

logit models - simple regression

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

tinytable_jo3yrwq0aegw8leic0x1

	overall	safe	effective	caring	well-led	responsive
Inadequate\|Req improv	-4.564	-4.405	-5.724	-6.250	-4.186	-6.042
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Req improv\|Good	-1.498	-1.387	-2.236	-3.241	-1.157	-2.603
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Good\|Outstanding	3.183	5.665	4.295	3.089	3.202	2.980
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
formCIC	0.631	0.817	0.266	0.961	0.581	0.496
	(0.024) *	(0.018) *	(0.475)	(0.004) **	(0.027) *	(0.129)
formGOV	0.622	0.698	0.534	0.339	0.573	0.387
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(0.013) *	(<0.001) ***	(0.001) **
formIND	-0.152	-0.107	-0.260	-0.023	-0.250	-0.158
	(0.023) *	(0.125)	(0.003) **	(0.830)	(<0.001) ***	(0.077) +
formNPO	0.376	0.438	0.338	0.298	0.337	0.544
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Num.Obs.	25414	25388	24627	24771	25153	24745
AIC	34004.5	28127.3	19581.9	17201.0	38017.8	22379.1
BIC	34061.5	28184.3	19638.7	17257.8	38074.7	22435.9
RMSE	2.25	2.13	2.13	2.20	2.26	2.24

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

tinytable_xwhmq8v9ci959wk815gq

	overall	safe	effective	caring	well-led	responsive
Inadequate\|Req improv	0.010	0.012	0.003	0.002	0.015	0.002
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Req improv\|Good	0.224	0.250	0.107	0.039	0.315	0.074
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Good\|Outstanding	24.119	288.639	73.336	21.944	24.574	19.690
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
formCIC	1.880	2.263	1.305	2.614	1.789	1.642
	(0.024) *	(0.018) *	(0.475)	(0.004) **	(0.027) *	(0.129)
formGOV	1.862	2.010	1.706	1.404	1.774	1.473
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(0.013) *	(<0.001) ***	(0.001) **
formIND	0.859	0.899	0.771	0.977	0.779	0.854
	(0.023) *	(0.125)	(0.003) **	(0.830)	(<0.001) ***	(0.077) +
formNPO	1.456	1.549	1.402	1.348	1.401	1.723
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Num.Obs.	25414	25388	24627	24771	25153	24745
AIC	34004.5	28127.3	19581.9	17201.0	38017.8	22379.1
BIC	34061.5	28184.3	19638.7	17257.8	38074.7	22435.9
RMSE	2.25	2.13	2.13	2.20	2.26	2.24

model with controls + adding service band, FPO as the reference group

# Define service user bands to include as controls
service_bands <- c(
  "service_user_band_sensory_impairment",
  "service_user_band_mental_health",
  "service_user_band_learning_disabilities_or_autistic_spectrum_disorder",
  "service_user_band_physical_disability",
  "service_user_band_dementia",
  "service_user_band_whole_population",
  "service_user_band_younger_adults",
  "service_user_band_older_people"
)

# ------------------------------------------------------------------------------
# Fixed Effects Models with Service Band Controls
# ------------------------------------------------------------------------------

# Overall domain
model_order_overall_covid_band <- clm(
  as.formula(paste(
    "rating ~ form + COVID + category + region + inherited +",
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Overall"),
  link = "logit"
)

# Safe domain
model_order_safe_covid_band <- clm(
  as.formula(paste(
    "rating ~ form + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Safe"),
  link = "logit"
)

# Effective domain
model_order_effective_covid_band <- clm(
  as.formula(paste(
    "rating ~ form + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Effective"),
  link = "logit"
)

# Caring domain
model_order_caring_covid_band <- clm(
  as.formula(paste(
    "rating ~ form + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Caring"),
  link = "logit"
)

# Well-led domain
model_order_well_led_covid_band <- clm(
  as.formula(paste(
    "rating ~ form + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Well-led"),
  link = "logit"
)

# Responsive domain
model_order_responsive_covid_band <- clm(
  as.formula(paste(
    "rating ~ form + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Responsive"),
  link = "logit"
)

Log-odds table

# ------------------------------------------------------------------------------
# Generate Summary Table (Hiding Service Bands and Region)
# ------------------------------------------------------------------------------

ordinal_models_covid_band_summary <-
  modelsummary(
    list(
      "Overall" = model_order_overall_covid_band,
      "Safe" = model_order_safe_covid_band,
      "Effective" = model_order_effective_covid_band,
      "Caring" = model_order_caring_covid_band,
      "Well-led" = model_order_well_led_covid_band,
      "Responsive" = model_order_responsive_covid_band
    ),
    coef_omit = "region|service_user_band_.*",  # Hide region and all service bands
    exponentiate = F,  # Odds ratios for interpretation
    statistic = "({p.value}) {stars}",
    title = "Ordinal Regression Models with Service Band Controls"
  )

# Display table
ordinal_models_covid_band_summary

tinytable_31l0kwwcvr4exg10v1o1

Ordinal Regression Models with Service Band Controls
	Overall	Safe	Effective	Caring	Well-led	Responsive
Inadequate\|Req improv	-5.277	-5.120	-6.267	-7.018	-4.775	-6.628
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Req improv\|Good	-2.189	-2.075	-2.768	-4.001	-1.719	-3.181
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Good\|Outstanding	2.592	5.101	3.878	2.510	2.723	2.517
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
formCIC	0.466	0.614	0.063	0.773	0.438	0.424
	(0.096) +	(0.079) +	(0.865)	(0.020) *	(0.096) +	(0.196)
formGOV	0.525	0.607	0.431	0.213	0.482	0.293
	(<0.001) ***	(<0.001) ***	(0.002) **	(0.122)	(<0.001) ***	(0.015) *
formIND	-0.179	-0.134	-0.288	-0.082	-0.271	-0.259
	(0.009) **	(0.062) +	(0.001) **	(0.449)	(<0.001) ***	(0.004) **
formNPO	0.218	0.256	0.197	0.135	0.199	0.307
	(<0.001) ***	(<0.001) ***	(0.009) **	(0.082) +	(<0.001) ***	(<0.001) ***
COVIDTRUE	-0.235	-0.163	-0.243	-0.249	-0.220	-0.248
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
categoryresidential	-0.446	-0.501	-0.553	-0.642	-0.329	-0.247
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
inheritedTRUE	0.006	-0.018	-0.027	-0.100	0.016	0.003
	(0.926)	(0.794)	(0.761)	(0.318)	(0.802)	(0.972)
Num.Obs.	25410	25384	24623	24767	25149	24741
AIC	33499.0	27594.2	19332.2	16863.0	37485.2	22085.7
BIC	33718.8	27814.0	19551.2	17082.2	37704.8	22304.8
RMSE	2.25	2.13	2.13	2.20	2.25	2.24

Exponentiated odds ratio model

# ------------------------------------------------------------------------------
# Generate Summary Table (Hiding Service Bands and Region)
# ------------------------------------------------------------------------------

ordinal_models_covid_band_exp <-
  modelsummary(
    list(
      "Overall" = model_order_overall_covid_band,
      "Safe" = model_order_safe_covid_band,
      "Effective" = model_order_effective_covid_band,
      "Caring" = model_order_caring_covid_band,
      "Well-led" = model_order_well_led_covid_band,
      "Responsive" = model_order_responsive_covid_band
    ),
    coef_omit = "region|service_user_band_.*",  # Hide region and all service bands
    exponentiate = T,  # Odds ratios for interpretation
    statistic = "({p.value}) {stars}",
    title = "Ordinal Regression Models with Service Band Controls"
  )

# Display table
ordinal_models_covid_band_exp

tinytable_vqhmp23eo0mkx5xf8cxv

Ordinal Regression Models with Service Band Controls
	Overall	Safe	Effective	Caring	Well-led	Responsive
Inadequate\|Req improv	0.005	0.006	0.002	0.001	0.008	0.001
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Req improv\|Good	0.112	0.126	0.063	0.018	0.179	0.042
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Good\|Outstanding	13.361	164.108	48.323	12.311	15.219	12.390
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
formCIC	1.593	1.848	1.065	2.167	1.549	1.529
	(0.096) +	(0.079) +	(0.865)	(0.020) *	(0.096) +	(0.196)
formGOV	1.690	1.834	1.539	1.237	1.619	1.341
	(<0.001) ***	(<0.001) ***	(0.002) **	(0.122)	(<0.001) ***	(0.015) *
formIND	0.836	0.875	0.750	0.921	0.762	0.772
	(0.009) **	(0.062) +	(0.001) **	(0.449)	(<0.001) ***	(0.004) **
formNPO	1.243	1.292	1.218	1.144	1.220	1.360
	(<0.001) ***	(<0.001) ***	(0.009) **	(0.082) +	(<0.001) ***	(<0.001) ***
COVIDTRUE	0.790	0.849	0.784	0.780	0.802	0.780
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
categoryresidential	0.640	0.606	0.575	0.526	0.719	0.781
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
inheritedTRUE	1.006	0.982	0.974	0.905	1.016	1.003
	(0.926)	(0.794)	(0.761)	(0.318)	(0.802)	(0.972)
Num.Obs.	25410	25384	24623	24767	25149	24741
AIC	33499.0	27594.2	19332.2	16863.0	37485.2	22085.7
BIC	33718.8	27814.0	19551.2	17082.2	37704.8	22304.8
RMSE	2.25	2.13	2.13	2.20	2.25	2.24

Run residential and community based seperately

## models for residential

# ------------------------------------------------------------------------------
# Fixed Effects Models with Service Band Controls
# ------------------------------------------------------------------------------

# Overall domain
model_order_overall_covid_bandr <- clm(
  as.formula(paste(
    "rating ~ form + COVID + region + inherited +",
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Overall" & category == "residential"),
  link = "logit"
)

# Safe domain
model_order_safe_covid_bandr <- clm(
  as.formula(paste(
    "rating ~ form + COVID + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Safe" & category == "residential"),
  link = "logit"
)

# Effective domain
model_order_effective_covid_bandr <- clm(
  as.formula(paste(
    "rating ~ form + COVID + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Effective" & category == "residential"),
  link = "logit"
)

# Caring domain
model_order_caring_covid_bandr <- clm(
  as.formula(paste(
    "rating ~ form + COVID + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Caring" & category == "residential"),
  link = "logit"
)

# Well-led domain
model_order_well_led_covid_bandr <- clm(
  as.formula(paste(
    "rating ~ form + COVID + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Well-led" & category == "residential"),
  link = "logit"
)

# Responsive domain
model_order_responsive_covid_bandr <- clm(
  as.formula(paste(
    "rating ~ form + COVID + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Responsive"& category == "residential"),
  link = "logit"
)

# ------------------------------------------------------------------------------
# Generate Summary Table (Hiding Service Bands and Region)
# ------------------------------------------------------------------------------

ordinal_models_covid_bandr_exp <-
  modelsummary(
    list(
      "Overall" = model_order_overall_covid_bandr,
      "Safe" = model_order_safe_covid_bandr,
      "Effective" = model_order_effective_covid_bandr,
      "Caring" = model_order_caring_covid_bandr,
      "Well-led" = model_order_well_led_covid_bandr,
      "Responsive" = model_order_responsive_covid_bandr
    ),
    coef_omit = "region|service_user_band_.*",  # Hide region and all service bands
    exponentiate = T,  # Odds ratios for interpretation
    statistic = "({p.value}) {stars}",
    title = "Ordinal Regression Models with Service Band Controls"
  )

# Display table
ordinal_models_covid_bandr_exp

tinytable_pa6hzsv5vqi0kikarr74

Ordinal Regression Models with Service Band Controls
	Overall	Safe	Effective	Caring	Well-led	Responsive
Inadequate\|Req improv	0.009	0.010	0.004	0.002	0.013	0.002
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Req improv\|Good	0.164	0.186	0.106	0.040	0.247	0.059
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Good\|Outstanding	17.328	233.342	61.760	23.439	20.368	13.491
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
formCIC	1.612	4.503	0.474	3.017	2.594	0.898
	(0.416)	(0.093) +	(0.187)	(0.176)	(0.120)	(0.875)
formGOV	1.409	1.391	1.366	1.367	1.272	1.025
	(0.012) *	(0.025) *	(0.088) +	(0.119)	(0.059) +	(0.882)
formIND	0.784	0.826	0.698	0.893	0.726	0.721
	(0.001) **	(0.016) *	(<0.001) ***	(0.351)	(<0.001) ***	(<0.001) ***
formNPO	1.239	1.262	1.265	1.207	1.249	1.259
	(<0.001) ***	(<0.001) ***	(0.007) **	(0.044) *	(<0.001) ***	(0.003) **
COVIDTRUE	0.800	0.832	0.780	0.882	0.798	0.783
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(0.092) +	(<0.001) ***	(<0.001) ***
inheritedTRUE	1.039	0.998	0.992	0.938	1.084	1.009
	(0.642)	(0.981)	(0.939)	(0.610)	(0.305)	(0.935)
Num.Obs.	15217	15197	14681	14760	15036	14730
AIC	21244.7	17828.5	12836.4	10591.8	23129.7	14354.5
BIC	21443.1	18026.8	13033.8	10789.4	23327.8	14552.0
RMSE	2.25	2.13	2.14	2.19	2.24	2.26

## models for community based services 

# ------------------------------------------------------------------------------
# Fixed Effects Models with Service Band Controls
# ------------------------------------------------------------------------------

# Overall domain
model_order_overall_covid_bandc <- clm(
  as.formula(paste(
    "rating ~ form + COVID + region + inherited +",
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Overall" & category == "community"),
  link = "logit"
)

# Safe domain
model_order_safe_covid_bandc <- clm(
  as.formula(paste(
    "rating ~ form + COVID + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Safe" & category == "community"),
  link = "logit"
)

# Effective domain
model_order_effective_covid_bandc <- clm(
  as.formula(paste(
    "rating ~ form + COVID + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Effective" & category == "community"),
  link = "logit"
)

# Caring domain
model_order_caring_covid_bandc <- clm(
  as.formula(paste(
    "rating ~ form + COVID + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Caring" & category == "community"),
  link = "logit"
)

# Well-led domain
model_order_well_led_covid_bandc <- clm(
  as.formula(paste(
    "rating ~ form + COVID + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Well-led" & category == "community"),
  link = "logit"
)

# Responsive domain
model_order_responsive_covid_bandc <- clm(
  as.formula(paste(
    "rating ~ form + COVID + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Responsive"& category == "community"),
  link = "logit"
)

# ------------------------------------------------------------------------------
# Generate Summary Table (Hiding Service Bands and Region)
# ------------------------------------------------------------------------------

ordinal_models_covid_bandc_exp <-
  modelsummary(
    list(
      "Overall" = model_order_overall_covid_bandc,
      "Safe" = model_order_safe_covid_bandc,
      "Effective" = model_order_effective_covid_bandc,
      "Caring" = model_order_caring_covid_bandc,
      "Well-led" = model_order_well_led_covid_bandc,
      "Responsive" = model_order_responsive_covid_bandc
    ),
    coef_omit = "region|service_user_band_.*",  # Hide region and all service bands
    exponentiate = T,  # Odds ratios for interpretation
    statistic = "({p.value}) {stars}",
    title = "Ordinal Regression Models with Service Band Controls"
  )

# Display table
ordinal_models_covid_bandc_exp

tinytable_yc52y9np4ygjlt4hw0wv

Ordinal Regression Models with Service Band Controls
	Overall	Safe	Effective	Caring	Well-led	Responsive
Inadequate\|Req improv	0.004	0.005	0.001	0.000	0.006	0.001
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Req improv\|Good	0.119	0.157	0.065	0.012	0.167	0.032
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Good\|Outstanding	17.451	229.337	81.995	11.131	15.175	16.297
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
formCIC	1.665	1.476	1.765	2.128	1.406	1.956
	(0.119)	(0.307)	(0.259)	(0.044) *	(0.246)	(0.082) +
formGOV	2.158	2.945	1.925	1.158	2.147	1.911
	(<0.001) ***	(<0.001) ***	(0.004) **	(0.448)	(<0.001) ***	(<0.001) ***
formIND	1.030	1.047	1.074	0.950	0.871	1.055
	(0.849)	(0.787)	(0.754)	(0.822)	(0.318)	(0.796)
formNPO	1.236	1.346	1.065	1.008	1.151	1.594
	(0.035) *	(0.012) *	(0.667)	(0.953)	(0.120)	(<0.001) ***
COVIDTRUE	0.782	0.895	0.796	0.658	0.816	0.778
	(<0.001) ***	(0.054) +	(0.003) **	(<0.001) ***	(<0.001) ***	(<0.001) ***
inheritedTRUE	0.960	0.983	0.957	0.852	0.918	0.999
	(0.718)	(0.890)	(0.783)	(0.345)	(0.397)	(0.996)
Num.Obs.	10193	10187	9942	10007	10113	10011
AIC	12196.9	9721.7	6442.3	6243.1	14319.6	7634.2
BIC	12384.8	9909.6	6629.6	6430.6	14507.4	7821.7
RMSE	2.25	2.12	2.11	2.21	2.27	2.20

Average Marginal Effect

#testing with overall model
overall_ame <- model_order_overall_covid_band |>
  avg_slopes(variables = "form")  
overall_ame

## 
##        Group Term              Contrast Estimate Std. Error     z Pr(>|z|)    S
##  Inadequate  form mean(CIC) - mean(FPO) -0.00374   0.001789 -2.09  0.03673  4.8
##  Inadequate  form mean(GOV) - mean(FPO) -0.00410   0.000655 -6.26  < 0.001 31.3
##  Inadequate  form mean(IND) - mean(FPO)  0.00195   0.000819  2.38  0.01715  5.9
##  Inadequate  form mean(NPO) - mean(FPO) -0.00196   0.000460 -4.26  < 0.001 15.6
##  Req improv  form mean(CIC) - mean(FPO) -0.05431   0.027791 -1.95  0.05066  4.3
##  Req improv  form mean(GOV) - mean(FPO) -0.06002   0.009575 -6.27  < 0.001 31.4
##  Req improv  form mean(IND) - mean(FPO)  0.02536   0.010188  2.49  0.01281  6.3
##  Req improv  form mean(NPO) - mean(FPO) -0.02747   0.006431 -4.27  < 0.001 15.7
##  Good        form mean(CIC) - mean(FPO)  0.03567   0.013242  2.69  0.00708  7.1
##  Good        form mean(GOV) - mean(FPO)  0.03820   0.004300  8.88  < 0.001 60.4
##  Good        form mean(IND) - mean(FPO) -0.02092   0.008716 -2.40  0.01639  5.9
##  Good        form mean(NPO) - mean(FPO)  0.02010   0.004419  4.55  < 0.001 17.5
##  Outstanding form mean(CIC) - mean(FPO)  0.02239   0.016368  1.37  0.17141  2.5
##  Outstanding form mean(GOV) - mean(FPO)  0.02592   0.006052  4.28  < 0.001 15.7
##  Outstanding form mean(IND) - mean(FPO) -0.00639   0.002294 -2.79  0.00533  7.6
##  Outstanding form mean(NPO) - mean(FPO)  0.00933   0.002484  3.75  < 0.001 12.5
##      2.5 %    97.5 %
##  -0.007245 -0.000230
##  -0.005382 -0.002816
##   0.000347  0.003557
##  -0.002862 -0.001059
##  -0.108782  0.000155
##  -0.078786 -0.041253
##   0.005390  0.045328
##  -0.040073 -0.014865
##   0.009710  0.061620
##   0.029774  0.046630
##  -0.038002 -0.003837
##   0.011442  0.028765
##  -0.009694  0.054466
##   0.014055  0.037778
##  -0.010889 -0.001895
##   0.004457  0.014194
## 
## Columns: term, group, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted 
## Type:  prob

# Compute AMEs for "form" only across domains
overall_ame <- avg_slopes(model_order_overall_covid_band, variables = "form")
safe_ame <- avg_slopes(model_order_safe_covid_band, variables = "form")
effective_ame <- avg_slopes(model_order_effective_covid_band, variables = "form")
caring_ame <- avg_slopes(model_order_caring_covid_band, variables = "form")
well_led_ame <- avg_slopes(model_order_well_led_covid_band, variables = "form")
responsive_ame <- avg_slopes(model_order_responsive_covid_band, variables = "form")

# Combine into one data frame with domain labels
ame_results <- bind_rows(
  mutate(overall_ame, domain = "Overall"),
  mutate(safe_ame, domain = "Safe"),
  mutate(effective_ame, domain = "Effective"),
  mutate(caring_ame, domain = "Caring"),
  mutate(well_led_ame, domain = "Well-led"),
  mutate(responsive_ame, domain = "Responsive")
)

ame_results

## 
##        Group Term              Contrast  Estimate Std. Error      z Pr(>|z|)
##  Inadequate  form mean(CIC) - mean(FPO) -0.003738   0.001789 -2.089  0.03673
##  Inadequate  form mean(GOV) - mean(FPO) -0.004099   0.000655 -6.261  < 0.001
##  Inadequate  form mean(IND) - mean(FPO)  0.001952   0.000819  2.384  0.01715
##  Inadequate  form mean(NPO) - mean(FPO) -0.001960   0.000460 -4.261  < 0.001
##  Req improv  form mean(CIC) - mean(FPO) -0.054313   0.027791 -1.954  0.05066
##     S    CI low   CI high
##   4.8 -0.007245 -0.000230
##  31.3 -0.005382 -0.002816
##   5.9  0.000347  0.003557
##  15.6 -0.002862 -0.001059
##   4.3 -0.108782  0.000155
## --- 86 rows omitted. See ?print.marginaleffects --- 
##  Good        form mean(NPO) - mean(FPO)  0.000127   0.000922  0.138  0.89038
##  Outstanding form mean(CIC) - mean(FPO)  0.024194   0.022278  1.086  0.27747
##  Outstanding form mean(GOV) - mean(FPO)  0.015775   0.007269  2.170  0.03000
##  Outstanding form mean(IND) - mean(FPO) -0.010910   0.003454 -3.159  0.00158
##  Outstanding form mean(NPO) - mean(FPO)  0.016611   0.003977  4.177  < 0.001
##     S    CI low   CI high
##   0.2 -0.001680  0.001934
##   1.8 -0.019469  0.067857
##   5.1  0.001527  0.030023
##   9.3 -0.017679 -0.004141
##  15.0  0.008816  0.024405
## Columns: term, group, contrast, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted, domain

Full model with controls, CIC as the reference group

# ------------------------------------------------------------------------------
# Fixed Effects Models with Service Band Controls
# ------------------------------------------------------------------------------

# Overall domain
model_order_overall_covid_band2 <- clm(
  as.formula(paste(
    "rating ~ form2 + COVID + category + region + inherited +",
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Overall"),
  link = "logit"
)

# Safe domain
model_order_safe_covid_band2 <- clm(
  as.formula(paste(
    "rating ~ form2 + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Safe"),
  link = "logit"
)

# Effective domain
model_order_effective_covid_band2 <- clm(
  as.formula(paste(
    "rating ~ form2 + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Effective"),
  link = "logit"
)

# Caring domain
model_order_caring_covid_band2 <- clm(
  as.formula(paste(
    "rating ~ form2 + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Caring"),
  link = "logit"
)

# Well-led domain
model_order_well_led_covid_band2 <- clm(
  as.formula(paste(
    "rating ~ form2 + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Well-led"),
  link = "logit"
)

# Responsive domain
model_order_responsive_covid_band2 <- clm(
  as.formula(paste(
    "rating ~ form2 + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Responsive"),
  link = "logit"
)

Log-odds table

# ------------------------------------------------------------------------------
# Generate Summary Table (Hiding Service Bands and Region)
# ------------------------------------------------------------------------------

ordinal_models_covid_band_summary2 <-
  modelsummary(
    list(
      "Overall" = model_order_overall_covid_band2,
      "Safe" = model_order_safe_covid_band2,
      "Effective" = model_order_effective_covid_band2,
      "Caring" = model_order_caring_covid_band2,
      "Well-led" = model_order_well_led_covid_band2,
      "Responsive" = model_order_responsive_covid_band2
    ),
    coef_omit = "region|service_user_band_.*",  # Hide region and all service bands
    exponentiate = F,  # Odds ratios for interpretation
    statistic = "({p.value}) {stars}",
    title = "Ordinal Regression Models with Service Band Controls"
  )

# Display table
ordinal_models_covid_band_summary2

tinytable_owij0375jdd7scfm3l5g

Ordinal Regression Models with Service Band Controls
	Overall	Safe	Effective	Caring	Well-led	Responsive
Inadequate\|Req improv	-5.743	-5.734	-6.331	-7.791	-5.213	-7.052
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Req improv\|Good	-2.655	-2.689	-2.831	-4.774	-2.157	-3.606
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Good\|Outstanding	2.127	4.487	3.815	1.737	2.285	2.093
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
form2FPO	-0.466	-0.614	-0.063	-0.773	-0.438	-0.424
	(0.096) +	(0.079) +	(0.865)	(0.020) *	(0.096) +	(0.196)
form2GOV	0.059	-0.007	0.368	-0.560	0.044	-0.131
	(0.842)	(0.984)	(0.352)	(0.117)	(0.874)	(0.707)
form2IND	-0.645	-0.748	-0.352	-0.855	-0.709	-0.684
	(0.025) *	(0.035) *	(0.356)	(0.014) *	(0.009) **	(0.044) *
form2NPO	-0.248	-0.358	0.134	-0.639	-0.239	-0.117
	(0.382)	(0.311)	(0.723)	(0.060) +	(0.370)	(0.725)
COVIDTRUE	-0.235	-0.163	-0.243	-0.249	-0.220	-0.248
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
categoryresidential	-0.446	-0.501	-0.553	-0.642	-0.329	-0.247
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
inheritedTRUE	0.006	-0.018	-0.027	-0.100	0.016	0.003
	(0.926)	(0.794)	(0.761)	(0.318)	(0.802)	(0.972)
Num.Obs.	25410	25384	24623	24767	25149	24741
AIC	33499.0	27594.2	19332.2	16863.0	37485.2	22085.7
BIC	33718.8	27814.0	19551.2	17082.2	37704.8	22304.8
RMSE	2.25	2.13	2.13	2.20	2.25	2.24

Exponentiated odds ratio model

# ------------------------------------------------------------------------------
# Generate Summary Table (Hiding Service Bands and Region)
# ------------------------------------------------------------------------------

ordinal_models_covid_band_exp2 <-
  modelsummary(
    list(
      "Overall" = model_order_overall_covid_band2,
      "Safe" = model_order_safe_covid_band2,
      "Effective" = model_order_effective_covid_band2,
      "Caring" = model_order_caring_covid_band2,
      "Well-led" = model_order_well_led_covid_band2,
      "Responsive" = model_order_responsive_covid_band2
    ),
    coef_omit = "region|service_user_band_.*",  # Hide region and all service bands
    exponentiate = T,  # Odds ratios for interpretation
    statistic = "({p.value}) {stars}",
    title = "Ordinal Regression Models with Service Band Controls"
  )

# Display table
ordinal_models_covid_band_exp2

tinytable_bxev039fq15d904djs8u

Ordinal Regression Models with Service Band Controls
	Overall	Safe	Effective	Caring	Well-led	Responsive
Inadequate\|Req improv	0.003	0.003	0.002	0.000	0.005	0.001
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Req improv\|Good	0.070	0.068	0.059	0.008	0.116	0.027
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Good\|Outstanding	8.386	88.820	45.360	5.682	9.823	8.106
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
form2FPO	0.628	0.541	0.939	0.462	0.645	0.654
	(0.096) +	(0.079) +	(0.865)	(0.020) *	(0.096) +	(0.196)
form2GOV	1.061	0.993	1.444	0.571	1.045	0.877
	(0.842)	(0.984)	(0.352)	(0.117)	(0.874)	(0.707)
form2IND	0.525	0.473	0.704	0.425	0.492	0.505
	(0.025) *	(0.035) *	(0.356)	(0.014) *	(0.009) **	(0.044) *
form2NPO	0.780	0.699	1.143	0.528	0.787	0.889
	(0.382)	(0.311)	(0.723)	(0.060) +	(0.370)	(0.725)
COVIDTRUE	0.790	0.849	0.784	0.780	0.802	0.780
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
categoryresidential	0.640	0.606	0.575	0.526	0.719	0.781
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
inheritedTRUE	1.006	0.982	0.974	0.905	1.016	1.003
	(0.926)	(0.794)	(0.761)	(0.318)	(0.802)	(0.972)
Num.Obs.	25410	25384	24623	24767	25149	24741
AIC	33499.0	27594.2	19332.2	16863.0	37485.2	22085.7
BIC	33718.8	27814.0	19551.2	17082.2	37704.8	22304.8
RMSE	2.25	2.13	2.13	2.20	2.25	2.24

Full model with controls, GOV as the reference group

# ------------------------------------------------------------------------------
# Fixed Effects Models with Service Band Controls
# ------------------------------------------------------------------------------

# Overall domain
model_order_overall_covid_band3 <- clm(
  as.formula(paste(
    "rating ~ form3 + COVID + category + region + inherited +",
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Overall"),
  link = "logit"
)

# Safe domain
model_order_safe_covid_band3 <- clm(
  as.formula(paste(
    "rating ~ form3 + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Safe"),
  link = "logit"
)

# Effective domain
model_order_effective_covid_band3 <- clm(
  as.formula(paste(
    "rating ~ form3 + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Effective"),
  link = "logit"
)

# Caring domain
model_order_caring_covid_band3 <- clm(
  as.formula(paste(
    "rating ~ form3 + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Caring"),
  link = "logit"
)

# Well-led domain
model_order_well_led_covid_band3 <- clm(
  as.formula(paste(
    "rating ~ form3 + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Well-led"),
  link = "logit"
)

# Responsive domain
model_order_responsive_covid_band3 <- clm(
  as.formula(paste(
    "rating ~ form3 + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Responsive"),
  link = "logit"
)

Log-odds table

# ------------------------------------------------------------------------------
# Generate Summary Table (Hiding Service Bands and Region)
# ------------------------------------------------------------------------------

ordinal_models_covid_band_summary3 <-
  modelsummary(
    list(
      "Overall" = model_order_overall_covid_band3,
      "Safe" = model_order_safe_covid_band3,
      "Effective" = model_order_effective_covid_band3,
      "Caring" = model_order_caring_covid_band3,
      "Well-led" = model_order_well_led_covid_band3,
      "Responsive" = model_order_responsive_covid_band3
    ),
    coef_omit = "region|service_user_band_.*",  # Hide region and all service bands
    exponentiate = F,  # Odds ratios for interpretation
    statistic = "({p.value}) {stars}",
    title = "Ordinal Regression Models with Service Band Controls"
  )

# Display table
ordinal_models_covid_band_summary3

tinytable_q2qj2tg8hpo9hgo0kvhf

Ordinal Regression Models with Service Band Controls
	Overall	Safe	Effective	Caring	Well-led	Responsive
Inadequate\|Req improv	-5.801	-5.727	-6.698	-7.230	-5.257	-6.921
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Req improv\|Good	-2.714	-2.682	-3.199	-4.214	-2.201	-3.475
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Good\|Outstanding	2.068	4.494	3.447	2.298	2.241	2.223
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
form3FPO	-0.525	-0.607	-0.431	-0.213	-0.482	-0.293
	(<0.001) ***	(<0.001) ***	(0.002) **	(0.122)	(<0.001) ***	(0.015) *
form3CIC	-0.059	0.007	-0.368	0.560	-0.044	0.131
	(0.842)	(0.984)	(0.352)	(0.117)	(0.874)	(0.707)
form3IND	-0.703	-0.741	-0.719	-0.295	-0.753	-0.553
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(0.085) +	(<0.001) ***	(<0.001) ***
form3NPO	-0.307	-0.351	-0.234	-0.078	-0.283	0.014
	(0.005) **	(0.008) **	(0.132)	(0.608)	(0.006) **	(0.918)
COVIDTRUE	-0.235	-0.163	-0.243	-0.249	-0.220	-0.248
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
categoryresidential	-0.446	-0.501	-0.553	-0.642	-0.329	-0.247
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
inheritedTRUE	0.006	-0.018	-0.027	-0.100	0.016	0.003
	(0.926)	(0.794)	(0.761)	(0.318)	(0.802)	(0.972)
Num.Obs.	25410	25384	24623	24767	25149	24741
AIC	33499.0	27594.2	19332.2	16863.0	37485.2	22085.7
BIC	33718.8	27814.0	19551.2	17082.2	37704.8	22304.8
RMSE	2.25	2.13	2.13	2.20	2.25	2.24

Exponentiated odds ratio model

# ------------------------------------------------------------------------------
# Generate Summary Table (Hiding Service Bands and Region)
# ------------------------------------------------------------------------------

ordinal_models_covid_band_exp3 <-
  modelsummary(
    list(
      "Overall" = model_order_overall_covid_band3,
      "Safe" = model_order_safe_covid_band3,
      "Effective" = model_order_effective_covid_band3,
      "Caring" = model_order_caring_covid_band3,
      "Well-led" = model_order_well_led_covid_band3,
      "Responsive" = model_order_responsive_covid_band3
    ),
    coef_omit = "region|service_user_band_.*",  # Hide region and all service bands
    exponentiate = T,  # Odds ratios for interpretation
    statistic = "({p.value}) {stars}",
    title = "Ordinal Regression Models with Service Band Controls"
  )

# Display table
ordinal_models_covid_band_exp3

tinytable_tluspflkc9loqzchaels

Ordinal Regression Models with Service Band Controls
	Overall	Safe	Effective	Caring	Well-led	Responsive
Inadequate\|Req improv	0.003	0.003	0.001	0.001	0.005	0.001
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Req improv\|Good	0.066	0.068	0.041	0.015	0.111	0.031
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
Good\|Outstanding	7.908	89.470	31.407	9.952	9.399	9.239
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
form3FPO	0.592	0.545	0.650	0.808	0.618	0.746
	(<0.001) ***	(<0.001) ***	(0.002) **	(0.122)	(<0.001) ***	(0.015) *
form3CIC	0.943	1.007	0.692	1.751	0.957	1.140
	(0.842)	(0.984)	(0.352)	(0.117)	(0.874)	(0.707)
form3IND	0.495	0.477	0.487	0.745	0.471	0.575
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(0.085) +	(<0.001) ***	(<0.001) ***
form3NPO	0.736	0.704	0.792	0.925	0.753	1.014
	(0.005) **	(0.008) **	(0.132)	(0.608)	(0.006) **	(0.918)
COVIDTRUE	0.790	0.849	0.784	0.780	0.802	0.780
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
categoryresidential	0.640	0.606	0.575	0.526	0.719	0.781
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
inheritedTRUE	1.006	0.982	0.974	0.905	1.016	1.003
	(0.926)	(0.794)	(0.761)	(0.318)	(0.802)	(0.972)
Num.Obs.	25410	25384	24623	24767	25149	24741
AIC	33499.0	27594.2	19332.2	16863.0	37485.2	22085.7
BIC	33718.8	27814.0	19551.2	17082.2	37704.8	22304.8
RMSE	2.25	2.13	2.13	2.20	2.25	2.24

Linear model for appendix

# Overall domain
model_linear_overall_covid_band <- lm(
  as.formula(paste(
    "rating_num ~ form  + COVID + category + region + inherited +",
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Overall")
)

# Safe domain
model_linear_safe_covid_band <- lm(
  as.formula(paste(
    "rating_num ~ form  + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Safe")
)

# Effective domain
model_linear_effective_covid_band <- lm(
  as.formula(paste(
    "rating_num ~ form + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Effective")
)

# Caring domain
model_linear_caring_covid_band <- lm(
  as.formula(paste(
    "rating_num ~ form + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Caring")
)

# Well-led domain
model_linear_well_led_covid_band <- lm(
  as.formula(paste(
    "rating_num ~ form  + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Well-led")
)

# Responsive domain
model_linear_responsive_covid_band <- lm(
  as.formula(paste(
    "rating_num ~ form  + COVID + category + region + inherited +", 
    paste(service_bands, collapse = " + ")
  )),
  data = filter(socialcare, domain == "Responsive")
)

# ------------------------------------------------------------------------------
# Generate Summary Table (Hiding Service Bands and Region) - Linear Models
# ------------------------------------------------------------------------------

linear_models_covid_band_summary <-
  modelsummary(
    list(
      "Overall" = model_linear_overall_covid_band,
      "Safe" = model_linear_safe_covid_band,
      "Effective" = model_linear_effective_covid_band,
      "Caring" = model_linear_caring_covid_band,
      "Well-led" = model_linear_well_led_covid_band,
      "Responsive" = model_linear_responsive_covid_band
    ),
    coef_omit = "region|service_user_band_.*",  # Hide region and all service bands
    exponentiate = FALSE,
    statistic = "({p.value}) {stars}",
    title = "Linear Models (DV = Rating_num) with Service Band Controls"
  )

# Display table
linear_models_covid_band_summary

tinytable_m9vcisqedrrdt5bn1hv8

Linear Models (DV = Rating_num) with Service Band Controls
	Overall	Safe	Effective	Caring	Well-led	Responsive
(Intercept)	2.970	2.893	2.955	3.060	2.904	3.034
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
formCIC	0.084	0.082	0.008	0.071	0.101	0.049
	(0.093) +	(0.073) +	(0.827)	(0.024) *	(0.067) +	(0.192)
formGOV	0.091	0.079	0.039	0.020	0.100	0.035
	(<0.001) ***	(<0.001) ***	(0.003) **	(0.086) +	(<0.001) ***	(0.011) *
formIND	-0.039	-0.029	-0.035	-0.004	-0.069	-0.027
	(0.005) **	(0.024) *	(<0.001) ***	(0.635)	(<0.001) ***	(0.009) **
formNPO	0.040	0.036	0.016	0.010	0.044	0.034
	(<0.001) ***	(<0.001) ***	(0.026) *	(0.101)	(<0.001) ***	(<0.001) ***
COVIDTRUE	-0.036	-0.017	-0.021	-0.017	-0.043	-0.025
	(<0.001) ***	(0.004) **	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
categoryresidential	-0.086	-0.082	-0.054	-0.053	-0.076	-0.028
	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***	(<0.001) ***
inheritedTRUE	-0.003	-0.006	-0.004	-0.007	0.000	-0.002
	(0.834)	(0.596)	(0.690)	(0.393)	(0.977)	(0.825)
Num.Obs.	25410	25384	24623	24767	25149	24741
R2	0.024	0.025	0.013	0.016	0.025	0.016
R2 Adj.	0.023	0.024	0.012	0.015	0.025	0.015
AIC	33952.1	29297.1	15411.7	9625.5	38328.4	18317.1
BIC	34163.9	29508.8	15622.6	9836.6	38539.8	18528.1
Log.Lik.	-16950.067	-14622.560	-7679.873	-4786.756	-19138.190	-9132.535
RMSE	0.47	0.43	0.33	0.29	0.52	0.35

Social Care CQC data analysis full data set 2025 - for MS revision and tables

Meng Ye

09/01/2025

Data Preparation

Descriptive data

Summerize at different Levels

Regression analysis

logit models - simple regression

model with controls + adding service band, FPO as the reference group

Log-odds table

Exponentiated odds ratio model

Run residential and community based seperately

Average Marginal Effect

Full model with controls, CIC as the reference group

Log-odds table

Exponentiated odds ratio model

Full model with controls, GOV as the reference group

Log-odds table

Exponentiated odds ratio model

Linear model for appendix