On the flight back from South Africa, a headline from McKnights popped into my email:
Immediately, my BS detector started buzzing. As Carl Sagan said, “Extraordinary claims require extraordinary evidence”. So, I went to the LTCCC’s site, and looked closer at the claim:
The claim still seemed too good to be true. So, I decided to replicate their findings.
First I loaded up the data:
import polars as pl
from plotnine import *
from plotnine.data import *
pbj_data = 'PBJ_Daily_Non-Nurse_Staffing_CY_*.csv'
target_columns = ['PROVNUM', 'STATE', 'CY_Qtr', 'WorkDate', 'MDScensus', 'Hrs_Admin'] # the columns of interest
convert_columns_to_integer = pl.col('MDScensus').cast(pl.Int16)
convert_columns_to_floats = pl.col('Hrs_Admin').cast(pl.Float64)
convert_columns_to_date = pl.col('WorkDate').str.strptime(pl.Date, fmt='%Y%m%d')
df_pbj = (pl
.scan_csv(pbj_data, infer_schema_length=0)
.select(target_columns)
.with_columns([convert_columns_to_integer, convert_columns_to_floats, convert_columns_to_date])
.collect()
)
df_pbj.head()
Ok, so let’s write the code to replicate their findings:
calculate_total_resident_days = pl.col('MDScensus').sum().alias('total_resident_days')
calculate_total_admin_hours = pl.col('Hrs_Admin').sum().alias('total_admin_hours')
calculate_avg_hrs = pl.col('Hrs_Admin').mean().round(4).alias('avg_hrs')
create_admin_hrs_per_resident_column = (pl.col('total_admin_hours') / pl.col('total_resident_days')).alias('admin_hrs_per_resident')
create_avg_hrs_pct_change_column = pl.col('avg_hrs').pct_change().round(4).alias('avg_hrs_pct_change')
create_admin_hrs_per_resident_pct_change_column = pl.col('admin_hrs_per_resident').pct_change().round(4).alias('admin_hrs_per_resident_pct_change')
df_ltccc = (df_pbj
.groupby('CY_Qtr')
.agg([calculate_avg_hrs, calculate_total_resident_days, calculate_total_admin_hours])
.with_columns([create_admin_hrs_per_resident_column])
.sort(by=['CY_Qtr'])
.with_columns([create_avg_hrs_pct_change_column, create_admin_hrs_per_resident_pct_change_column])
.drop(['total_resident_days', 'total_admin_hours'])
)
df_ltccc
We are able to replicate the LTCCC’s findings. We can see the change from 8.44 to 6.17 hrs, from 2019-Q3 to 2022-Q3. I was also able to replicate the -20% time per resident claim (the last column). Fair enough.
But something is strange about the “Hrs_Admin” column, which represents the hours worked. Let’s take a look at the distribution:
(ggplot(df_pbj, aes(x='Hrs_Admin', color='CY_Qtr')) +
geom_density() +
theme(dpi=100)
)
Ok. This is strange. There are data that are beyond 24 hrs! How can someone work more than 24 hours in a day?! There are values up to just under 500 hours. We might have a data problem here. Let’s break out the data by quarter and hours that are over 24 hours:
df_pbj.groupby('CY_Qtr')\
.agg([pl.col('Hrs_Admin').count().alias('total_count'),
pl.col('Hrs_Admin').where(pl.col('Hrs_Admin') > 24).count().alias('count_over_24hrs'),
(pl.col('Hrs_Admin').where(pl.col('Hrs_Admin') > 24).count() / pl.col('Hrs_Admin').count()).round(3).alias('prop_over_24hrs')])\
.sort(by='CY_Qtr')
Ok, so we see that 6.7% of the data in 2019-Q3 is over 24 hours, whereas it’s only 1.9% in 2022-Q3. This suggests that data quality improved over time, which is what we would expect. It also suggests that the LTCCC did not account for this data quality issue, which is extraordinarily naive and reckless, considering the claims that it is making.
Let’s look at the data split at 24 hours or under, and over 24 hours:
(ggplot(df_pbj.filter(pl.col('Hrs_Admin') <=24), aes(x='Hrs_Admin', color='CY_Qtr')) +
annotate(geom_vline, xintercept=24, color='purple', size=1, linetype='dashed') +
annotate(geom_text, label='24-hour mark', x=20, y=.6, ha='right', size=10, color='purple') +
annotate(geom_segment, x=23.5, xend=20, y=.6, yend=.6, arrow=arrow(ends='first', type='closed'), color='purple') +
geom_density() +
theme(dpi=100)
)
The distribution of the Hrs_Admin at or under 24 hours looks very similar. We see the peaking at 8 hours, which is a normal shift, and the subpeaks at 4 hours (half-shift), 6 hours, and 16 hours (double-shift). The data is very similar across quarters. Now let’s look at greater than 24 hours across quarters:
(ggplot(df_pbj.filter(pl.col('Hrs_Admin') > 24).select(['CY_Qtr', 'Hrs_Admin']), aes(x='Hrs_Admin', fill='CY_Qtr')) +
geom_histogram(bins=50, alpha=.5) +
annotate(geom_vline, xintercept=24, color='purple', size=1, linetype='dashed') +
annotate(geom_text, label='24-hour mark', x=100, y=35_000, ha='left', size=10, color='purple') +
annotate(geom_segment, x=30, xend=100, y=35_000, yend=35_000, arrow=arrow(ends='first', type='closed'), color='purple') +
annotate(geom_segment, x=80, xend=200, y=10_000, yend=15_000, arrow=arrow(ends='first', type='closed'), color='red') +
annotate(geom_text, label='2019-Q3 has more than 3.5x \n the number of values \n than 2022-Q3!', x=210, y=15_000, ha='left', size=10, color='red') +
theme(dpi=100)
)
Yeah, so definitely have a problem here. The 2019-Q3 data has a lot more data indicating that more than 24 hours were worked. And as we can see, the 2019-Q3 errata are more numerous and more right-tailed. This probably explains most of the change in administrator hours that LTCCC claims.
There is a formal statistical way to deal with errata like this in the data (via a process called Bayesian imputation and regularization) but it is complicated and beyond the scope of this article. But, let’s apply a crude common-sense adjustment to see its effect. For all the data greater than 24 hours, we are going to set the value to a full shift (8 hours) and recalculate the LTCCC’s numbers, using the very same calculation methodology.
pbj_data = 'PBJ_Daily_Non-Nurse_Staffing_CY_*.csv'
target_columns = ['PROVNUM', 'STATE', 'CY_Qtr', 'WorkDate', 'MDScensus', 'Hrs_Admin'] # the columns of interest
convert_columns_to_integer = pl.col('MDScensus').cast(pl.Int16)
convert_columns_to_floats = pl.col('Hrs_Admin').cast(pl.Float64)
convert_columns_to_date = pl.col('WorkDate').str.strptime(pl.Date, fmt='%Y%m%d')
df_pbj = (pl
.scan_csv(pbj_data, infer_schema_length=0)
.select(target_columns)
.with_columns([convert_columns_to_integer, convert_columns_to_floats, convert_columns_to_date])
.with_columns([pl.when(pl.col('Hrs_Admin') > 24).then(pl.lit(8)).otherwise(pl.col('Hrs_Admin')).alias('Hrs_Admin')]) <<--- # making the crude adjustment for the errata!
.collect()
)
calculate_total_resident_days = pl.col('MDScensus').sum().alias('total_resident_days')
calculate_total_admin_hours = pl.col('Hrs_Admin').sum().alias('total_admin_hours')
calculate_avg_hrs = pl.col('Hrs_Admin').mean().round(4).alias('avg_hrs')
create_admin_hrs_per_resident_column = (pl.col('total_admin_hours') / pl.col('total_resident_days')).alias('admin_hrs_per_resident')
create_avg_hrs_pct_change_column = pl.col('avg_hrs').pct_change().round(4).alias('avg_hrs_pct_change')
create_admin_hrs_per_resident_pct_change_column = pl.col('admin_hrs_per_resident').pct_change().round(4).alias('admin_hrs_per_resident_pct_change')
df_ltccc = (df_pbj
.groupby('CY_Qtr')
.agg([calculate_avg_hrs, calculate_total_resident_days, calculate_total_admin_hours])
.with_columns([create_admin_hrs_per_resident_column])
.sort(by=['CY_Qtr'])
.with_columns([create_avg_hrs_pct_change_column, create_admin_hrs_per_resident_pct_change_column])
.drop(['total_resident_days', 'total_admin_hours'])
)
Once that adjustment is made, the average hours go from 5.52 in 2019-Q3 to 5.39 in 2022-Q3. This is a -2.4%, less than -10x the average hours delta, without the adjustment.
And, with this adjustment, Administrator hours per resident stay day actually INCREASE by 6.76%, not decrease by 20%!
Simply put, the LTCCC’s claim is bogus. And it is bogus because they did not even attempt to account for obvious errata in the data. Granted, my adjustment is crude but it clearly shows, via negativa, that the LTCCC’s claims are not credible!
Before making claims that impugn an industry and its workers, I suggest a modicum of humility and basic care, and common sense in working with data. The LTCCC has an agenda, and I respect that, but, as the kids say, “clout-chasing” undercuts that agenda and besmirches people who are actually working hard to care for residents.
