R: Nonparametric Bayesian inference for the covariate-specific...

cROC.bnp {ROCnReg}

R Documentation

Nonparametric Bayesian inference for the covariate-specific ROC curve (cROC).

Description

This function estimates the covariate-specific ROC curve (cROC) using the nonparametric Bayesian approach proposed by Inacio de Carvalho et al. (2013).

Usage

cROC.bnp(formula.h, formula.d, group, tag.h, data, 
    newdata, standardise = TRUE, p = seq(0, 1, l = 101), ci.level = 0.95,
    compute.lpml = FALSE, compute.WAIC = FALSE, compute.DIC = FALSE, 
    pauc = pauccontrol(), density = densitycontrol(),
    prior.h = priorcontrol.bnp(), prior.d = priorcontrol.bnp(), 
    mcmc = mcmccontrol(),
    parallel = c("no", "multicore", "snow"), ncpus = 1, cl = NULL)

Arguments

`formula.h`	A `formula` object specifying the regression function associated to each component of the single-weights linear dependent Dirichlet process mixture of normals model used to estimate the conditional distribution function of the diagnostic test outcome in the healthy population. Regarding the modelling of continuous covariates, both linear and nonlinear effects are allowed, with nonlinear effects being modelled through B-spline basis expansions (see Note).
`formula.d`	A `formula` object specifying the regression function associated to each component of the single weights linear dependent Dirichlet process mixture model used to estimate the conditional distribution function of the diagnostic test outcome in the diseased population. Both linear and nonlinear (through the use of smooth functions approximated by linear combinations of cubic B-splines basis functions) covariate effects are allowed (see Note).
`group`	A character string with the name of the variable that distinguishes healthy from diseased individuals.
`tag.h`	The value codifying healthy individuals in the variable `group`.
`data`	A data frame representing the data and containing all needed variables.
`newdata`	Optional data frame containing the values of the covariates at which the covariate-specific ROC curve (AUC and pAUC, if computed) will be computed. If not supplied, the function `cROCData` is used to build a default dataset.
`standardise`	A logical value. If TRUE both the test outcomes and the continuous covariates assumed to have a linear effect are standardised (i.e., the resulting variables have mean zero and standard deviation of one). The default is TRUE.
`p`	Set of false positive fractions (FPF) at which to estimate the covariate-specific ROC curve.
`ci.level`	An integer value (between 0 and 1) specifying the level for the credible interval. The default is 0.95.
`compute.lpml`	A logical value. If TRUE, the log pseudo marginal likelihood (LPML, Geisser and Eddy, 1979) and the conditional predictive ordinates (CPO) are computed.
`compute.WAIC`	A logical value. If TRUE, the widely applicable information criterion (WAIC, Gelman et al., 2014; Watanabe, 2010) is computed.
`compute.DIC`	A logical value. If TRUE, the deviance information criterion (DIC)(Celeux et al., 2006, Spiegelhalter et al., 2002) is computed.
`pauc`	A list of control values to replace the default values returned by the function `pauccontrol`. This argument is used to indicate whether the partial area under the covariate-adjusted ROC curve should be computed, and in case it is computed, whether the focus should be placed on restricted false positive fractions (FPFs) or on restricted true positive fractions (TPFs), and the upper bound for the FPF (if focus is FPF) or the lower bound for the TPF (if focus is TPF).
`density`	A list of control values to replace the default values returned by the function `densitycontrol`. This argument is used to indicate whether the conditional densities of the marker in the healthy and diseased population should be computed, and in case it is to be computed, at which grid of test outcomes in each of the populations the conditional densities should be evaluated.
`prior.h`	Hyparameter specification for the healthy population. A list of control values to replace the default values returned by the function `priorcontrol.bnp`. See `priorcontrol.bnp` for details.
`prior.d`	Hyparameter specification for the diseased population. A list of control values to replace the default values returned by the function `priorcontrol.bnp`. See `priorcontrol.bnp` for details.
`mcmc`	A list of control values to replace the default values returned by the function `mcmccontrol`. See `mcmccontrol` for details.
`parallel`	A characters string with the type of parallel operation: either "no" (default), "multicore" (not available on Windows) or "snow".
`ncpus`	An integer with the number of processes to be used in parallel operation. Defaults to 1.
`cl`	An object inheriting from class `cluster` (from the `parallel` package), specifying an optional parallel or snow cluster if parallel = "snow". If not supplied, a cluster on the local machine is created for the duration of the call.

Details

Estimates the covariate-specific ROC curve (cROC) defined as

ROC(p|\mathbf{x}) = 1 - F_{D}\{F_{\bar{D}}^{-1}(1-p|\mathbf{x})|\mathbf{x}\},

where

F_{D}(y|\mathbf{x}) = Pr(Y_{D} \leq y | \mathbf{X}_{D} = \mathbf{x}),

F_{\bar{D}}(y|\mathbf{x}) = Pr(Y_{\bar{D}} \leq y | \mathbf{X}_{\bar{D}} = \mathbf{x}).

Note that, for the sake of clarity, we assume that the covariates of interest are the same in both the healthy and diseased populations. The method implemented in this function estimates F_{D}(\cdot|\mathbf{x}) and F_{\bar{D}}(\cdot|\mathbf{x}) by means of a single-weights linear dependent Dirichlet process mixture of normals model (De Iorio et al., 2009). More precisely, and letting \{(\mathbf{x}_{\bar{D}i},y_{\bar{D}i})\}_{i=1}^{n_{\bar{D}}} and \{(\mathbf{x}_{Dj},y_{Dj})\}_{j=1}^{n_{D}} be two independent random samples from the nondiseased and diseased populations, respectively, our postulated model for the conditional distribution in each group function takes the following form

F_{\bar{D}}(y_{\bar{D}i}|\mathbf{X}_{\bar{D}}=\mathbf{x}_{\bar{D}i}) = \sum_{l=1}^{L_{\bar{D}}}\omega_{l\bar{D}}\Phi(y_{\bar{D}i}\mid\mu_{l\bar{D}}(\mathbf{x}_{\bar{D}i}),\sigma_{l\bar{D}}^2),

F_{D}(y_{Dj}|\mathbf{X}_{D} = \mathbf{x}_{Dj}) = \sum_{l=1}^{L_{D}}\omega_{lD}\Phi(y_{Dj}\mid\mu_{lD}(\mathbf{x}_{\bar{D}i}),\sigma_{lD}^2),

where \Phi(y|\mu, \sigma^2) denotes the cumulative distribution function of the normal distribution, evaluated at y, with mean mu and variance \sigma^2. The regression function \mu_{ld}(\mathbf{x}_{di}) can incorportate both linear and nonlinear (through B-splines) effects of continuous covariates, categorical covariates (factors) as well as interactions. Interactions between categorical and (nonlinear) continuous covariates are also allowed (factor-by curve interactions). For the sake of simplicity we write \mu_{ld}(\mathbf{x}_{di}) = \mathbf{z}_{di}^{T}\mathbf{\beta}_{ld}, where \mathbf{z}_{di} is the ith column of the design matrix (possibly containing a basis representation of some/all continuous covariates), d \in \{D, \bar{D}\}. Here L_d is a pre-specified upper bound on the number of mixture components. The \omega_{ld}'s result from a truncated version of the stick-breaking construction (\omega_{1d} = v_{1d}; \omega_{ld} = v_{ld}\prod_{r<l}(1-v_{dr}), l=2,\ldots,L_{d}; v_{d1},\ldots,v_{L_{d}-1}\sim Beta (1,\alpha_{d}); v_{Ld} = 1, \alpha_d \sim \Gamma(a_{\alpha_d},b_{\alpha_d})), \mathbf{\beta}_{ld}\sim N_{Q_d}(\mathbf{m}_{d},\mathbf{S}_{d}), and \sigma_{ld}^{-2}\sim\Gamma(a_{d},b_{d}). It is further assumed that \mathbf{m}_{d} \sim N_{Q_k}(\mathbf{m}_{0d},\mathbf{S}_{0d}) and \mathbf{S}_{d}^{-1}\sim W(\nu,(\nu_k\Psi_d)^{-1}). Here W(\nu,(\nu\Psi)^{-1}) denotes a Wishart distribution with \nu degrees of freedom and expectation \Psi^{-1}, Here \Gamma(a,b) denotes a Gamma distribution with shape parameter a and rate parameter b, and Q_d denotes the dimension of the vector \mathbf{z}_{di}. It is worth mentioning that when L_d=1, the model for the conditional distribution of the test outcomes reduces to a normal regression model (where continuous covariates effects are modelled either parametrically or nonparametrically). For a detailed description, we refer to Inacio de Carvalho et al. (2013).

The covariate-specific area under the curve is

AUC(\mathbf{x})=\int_{0}^{1}ROC(p|\mathbf{x})dp.

When the upper bound on the number of mixture components is 1, i.e., L_d = 1 (d \in \{D, \bar{D}\}), there is a closed-form expression for the covariate-specific AUC (binormal model), which is used in the package. In contrast, when L_D > 1 or L_{\bar{D}} > 1, the integral is computed numerically using Simpson's rule. With regard to the partial area under the curve, when focus = "FPF" and assuming an upper bound u_1 for the FPF, what it is computed is

pAUC_{FPF}(u_1|\mathbf{x})=\int_0^{u_1} ROC(p|\mathbf{x})dp.

As for the AUC, when L_d = 1 (d \in \{D, \bar{D}\}), there is a closed-form expression for the pAUC_{FPF} (Hillis and Metz, 2012), and when L_D > 1 or L_{\bar{D}} > 1 the integral is approximated numerically using Simpson's rule. The returned value is the normalised pAUC, pAUC_{FPF}(u_1|\mathbf{x})/u_1 so that it ranges from u_1/2 (useless test) to 1 (perfect marker). Conversely, when focus = "TPF", and assuming a lower bound for the TPF of u_2, the partial area corresponding to TPFs lying in the interval (u_2,1) is computed as

pAUC_{TPF}(u_2|\mathbf{x})=\int_{u_2}^{1}ROC_{TNF}(p|\mathbf{x})dp,

where ROC_{TNF}(p|\mathbf{x}) is a 270^\circ rotation of the ROC curve, and it can be expressed as ROC_{TNF}(p|\mathbf{x}) = F_{\bar{D}}\{F_{D}^{-1}(1-p|\mathbf{x})|\mathbf{x}\}. Again, when L_d = 1 (d \in \{D, \bar{D}\}), there is a closed-form expression for the pAUC_{TNF} (Hillis and Metz, 2012), and when L_D > 1 or L_{\bar{D}} > 1 the integral is approximated numerically using Simpson's rule. The returned value is the normalised pAUC, pAUC_{TPF}(u_2|\mathbf{x})/(1-u_2), so that it ranges from (1-u_2)/2 (useless test) to 1 (perfect test).

It is worth referring that with respect to the computation of the DIC, when L=1, it is computed as in Spiegelhalter et al. (2002), and when L>1, DIC3 as described in Celeux et al. (2006) is computed. Also, for the computation of the conditional predictive ordinates (CPO) we follow the stable version proposed by Gelman et al. (2014).

Value

As a result, the function provides a list with the following components:

`call`	The matched call.
`newdata`	A data frame containing the values of the covariates at which the covariate-specific ROC curve (as well as the AUC, pAUC, dens and reg.fun, if required) was computed.
`data`	The original supplied data argument.
`missing.ind`	A logical value indicating whether for each pair of observations (test outcomes and covariates) missing values occur.
`marker`	The name of the diagnostic test variable in the dataframe.
`group`	The value of the argument `group` used in the call.
`tag.h`	The value of the argument `tag.h` used in the call.
`p`	Set of false positive fractions (FPF) at which the covariate-specific ROC curve has been estimated.
`ci.level`	The value of the argument `ci.level` used in the call.
`prior`	A list returning the hyperparameter values in the healthy and diseased populations.
`ROC`	A list returning the `np` (length of the vector `p`) by npred predicted covariate-specific ROC curves (cROC) (posterior mean) and `ci.level`*100% pointwise posterior credible bands.
`AUC`	Estimated area under the covariate-specific ROC curve (posterior mean) and `ci.level`*100% posterior credible band.
`pAUC`	If computed, estimated partial area under the covariate-adjusted ROC curve (posterior mean) and `ci.level`*100% credible band. Note that the returned values are normalised, so that the maximum value is one (see more on Details).
`dens`	Named list of length two, with components 'h' (healthy) and 'd' (diseased). Each component is a list with two components: `grid` (grid of test outcomes where the conditional densities were evaluated) and `dens` (MCMC realisations of the corresponding conditional densities).
`reg.fun`	Named list of length two, with components 'h' (healthy) and 'd' (diseased). Each component is a data frame containing the predicted regression function (posterior mean) and `ci.level`*100% credible band.
`lpml`	If computed, named list of length two, with components 'h' (healthy) and 'd' (diseased). Each component is a list with two components: the log pseudo marginal likelihood (LPML) and the conditional predictive ordinates (CPO).
`WAIC`	If computed, named list of length two, with components 'h' (healthy) and 'd' (diseased). Each component is a list with two components: widely applicable information criterion (WAIC) and associated complexity penalty (pW).
`DIC`	If computed, named list of length two, with components 'h' (healthy) and 'd' (diseased). Each component is a list with two components: deviance information criterion (DIC) and associated complexity penalty (pD).
`fit`	Named list of length two, with components `h` (healthy) and `d` (diseased). Each component is a list with the following information: (1) `formula`: the value of the argument `formula.h` or `formula.d` used in the call. (2) `mm`: information needed to construct the model matrix associated with single-weights linear dependent Dirichlet process mixture of normals model. (3) `beta`: array of dimension `nsave`x`L`x`Q` with the sampled regression coefficients. (4) `sd`: matrix of dimension `nsave`x`L` with the sampled variances. (4) `probs`: matrix of dimension `nsave`x`L` with the sampled components' weights. Here, `nsave` is the number of Gibbs sampler iterations saved, `L` is the maximum number of mixture components, and `Q` is the dimension of vector `\mathbf{z}_{d}`, `d \in \{D, \bar{D}\}`. (see also Details). (see also Details).
`data_model`	A list with the data used in the fit: observed diagnostic test outcome and design matrices, separately for the healthy and diseased groups.

Note

The input arguments formula.h and formula.d are similar to that used for the glm function, except that flexible specifications can be added by means of the function f(). For instance, specification y \sim x1 + f(x2, K = 3) would assume a linear effect of x1 (if x1 continuous) and the effect of x2 would be modeled using B-splines basis functions. The argument K = 3 indicates that 3 internal knots will be used, with the quantiles of x2 used for their location. Categorical variables (factors) can be also incorporated, as well as interaction terms. For example, to include the factor-by-curve interaction between age and gender we need to specify, e.g., y \sim gender + f(age, by = gender, K = c(3, 5)). Note that, in this case, the number of knots can be different for each level of the factor. The order of the vector K of knots should match the levels of the factor.

References

Celeux, G., Forbes, F., Robert C. P., and Titerrington, D. M. (2006). Deviance information criteria for missing data models. Bayesian Analysis, 1, 651–674.

De Iorio, M., Johnson, W. O., Muller, P., and Rosner, G. L. (2009). Bayesian nonparametric nonproportional hazards survival modeling. Biometrics, 65, 762–775.

Geisser, S. and Eddy, W.F. (1979) A Predictive Approach to Model Selection, Journal of the American Statistical Association, 74, 153–160.

Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., and Rubin, D.B. (2014). Bayesian Data Analysis, 3rd ed. CRC Press: Boca Raton, FL.

Gelman, A., Hwang, J., and Vehtari, A. (2014). Understanding predictive information criteria for Bayesian models. Statistics and Computing, 24, 997–1010.

Hillis, S. L. and Metz, C.E. (2012). An Analytic Expression for the Binormal Partial Area under the ROC Curve. Academic Radiology, 19, 1491–1498.

Inacio de Carvalho, V., Jara, A., Hanson, T. E., and de Carvalho, M. (2013). Bayesian nonparametric ROC regression modeling. Bayesian Analysis, 8, 623–646.

Speigelhalter, D. J., Best, N. G., Carlin, B. P., and van der Linde, A. (2002). Bayesian measures of model comparison and fit. Journal of the Royal Statistical Society, Ser. B, 64, 583–639.

Watanabe, S. (2010). Asymptotic Equivalence of Bayes Cross Validation and Widely Applicable Information Criterion in Singular Learning Theory. Journal of Machine Learning Research, 11, 3571–3594.

Examples

library(ROCnReg)
data(psa)
# Select the last measurement
newpsa <- psa[!duplicated(psa$id, fromLast = TRUE),]

# Log-transform the biomarker
newpsa$l_marker1 <- log(newpsa$marker1)

cROC_bnp <- cROC.bnp(formula.h = l_marker1 ~ f(age, K = 0),
               formula.d = l_marker1 ~ f(age, K = 0),
               group = "status", 
               tag.h = 0,
               data = newpsa,
               standardise = TRUE, 
               p = seq(0, 1, len = 101),
               compute.lpml = TRUE, 
               compute.WAIC = TRUE,
               compute.DIC = TRUE, 
               pauc = pauccontrol(compute = TRUE, value = 0.5, focus = "FPF"),
               density = densitycontrol(compute = TRUE, grid.h = NA, grid.d = NA),
               mcmc = mcmccontrol(nsave = 500, nburn = 100, nskip = 1))
summary(cROC_bnp)

plot(cROC_bnp)

[Package ROCnReg version 1.0-9 Index]