Multivariable Linear Model for Aged Care Star Ratings

Analysis
MOA-Benchmarking
Author

Filip Reierson

Published

May 26, 2023

I have previously performed descriptive analysis of the star ratings for Q2 FY22-23. In this article, I used the same dataset, and a linear model to assess the importance of various factors in explaining the overall star rating. The data is available on the GEN Aged Care Data website (Australian Government, Department of Health and Aged Care 2023). I have also shared the R code here for those who would like to understand the details or follow along. The results of the analysis were shared on LinkedIn in a slide deck which can be seen in Figure 1.

Figure 1: The findings were shared in a slide deck.

Model fitting

Code 
library(dplyr)
library(tidyr)
library(tibble)
library(readxl)
stars_raw <-
  read_excel(
    'star-ratings-extract.xlsx',
    sheet = 2
  ) |>
  select(
    resi_experience = `Residents' Experience rating`,
    compliance = `Compliance rating`,
    staffing = `Staffing rating`,
    quality_measures = `Quality Measures rating`,
    mmm_code = `MMM Code`,
    service_size = Size,
    service = `Service Name`,
    provider = `Provider Name`,
    purpose = Purpose
  )

stars <- stars_raw |>
  mutate(
    rating = resi_experience * .33 +
      compliance * .3 +
      staffing * .22 +
      quality_measures * .15,
    mmm_group = factor(
      case_when(
        mmm_code == 'MM1' ~ 'Metropolitan',
        mmm_code == 'MM2' ~ 'Regional centres',
        T ~ 'Rural and remote'
      ),
      levels = c('Metropolitan', 'Regional centres', 'Rural and remote')
    ),
    service_size = factor(service_size, levels = c('Small', 'Medium', 'Large'))
  ) |>
  group_by(provider) |>
  mutate(provider_services = n(),
         provider_size = factor(
           case_when(
             provider_services < 5 ~ 'Small',
             provider_services < 20 ~ 'Medium',
             provider_services >= 20 ~ 'Large',
             T ~ NA_character_
           ),
           levels = c('Small', 'Medium', 'Large')
         )) |>
  ungroup() |>
  drop_na()

The regression model was fitted on the unrounded overall star ratings, without the compliance rule, which otherwise sometimes caps the overall star rating. This makes the residuals for the model follow a distribution closer to a normal distribution. The ratings themselves are also approximately normal, as Figure 2 shows.

Code 
x <- seq(1,5,length.out=1000)
norm_density <- dnorm(x,
                      mean = mean(stars$rating),
                      sd = sd(stars$rating))
empirical_density <- density(stars$rating)
hist(stars$rating, n=20, prob = TRUE, ylim = c(0, max(empirical_density$y)), main = "", xlab = 'Star rating')
lines(empirical_density, lwd = 2, col = 'cornflowerblue')
lines(x, norm_density, lwd = 2, col='red')
Figure 2: Histogram of unrounded star ratings. Normal density, based on the sample mean and standard deviation, shown in red. Empirical density shown in blue.

Next we fit the linear model.

Code 
mdl <- lm(rating ~ purpose + service_size +
            provider_size + mmm_group, data = stars)
mdl |> summary()

Call:
lm(formula = rating ~ purpose + service_size + provider_size + 
    mmm_group, data = stars)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.73366 -0.30386  0.01325  0.31298  1.52157 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                3.464489   0.025236 137.282   <2e-16 ***
purposeGovernment          0.567326   0.040127  14.138   <2e-16 ***
purposeNot for profit      0.029961   0.020844   1.437   0.1507    
service_sizeMedium        -0.190595   0.022034  -8.650   <2e-16 ***
service_sizeLarge         -0.277905   0.025525 -10.888   <2e-16 ***
provider_sizeMedium        0.003260   0.022631   0.144   0.8855    
provider_sizeLarge         0.028753   0.022671   1.268   0.2048    
mmm_groupRegional centres -0.004366   0.034559  -0.126   0.8995    
mmm_groupRural and remote -0.038159   0.022839  -1.671   0.0949 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.457 on 2491 degrees of freedom
Multiple R-squared:  0.1737,    Adjusted R-squared:  0.1711 
F-statistic: 65.47 on 8 and 2491 DF,  p-value: < 2.2e-16
Code 
mdl |> anova()
Analysis of Variance Table

Response: rating
                Df Sum Sq Mean Sq  F value Pr(>F)    
purpose          2  80.29  40.146 192.2546 <2e-16 ***
service_size     2  28.09  14.044  67.2579 <2e-16 ***
provider_size    2   0.39   0.195   0.9316 0.3941    
mmm_group        2   0.59   0.297   1.4228 0.2412    
Residuals     2491 520.16   0.209                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Diagnostics

To confirm that we can interpret the results directly, we should check the diagnostic plots in Figure 3. Diagnostics appear good, especially given the high number of observations. Residuals are approximately normal based on the Q-Q plot. The residuals versus fitted values shows a mostly even spread along the horizontal, indicating we have homoscedasticity.

Code 
par(mfrow=c(2,2))
plot(mdl)
Figure 3: Diagnostics of the linear model treating unrounded star ratings as the dependent variable.

Sensitivity analysis

Code 
library(gt)
mdl2 <- lm((rating-staffing*0.22)/(1-0.22) ~ purpose + service_size +
            provider_size + mmm_group, data = stars)
gov_v_fp <- round(coef(mdl)[['purposeGovernment']], 2)
gov_v_fp_staffing_excl <- round(coef(mdl2)[['purposeGovernment']], 2)

Based on previous explorations of the data we know that government homes generally had a higher star rating for staffing than non-government homes. To determine if this is the reason that government homes performed better overall, we calculated a new star rating which excluded the staffing star rating, then fitted the same model.

Based on the above model output, government homes have an average star rating that is 0.57 higher than for profit homes, controlling for the other factors in the model. However, if we exclude the staffing star rating we instead get an effect size of 0.27. We conclude that government homes generally received better ratings than for profit homes, but the magnitude of the difference is partly driven by staffing.

Table of coefficients

In Table 1 we show the model coefficients, their confidence intervals, and associated p-values.

Code 
library(broom)
library(gt)
tidy(mdl) |>
  mutate(
    term = case_when(
      term == '(Intercept)' ~ 'Intercept',
      term == 'purposeGovernment' ~ 'Government',
      term == 'purposeNot for profit' ~ 'Not for profit',
      term == 'service_sizeMedium' ~ '61-100 residents',
      term == 'service_sizeLarge' ~ '101+ residents',
      term == 'provider_sizeMedium' ~ '5-19 services',
      term == 'provider_sizeLarge' ~ '20+ services',
      term == 'mmm_groupRegional centres' ~ 'Regional centres',
      term == 'mmm_groupRural and remote' ~ 'Rural and remote'
    )
  ) |>
  mutate(LowerCI = estimate - 1.96*std.error,
         UpperCI = estimate + 1.96*std.error) |>
  select(term, estimate, LowerCI, UpperCI, p.value) |>
  rename(`p-value` = p.value) |>
  gt() |>
  cols_label('LowerCI'='95%l',
             'UpperCI'='95%u'
             ) |>
  fmt_number(c(estimate, LowerCI, UpperCI)) |>
  fmt(`p-value`, fns = scales::label_pvalue()) |>
  tab_footnote('Star rating relative to for profit homes.', locations = cells_body(1,2:3)) |>
  tab_footnote(paste0('If we excluded staffing, the estimate would be ', round(gov_v_fp_staffing_excl,2), 
               ' instead.'), locations = cells_body(1,2)) |>
  tab_footnote('Star rating relative to a service with 1-60 residents.', locations = cells_body(1,4:5))|>
  tab_footnote('Star rating relative to a provider with 1-4 services.', locations = cells_body(1,6:7)) |>
  tab_footnote('Star rating relative to a service in a metropolitan area.', locations = cells_body(1,8:9)) |>
  tab_footnote('Predicted star rating for a service that is for profit, metropolitan, has 1-60 residents, and the provider has 1-4 services.', locations = cells_body(1,1))
Table 1: Table of coefficients for a multivariable linear model regressing unrounded star ratings on properties of the service.
term estimate 95%l 95%u p-value
Intercept1 3.46 3.42 3.51 <0.001
Government2,3 0.57 0.49 0.65 <0.001
Not for profit2 0.03 −0.01 0.07 0.151
61-100 residents4 −0.19 −0.23 −0.15 <0.001
101+ residents4 −0.28 −0.33 −0.23 <0.001
5-19 services5 0.00 −0.04 0.05 0.885
20+ services5 0.03 −0.02 0.07 0.205
Regional centres6 0.00 −0.07 0.06 0.899
Rural and remote6 −0.04 −0.08 0.01 0.095
1 Predicted star rating for a service that is for profit, metropolitan, has 1-60 residents, and the provider has 1-4 services.
2 Star rating relative to for profit homes.
3 If we excluded staffing, the estimate would be 0.27 instead.
4 Star rating relative to a service with 1-60 residents.
5 Star rating relative to a provider with 1-4 services.
6 Star rating relative to a service in a metropolitan area.

References

Australian Government, Department of Health and Aged Care. 2023. “Star Ratings Quarterly Data Extract – 2022-23 Financial Year – Quarter 2.”