# Tutorial of QualityMeasure functions

#### Kenneth Nieser

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.

    library(QualityMeasure)

# 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)