Package {posc}

posc: Probability of Outcome Superiority Curves

Description

A probability of outcome superiority curve (POSC) shows, for two people who differ by \Delta on a predictor, the probability that the person with the higher predictor score also has the higher outcome:

P(\Delta) = P(Y_1 > Y_2 \mid X_1 - X_2 = \Delta).

Details

Because the two people are exchangeable, every POSC satisfies P(-\Delta) = 1 - P(\Delta): it passes through 0.5 at \Delta = 0 and is point-symmetric about that point. The curve for \Delta \ge 0 therefore carries all of the information, and is what is plotted by default (use plot(x, full = TRUE) to see the full sigmoid).

Functions

• posc_uni(): POSC from a correlation, assuming bivariate normality.

• posc_multi(): POSC for an optimally weighted composite of several predictors, from a correlation matrix.

• posc_lm(): POSC for the fitted values of a linear model.

• posc_emp(): empirical POSC estimated from raw data without assuming normality, with bootstrap confidence bands.

• posc_comp(): overlay and compare several POSCs.

• posc_bins(): binned observed proportions for an empirical POSC.

posc_uni(), posc_multi(), posc_lm() and posc_emp() return a "posc" object with print(), summary(), plot(), ggplot2::autoplot(), predict() and as.data.frame() methods.

Author(s)

Maintainer: Matthew B. Jané matthewbjane@gmail.com

Authors:

• Matthew B. Jané matthewbjane@gmail.com

See Also

Useful links:

• https://matthewbjane.com/posc/

• https://github.com/MatthewBJane/posc

• Report bugs at https://github.com/MatthewBJane/posc/issues

Coerce a POSC to a data frame

Description

Coerce a POSC to a data frame

Usage

## S3 method for class 'posc'
as.data.frame(x, row.names = NULL, optional = FALSE, full = TRUE, ...)

Arguments

Value

A data frame with columns delta, p, lower and upper.

Examples

head(as.data.frame(posc_uni(r = 0.5, n = 100), full = FALSE))

Plot a probability of outcome superiority curve

Description

autoplot() returns a ggplot2::ggplot() object that can be modified further; plot() draws it.

Usage

## S3 method for class 'posc'
autoplot(
  object,
  full = FALSE,
  color = "#1F5A96",
  ci = TRUE,
  reference = FALSE,
  bins = FALSE,
  percent = FALSE,
  grid = TRUE,
  mark = NULL,
  mark_ci = FALSE,
  title = NULL,
  subtitle = NULL,
  ...
)

## S3 method for class 'posc'
plot(x, ...)

Arguments

Value

autoplot() returns a ggplot object; plot() returns it invisibly.

Examples

fit <- posc_uni(r = 0.5, n = 120, sd_x = 15,
                predictor_name = "IQ", outcome_name = "job performance")
plot(fit)
plot(fit, full = TRUE)
plot(fit, mark = c(10, 20, 30))
plot(fit, mark = 15, mark_ci = TRUE, percent = TRUE)
plot(fit, title = TRUE, subtitle = TRUE)

Plot a comparison of POSCs

Description

Plot a comparison of POSCs

Usage

## S3 method for class 'posc_comparison'
autoplot(
  object,
  full = FALSE,
  colors = NULL,
  ci = TRUE,
  percent = FALSE,
  grid = TRUE,
  mark = NULL,
  mark_ci = FALSE,
  title = NULL,
  subtitle = NULL,
  ...
)

## S3 method for class 'posc_comparison'
plot(x, ...)

Arguments

Value

autoplot() returns a ggplot object; plot() returns it invisibly.

Examples

cmp <- posc_comp(a = posc_uni(0.5, n = 100), b = posc_uni(0.3, n = 100))
plot(cmp)
plot(cmp, mark = 1, title = TRUE, subtitle = TRUE)

Binned observed proportions for an empirical POSC

Description

Returns the observed proportion of concordant pairs (the person higher on the predictor is also higher on the outcome) in bins of the absolute difference in the predictor, with simultaneous (family-wise) bootstrap confidence intervals. These are the points and whiskers drawn by plot(fit, bins = TRUE) for curves from posc_emp(); see its Details for the method.

Usage

posc_bins(object, level = object$level)

Arguments

Value

A data frame with one row per bin: delta_low and delta_high (bin limits), delta (mean absolute difference of the pairs in the bin), n_pairs, p (observed proportion), and lower and upper (simultaneous confidence limits; NA if the curve was fitted with n_boot = 0). Bins that contain no pairs have NA estimates.

Examples

set.seed(1)
x <- rnorm(100)
y <- x + rnorm(100)
fit <- posc_emp(x, y, n_boot = 200, bins = 6)
posc_bins(fit)

Compare probability of outcome superiority curves

Description

Overlays several POSCs and tests whether their underlying correlations differ.

Usage

posc_comp(..., labels = NULL)

Arguments

Details

Under bivariate normality a POSC is determined entirely by its correlation, so curves from independent samples can be compared by testing whether the correlations are equal. For two curves this is the usual Fisher-z test,

z = \frac{\tanh^{-1} r_1 - \tanh^{-1} r_2} {\sqrt{1/(n_1 - 3) + 1/(n_2 - 3)}};

for more than two it is the homogeneity test Q = \sum_i (n_i - 3)(z_i - \bar z)^2, referred to a \chi^2 distribution with one fewer degrees of freedom than curves, where \bar z is the weighted mean of the z_i.

Both tests assume the curves come from independent samples. For curves from posc_emp() the test compares the Pearson correlations, not the shapes of the curves. No test is reported if any curve lacks a sample size.

Value

An object of class "posc_comparison", a list with elements posc (the curves), table (a data frame of each curve's method, correlation and sample size) and test (a list with the statistic, degrees of freedom and p-value, or NULL). It has print(), plot() and ggplot2::autoplot() methods.

Examples

a <- posc_uni(r = 0.5, n = 100)
b <- posc_uni(r = 0.3, n = 120)
cmp <- posc_comp(`Structured interview` = a, `Unstructured interview` = b)
cmp
plot(cmp)

Empirical POSC from raw data

Description

Estimates the probability of outcome superiority curve directly from predictor and outcome scores, without assuming bivariate normality, with a bootstrap confidence band.

Usage

posc_emp(
  x,
  y,
  df = 3,
  level = 0.95,
  n_boot = 1000,
  max_pairs = 20000,
  max_delta = NULL,
  bins = 10,
  standardize = FALSE,
  predictor_name = "X",
  outcome_name = "Y"
)

Arguments

Details

Pairwise data. Every pair of people (i, j) with x_i \ne x_j contributes the difference \Delta_{ij} = x_i - x_j and an indicator of whether the person with the higher predictor score also has the higher outcome (ties in the outcome count as 1/2). For large samples a random subset of max_pairs pairs is used.

Model. The curve is modelled as

\mathrm{logit}\, P(\Delta) = \mathrm{sign}(\Delta)\, g(|\Delta|),

where g is a natural cubic spline with df degrees of freedom and no intercept, whose basis is zero at |\Delta| = 0. This builds in the two properties every POSC must have, P(0) = 0.5 and P(-\Delta) = 1 - P(\Delta), while leaving the shape free. The 0/1 responses are shrunk by 0.0005 towards 0.5 so that the fit exists even for near-perfect relationships; this changes the curve by a negligible amount. Increase df for a more flexible curve; df = 1 makes the logit of the curve linear in \Delta. Because a natural spline has zero curvature at its boundary, the curve is smooth through \Delta = 0. It is not forced to be monotone.

Confidence band. Pairs that share a person are not independent, so model-based standard errors from the pairwise regression would be far too small. Instead, people are resampled with replacement n_boot times, the curve is refitted to each resample (with the knots held fixed), and the band is formed from pointwise percentiles of the bootstrap curves. Pairs formed by a person and their own duplicate are excluded. Use set.seed() beforehand for reproducible results.

Binned observed proportions. To check the fitted curve against the data, the pairs with |\Delta| \le max_delta are grouped into bins of |\Delta| and the observed proportion of pairs in which the person higher on the predictor is also higher on the outcome is computed in each bin. plot(fit, bins = TRUE) shows these as points with whiskers, and posc_bins() returns them as a table. The whiskers are simultaneous (family-wise) confidence intervals: with probability level they cover all of the bin proportions at once. They are computed from the same person-level bootstrap as the curve, on the empirical-logit scale \ell_k = \mathrm{logit}\{(n_k \hat p_k + 0.5)/(n_k + 1)\} for a bin with n_k pairs: the interval is \mathrm{logit}^{-1}(\ell_k \pm c\,\widehat{se}_k), where \widehat{se}_k is the bootstrap standard error of \ell_k and c is the level quantile of \max_k |\ell^*_k - \ell_k| / \widehat{se}_k over the bootstrap resamples (the max-t or sup-t method). The intervals are therefore asymmetric and lie between 0 and 1.

Non-monotone relationships. A POSC averages over where in the predictor's distribution a pair of people sits. For a U-shaped relationship, the person higher on the predictor is more likely to have the higher outcome on one side of the minimum and less likely on the other, so the curve reflects the balance of the two. For a symmetric U centred on the middle of the predictor's distribution the true POSC is flat at 0.5, even though the predictor is strongly related to the outcome. The empirical curve is not constrained to be monotone and recovers such shapes; the dashed normal-theory reference, which depends only on the linear correlation, is not meaningful in that case.

Value

An object of class "posc": a list with elements

curve

data frame with columns delta, p, lower, upper, over a symmetric grid of differences.

method, r, n, level, sd_x, max_delta, labels, details

information about how the curve was computed.

Use plot(), summary() and predict() to work with it.

See Also

posc_uni() for the curve implied by bivariate normality, which plot(fit, reference = TRUE) overlays as a dashed line.

Examples

set.seed(1)
x <- rnorm(120)
y <- exp(0.8 * x + rnorm(120, sd = 0.6))   # skewed, non-normal outcome

fit <- posc_emp(x, y, n_boot = 200)   # use n_boot >= 1000 in practice
fit
plot(fit)
posc_bins(fit)
plot(fit, bins = TRUE, reference = TRUE)   # add data and normal curve

POSC for the fitted values of a linear model

Description

Computes the probability of outcome superiority curve for a linear model's predictions: for two people whose predicted outcomes differ by \Delta, the probability that the one predicted higher actually scores higher.

Usage

posc_lm(
  model,
  adjust = TRUE,
  level = 0.95,
  max_delta = NULL,
  predictor_name = NULL,
  outcome_name = NULL
)

Arguments

Details

The fitted values are treated as a single predictor whose correlation with the observed outcome is the multiple correlation R of the model, and the curve is computed as in posc_uni() under bivariate normality. Differences are expressed in the units of the outcome (that is, of the fitted values). With adjust = TRUE (the default) R is replaced by the square root of the adjusted R^2; see posc_multi() for why.

Value

An object of class "posc": a list with elements

curve

data frame with columns delta, p, lower, upper, over a symmetric grid of differences.

method, r, n, level, sd_x, max_delta, labels, details

information about how the curve was computed.

Use plot(), summary() and predict() to work with it.

Examples

set.seed(1)
d <- data.frame(x1 = rnorm(150), x2 = rnorm(150))
d$y <- d$x1 + 0.5 * d$x2 + rnorm(150)
fit <- posc_lm(lm(y ~ x1 + x2, data = d))
fit
plot(fit)

POSC for a composite of several predictors (bivariate normal)

Description

Computes the probability of outcome superiority curve for the optimally weighted (regression) composite of several predictors, from a correlation matrix.

Usage

posc_multi(
  r_mat,
  n,
  outcome = ncol(r_mat),
  adjust = TRUE,
  level = 0.95,
  max_delta = 3,
  predictor_name = "Predicted Y",
  outcome_name = NULL,
  outcome_idx = NULL
)

Arguments

Details

The composite that best predicts the outcome correlates with it at the multiple correlation

R = \sqrt{\mathbf{r}_{xy}^\top \mathbf{R}_{xx}^{-1} \mathbf{r}_{xy}},

where \mathbf{R}_{xx} is the predictor intercorrelation matrix and \mathbf{r}_{xy} the predictor-outcome correlations. The POSC is then computed as in posc_uni(), with differences in the composite expressed in its standard-deviation units.

A sample multiple correlation is biased upwards as an estimate of the population multiple correlation, because the weights are fitted to the same sample. With adjust = TRUE (the default) R is replaced by the square root of the adjusted R^2, 1 - (1 - R^2)(n - 1)/(n - k - 1) for k predictors, which corrects most of this bias. (How well weights estimated in this sample would predict in a new sample is lower still.)

The confidence band uses a Fisher-z interval for R with standard error 1/\sqrt{n - k - 2}, truncated at zero, and should be treated as approximate.

Value

An object of class "posc": a list with elements

curve

data frame with columns delta, p, lower, upper, over a symmetric grid of differences.

method, r, n, level, sd_x, max_delta, labels, details

information about how the curve was computed.

Use plot(), summary() and predict() to work with it.

See Also

posc_lm() to start from a fitted linear model instead.

Examples

r_mat <- matrix(c(1.0, 0.4, 0.3, 0.2,
                  0.4, 1.0, -0.3, 0.2,
                  0.3, -0.3, 1.0, 0.3,
                  0.2, 0.2, 0.3, 1.0), 4, 4,
                dimnames = list(c("X1", "X2", "X3", "Y"),
                                c("X1", "X2", "X3", "Y")))
fit <- posc_multi(r_mat, n = 100, outcome = "Y")
fit
plot(fit)

POSC from a single correlation (bivariate normal)

Description

Computes the probability of outcome superiority curve implied by a correlation between a predictor and an outcome, assuming the two are bivariate normal.

Usage

posc_uni(
  r,
  n = NULL,
  sd_x = 1,
  level = 0.95,
  max_delta = 3 * sd_x,
  predictor_name = "X",
  outcome_name = "Y"
)

Arguments

Details

If X and Y are bivariate normal with correlation r, then for two randomly chosen people whose predictor scores differ by \Delta

P(Y_1 > Y_2 \mid X_1 - X_2 = \Delta) = \Phi\left(\frac{r\,\Delta / \sigma_X}{\sqrt{2(1 - r^2)}}\right),

where \Phi is the standard normal distribution function and \sigma_X the standard deviation of the predictor. Set sd_x to put \Delta on the raw scale of the predictor; the default sd_x = 1 gives differences in standard-deviation units.

The curve is monotone in r, so the pointwise confidence band is obtained by evaluating it at the limits of the Fisher-z confidence interval for r.

Value

An object of class "posc": a list with elements

curve

data frame with columns delta, p, lower, upper, over a symmetric grid of differences.

method, r, n, level, sd_x, max_delta, labels, details

information about how the curve was computed.

Use plot(), summary() and predict() to work with it.

See Also

posc_emp() for a curve estimated from raw data without assuming normality; posc_comp() to compare curves.

Examples

# Differences in SD units
fit <- posc_uni(r = 0.5, n = 100)
fit
plot(fit)

# Raw units: an IQ test with SD = 15
iq <- posc_uni(r = 0.5, n = 100, sd_x = 15, predictor_name = "IQ",
               outcome_name = "job performance")
predict(iq, delta = c(5, 15, 30))

Predicted probabilities from a POSC

Description

Evaluates a probability of outcome superiority curve, with its confidence interval, at any differences in the predictor.

Usage

## S3 method for class 'posc'
predict(object, delta, level = object$level, ...)

Arguments

Value

A data frame with columns delta, p (probability that person 1 has the higher outcome), lower and upper (confidence limits; NA when no interval is available).

Examples

fit <- posc_uni(r = 0.5, n = 100, sd_x = 15)
predict(fit, delta = c(5, 10, 15))

Print and summarise a POSC

Description

print() gives a compact description of the curve and its probabilities at a few differences in the predictor. summary() returns those probabilities as a data frame, at differences of your choosing.

Usage

## S3 method for class 'posc'
print(x, digits = 3, ...)

## S3 method for class 'posc'
summary(object, delta = NULL, ...)

Arguments

Value

print() returns x invisibly. summary() returns a data frame as described in predict.posc().

Examples

fit <- posc_uni(r = 0.4, n = 250)
fit
summary(fit, delta = c(0.25, 0.5, 1))

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

ggplot2

autoplot()