This vignette shows how to use the functions in QualityMeasure to analyze quality measure performance and reliability. In the first half, we demonstrate how to use functions for measures without risk adjustment. In the second half, we show how to incorporate risk adjustment into the analyses.

Example #1: Measures without risk adjustment

Simulate data

We’ll simulate some data to use for the example analyses in this section.

### Data parameters
n.entity = 100  # number of accountable entities
n.obs = 50      # average number of observations per accountable entity
mu = .2         # marginal probability of the outcome
r = .7          # median reliability

### Simulate data
example_df_1 <- simulateData(n.entity = n.entity, n.obs = n.obs, mu = mu, r = r)

Provider profiling analyses

The calcPerformance() function can be used to obtain some basic provider profiling results, including confidence intervals for entity-level rates.

profiling.results <- calcPerformance(df = example_df_1)
perf.results <- profiling.results$perf.results

Several control parameters are available through the controlPerf() parameter in calcPerformance(), including:

min.n which filters out entities with sample sizes below the specified amount; default is 2.
alpha which specifies the statistical significance level; default is 0.05.

The estimated marginal probability of the outcome is 0.23. Recall that the true marginal probability according to our data-generation process is 0.2.

Using the perf.results dataframe, you can obtain summary statistics of the performance distribution. We see that the outcome rate varies from 0.02 to 0.542 with a median entity-level rate of 0.209.

Performance summary statistics across entities
Method	Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
Unadjusted	0.02	0.153	0.209	0.224	0.289	0.542

Number of observations per entity

plotN(perf.results$n)

Unadjusted rates across entities

plotPerformance(df = perf.results)

Categorization of entities according to whether rates are lower, higher, or no different from the overall average can also be found in the perf.results output from calcPerformance() function.

Categorization of entities based on their outcome rates
Category	Number of entities
Higher than average	15
Lower than average	12
No different than average	73

Plot of ORs comparing each entity with the average

Another way to assess whether entities have measure performance that differs from the national average is to examine the predicted random intercept values. These values likely will show fewer outliers, given that estimates are shrunken towards the mean.

Entities highlighted in red have estimated ORs that are statistically different from 1.0 at alpha = 0.05. There are 19 entities with outcome odds that are significantly different from the average entity performance.

plotPerformance(df = perf.results, plot.type = 'OR')

Beta-binomial-based reliability estimates

You can calculate reliability estimates based on a Beta-Binomial model by using the calcBetaBin() function.

tutorial_BB_results <- calcBetaBin(df = example_df_1)

The estimated alpha is 4.783 and beta is 16.494, which you can find in tutorial_BB_results.

Summary statistics of the reliability distribution:

summary(tutorial_BB_results$est.BB)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.6349  0.6837  0.6972  0.6995  0.7135  0.7694

Alternatively, if you have data that are aggregated by entity, you can calculate Beta-Binomial estimates by using

df.agg <- data.frame(n = aggregate(y ~ entity, data = example_df_1, length)$y,
                     x = aggregate(y ~ entity, data = example_df_1, sum)$y)

tutorial_BB_agg_results <- calcBetaBin(df = df.agg, df.aggregate = T, n = 'n', x = 'x')`

where n is the column of sample sizes for each entity and x is the column of counts of observations meeting measure numerator criteria for each entity.

Again, the estimated alpha is 4.783 and beta is 16.494, and the distribution of reliability estimates is summarzied by:

#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.6349  0.6837  0.6972  0.6995  0.7135  0.7694

These estimates match our results from analyzing the data before aggregation.

Calculate reliability from all methods

If you would like to estimate reliability using the split-sample method, Beta-Binomial method, and hierarchical logistic regression method at once, you can use the calcReliability() function.

Again, minimum sample size per entity can be specified through the min.n parameter in controlPerf(). In addition, several control parameters are available through the controlRel() function:

n.resamples specifies the number of resamples to use for the calculation of the permutation split-sample reliability estimate; default is 100.
n.cores specifies the number of cores to use for parallel processing when calculating the split-sample estimates; default is 2.
SSR.method specifies whether the permutation method (i.e., resampling without replacement) or bootstrap method (i.e., resampling with replacement) should be used for split-sample reliability estimates. Default is the permutation method.

tutorial_reliability_results_1 <- calcReliability(df = example_df_1, ctrRel = controlRel(n.resamples = 100))

Reliability estimates
Method	Reliability	Min Reliability	Max Reliability
Permutation split-sample	0.697	NA	NA
Hierarchical logit regression	0.696	0.634	0.768
Beta-Binomial	0.697	0.635	0.769

plotReliability(tutorial_reliability_results_1)

Example #2: Measures with risk adjustment

In this next part of the tutorial, we will work with an example measure where risk adjustment is required.

Simulate data

We can use the built-in function simulateData() to simulate data from a hierarchical logistic regression model with covariates for risk adjustment. The simulated data will include a continuous covariate x1 which is sampled from a standard Normal distribution.

### Data parameters 
n.entity = 100  # number of accountable entities
n.obs = 50 # average number of patients/cases per accountable entity
mu = .2 # marginal probability of the outcome
r = .7 # reliability for entity with an median number of patients
beta1 = log(1.5) # parameter for risk adjustment model---coefficient for x1 which is simulated from a standard Normal


### Simulate data 
example_df_2 <- simulateData(n.entity = n.entity, n.obs = n.obs, mu = mu, r = r, beta1 = beta1)

Simulated data with a covariate
entity	z	x1	lp	p	y
1	0.480334	0.8538955	-0.5597355	0.3636087	1
1	0.480334	0.4187967	-0.7361529	0.3238460	0
1	0.480334	0.3399565	-0.7681199	0.3168860	0
1	0.480334	0.5964251	-0.6641308	0.3398123	1
1	0.480334	1.8714180	-0.1471657	0.4632748	0
1	0.480334	0.6028704	-0.6615175	0.3403988	0
1	0.480334	-0.7661688	-1.2166151	0.2285327	0
1	0.480334	-0.6203265	-1.1574811	0.2391253	0
1	0.480334	0.7901903	-0.5855658	0.3576529	1
1	0.480334	-0.2419765	-1.0040734	0.2681413	0

Fit risk adjustment model to the data

To incorporate risk adjustment, we specify a model to use for risk adjustment:

model = 'y ~ x1 + (1 | entity)'

This is a model that adjusts for x1 and includes a random intercept for entity.

Estimates of model parameters

model.perf <- model_performance(df = example_df_2, model = model)
plotEstimates(model.perf)

The estimated value of the regression coefficient, after exponentiating is, 1.44 with a 95% confidence interval ranging from 1.34 to 1.55.

Discrimination

The estimated c-statistic is 0.69, and below is a plot of the distributions of predicted probabilities separately for observations with and without the outcome occurring.

plotPredictedDistribution(model.perf)

Calibration

Below is a plot of the observed rate of the outcome for quantiles of the predicted probabilities from the model. The number of quantiles can be specified through the quantiles argument of the plotCalibration() function.

plotCalibration(model.perf, quantiles = 5)

Provider profiling analyses

The calcPerformance() function can be used to obtain some basic provider profiling results.

This analysis generates confidence intervals for entity-level performance. Confidence intervals for the observed-to-expected (O/E) standardized rates are Wilson score intervals, while the confidence intervals for the predicted-to-expected (P/E) standardized rates use a parametric bootstrapping method. The number of bootstraps and the number of cores to use for parallel processing can be adjusted with the controlPerf() function within calcPerformance().

perf.out <- calcPerformance(df = example_df_2, model = model, ctrPerf = controlPerf(n.boots = 1000, n.cores = 2))
tutorial_profiling_results <- perf.out$perf.results

The estimated marginal probability of the outcome is 0.22.

Table: Performance summary statistics across entities

Number of observations per entity

plotN(tutorial_profiling_results$n)

Unadjusted rates across entities

plotPerformance(df = tutorial_profiling_results)

OE-risk-standardized rates across entities

plotPerformance(df = tutorial_profiling_results, plot.type = 'oe')

Categorization of entities according to whether their risk-standardized rates are lower, higher, or no different from the overall, unadjusted average can be found in the perf.results output from calcPerformance() function.

Categorization of entities based on OE-risk-standardized rates
Category	Number of entities
Higher than average	16
Lower than average	12
No different than average	72

PE-risk-standardized rates across entities

plotPerformance(df = tutorial_profiling_results, plot.type = 'pe')

Categorization of entities based on PE-risk-standardized rates
Category	Number of entities
Higher than average	1
No different than average	99

Plot of ORs comparing each entity with the average

Another way to assess whether entities have measure performance that differs from the average is to examine the predicted random intercept values. Entities highlighted in red have estimated ORs that are statistically different from 1.0 at alpha = 0.05. There are 19 entities with outcome adjusted odds that are significantly different from the average entity performance.

plotPerformance(df = tutorial_profiling_results, plot.type = 'OR')

Calculate reliability from all methods

Again, to calculate reliability estimates from all methods, we use the calcReliability() function. This time we specify the risk adjustment model as follows:

tutorial_reliability_results_2 <- calcReliability(df = example_df_2, model = model)

Note that while the results from the Beta-Binomial method are included, these estimates currently do not account for risk adjustment variables.

Reliability estimates
Method	Reliability	Min Reliability	Max Reliability
Permutation split-sample	0.682	NA	NA
Hierarchical logit regression	0.691	0.608	0.759
Beta-Binomial	0.686	0.622	0.760

plotReliability(tutorial_reliability_results_2)

Tutorial of QualityMeasure functions

Kenneth Nieser

Example #1: Measures without risk adjustment

Simulate data

Provider profiling analyses

Number of observations per entity

Unadjusted rates across entities

Plot of ORs comparing each entity with the average

Beta-binomial-based reliability estimates

Calculate reliability from all methods

Example #2: Measures with risk adjustment

Simulate data

Fit risk adjustment model to the data

Estimates of model parameters

Discrimination

Calibration

Provider profiling analyses

Number of observations per entity

Unadjusted rates across entities

OE-risk-standardized rates across entities

PE-risk-standardized rates across entities

Plot of ORs comparing each entity with the average

Calculate reliability from all methods