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.

plotN(perf.results$n)

plotPerformance(df = perf.results)

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

plotReliability(tutorial_reliability_results_1)

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)

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.

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

plotPredictedDistribution(model.perf)

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

plotN(tutorial_profiling_results$n)

plotPerformance(df = tutorial_profiling_results)

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

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

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.

plotReliability(tutorial_reliability_results_2)