Generate frequentist q-values by inverting multiple-test procedures
qqvalue varname [if] [in] [ , method(method_name) bestof(#) qvalue(newvarname) npvalue(newvarname) rank(newvarname) svalue( newvarname) rvalue(newvarname) float fast ]
where method_name is one of
bonferroni | sidak | holm | holland | hochberg | simes | yekutieli
by varlist: can be used with qqvalue. (See help for by.) If by varlist: is used, then all generated variables are calculated using the specified multiple-test procedure within each by-group defined by the variables in the varlist.
Description
qqvalue is similar to the R package p.adjust. It inputs a single variable, assumed to contain P-values calculated for multiple comparisons, in a dataset with 1 observation per comparison. It outputs a new variable, containing the frequentist q-values corresponding to these P-values, calculated by inverting a multiple-test procedure specified by the user. These q-values represent, for each corresponding P-value, the minimum uncorrected P-value threshold for which that P-value would be in the discovery set, assuming that the specified multiple-test procedure was used on the same set of input P-values to generate a corrected P-value threshold. These minimum uncorrected P-value thresholds may represent familywise error rates or false discovery rates, depending on the procedure used. Optionally, qqvalue may output other variables, containing the various intermediate results used in calculating the q-values. The multiple-test procedures available for qqvalue are a subset of those available using the multproc module of the smileplot package, which can be downloaded from SSC.
Options for qqvalue
method(method_name) specifies the multiple-test procedure method to be used for calculating the q-values from the input P-values. The method_name may be bonferroni, sidak, holm, holland, hochberg, simes, or yekutieli. These method names specify that the q-values will be calculated from the input P-values by inverting the multiple-test procedure specified by the method() option of the same name for the multproc option of the smileplot package, which can be downloaded from SSC. If method() is unset, then it is set to bonferroni.
bestof(#) specifies an integer number. If the bestof() option is specified (and is greater than the number of input P-values), then the q-values are calculated assuming that the input P-values are a subset (usually the smallest) of a superset of P-values. If the method() option specifies a one-step method (such as bonferroni or sidak), then the q-values do not depend on the other P-values in the superset, but only on the number of P-values in the superset. If the method() option specifies a step-down method (such as holm or holland), then it is assumed that all the other P-values in the superset are greater than the largest of the input P-values. If the method() option specifies a step-up method (such as hochberg, simes, or yekutieli), then it is assumed that all the other P-values in the superset are equal to 1, implying that the q-values will be conservative, and define an upper bound to the respective q-values that would have been calculated, if we knew the other P-values in the superset. If bestof() is unspecified (or non-positive), then the input P-values are assumed to be the full set of P-values calculated. The bestof() option is useful if the input P-values are known (or suspected) to be the smallest of a greater set of P-values, which we do not know. This often happens if the input P-values are from a genome scan reported in the literature.
qvalue(newvarname) specifies the name of a new output variable to be generated, containing the q-values calculated from the input P-values, using the multiple-test procedure specified by the method() option.
npvalue(newvarname) specifies the name of a new output variable to be generated, containing, in each observation, the total number of P-values in the sample of observations specified by the if and in qualifiers, or in the by-group containing that observation, if the by: prefix is specified.
rank(newvarname) is the name of a new variable to be generated, containing, in each observation, the rank of the corresponding P-value, from the lowest to the highest. Tied P-values are ranked according to their position in the input dataset. If the by: prefix is specified, then the ranks are defined within the by-group.
svalue(newvarname) specifies the name of a new output variable to be generated, containing the s-values calculated from the input P-values. The s-values are an intermediate result, calculated in the course of calculating the q-values, and are used mainly for validation. They are calculated from the input P-values by inverting the formulas used for the rank-specific critical P-value thresholds calculated by the multproc module of the smileplot package. These rank-specific P-value thresholds are returned in the generated variable specified by the critical() option of multproc. Note that the s-values may have values greater than 1.
rvalue(newvarname) specifies the name of a new output variable to be generated, containing the r-values calculated from the input P-values. The r-values are an intermediate result, calculated in the course of calculating the q-values, and are used mainly for validation. They are calculated from the s-values by truncating the s-values to a maximum of 1. The q-values are calculated from the r-values using a procedure dependent on the multiple-test procedure specified by the method() option. If the multiple-test procedure is a one-step procedure, such as bonferroni or sidak, then the q-values are equal to the corresponding r-values. If the multiple-test procedure is a step-down procedure, such as holm or holland, then the q-value for each P-value is equal to the cumulative maximum of all the r-values corresponding to P-values of rank equal to or less than that P-value. If the multiple-test procedure is a step-up procedure, such as hochberg, simes or yekutieli, then the q-value for each P-value is equal to the cumulative minimum of all the r-values corresponding to P-values of rank equal to or greater than that P-value.
float specifies that the output variables specified by the qvalue(), rvalue() and svalue() options will be created as variables of type float. If float is absent, then these variables are created as variables of type double. Whether or not float is specified, all generated variables are stored to the lowest precision possible without loss of information.
fast is an option for programmers. It specifies that qqvalue will not take any action so that it can restore the original data in the event of failure, or if the user presses Break.
Remarks
The methods and formulas for qqvalue are given in Newson (2010). Multiple-test procedures are reviewed in Newson et al. (2003), and described in the on-line help for multproc. All of these sources contain extensive references for further reading.
The qqplot package is similar to the R package p.adjust, which also calculates frequentist q-values corresponding to multiple-test procedures. Note that, in the on-line documentation for p.adjust in R, the q-values are referred to as "adjusted P-values", although a lot of users refer to them as "q-values". There is no clear consensus regarding the correct terminology to use, even among statisticians. The term "q-value" was introduced in Storey (2003) to describe a minimum positive false discovery rate (pFDR) under which a P-value will be included in a discovery set, assuming that this discovery set is defined to control the pFDR. The pFDR is a quantity defined for empirical Bayesian methods. By contrast, the multiple-test procedures used by qqvalue and p.adjust define the discovery set to control either the familywise error rate (FWER) or the false discovery rate (FDR), both of which are defined for purely frequentist methods. For this reason, I originally used the term "quasi-q-values" to denote frequentist q-values, and chose the name qqvalue for the package to compute these. However, I was later advised that the prefix "quasi-" was not really necessary. I therefore now simply use the term "q-values", or "frequentist q-values" if I need to distinguish them from Bayesian q-values.
qqvalue, multproc and smileplot all require input datasets with 1 observation for each of a set of P-values, usually corresponding to a set of estimated parameters. Such input datasets may be produced using the official Stata utilities statsby and postfile, or alternatively by the user-written Stata package parmest, which can be downloaded from SSC.
Technical note
If the user specifies method(sidak), then qqvalue uses the formula for method(bonferroni) to calculate the s-values, r-values and q-values corresponding to input P-values too small to be subtracted from 1 in double precision to give a result less than 1. Similarly, if the user specifies method(holland), then qqvalue uses the formula for method(holm) to calculate the s-values, r-values and q-values corresponding to input P-values too small to be subtracted from 1 in double precision to give a result less than 1. This practice ensures that the q-values will be sensible, and no smaller than the corresponding P-values, even if these P-values are too small to be subtracted from 1 to give a result less than 1. It works because the Sidak q-value converges in ratio to the Bonferroni q-value in the limit as the corresponding P-value tends to zero.
Examples
The following example uses the auto data, distributed with Stata. The somersd package is used to measure the Somers' D parameters for rank associations between a list of car-specific variables and non-US origin. The parmest package is then used to replace the dataset in memory with a new dataset, with 1 observation per estimated parameter and data on parameter estimates, confidence limits and P-values. We then use qqvalue to calculate q-values corresponding to the P-values, using the Simes procedure. The parmest and somersd packages can be downloaded from SSC.
. sysuse auto, clear . somersd foreign price mpg headroom trunk weight length turn displacement gear_ratio, tdist . parmest, norestore . qqvalue p, method(simes) qvalue(myqval) . list
The following example also uses the auto data. It first uses the somersd package, together with the parmby module of the parmest package, to create a new dataset in the memory, with 1 observation for each of a list of rank correlations involving car price in each car origin group (US and non-US cars). We then use qqvalue, with the by varlist: prefix, to demonstrate the calculation of 2 separate sets of q-values, one for US-made cars and one for non-US-made cars.
. sysuse auto, clear . parmby "somersd price mpg headroom trunk weight length turn displacement gear_ratio, tdist", by(foreign) norestore . by foreign: qqvalue p, method(simes) qvalue(myqval) . by foreign: list
Author
Roger B. Newson, National Heart and Lung Institute, Imperial College London, UK. Email: r.newson@imperial.ac.uk
References
Newson, R. B. 2010. Frequentist q-values for multiple-test procedures. The Stata Journal 10(4): 568-584. Download from THe Stata Journal website.
Newson, R. and the ALSPAC Study Team. 2003. Multiple-test procedures and smile plots. The Stata Journal 3(2): 109-132. Download from THe Stata Journal website.
Storey, J. D. 2003. The positive false discovery rate: a Bayesian interpretation and the q-value. The Annals of Statistics 31(6): 2013–2035.
Also see
Manual: [R] by, [R] statsby, [P] postfile. On-line: help for [D] by, [D] statsby, [P] postfile help for multproc, smileplot, parmest, somersd if installed