help for qqvalue (Roger Newson)

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, 2010a.
A presentation on qqvalue is given in Newson, 2010b.
Multiple-test procedures are reviewed in Newson et al, 2003a,
and described in the on-line help for multproc.
A presentation on multiple-test procedures
is given in Newson et al, 2003b.
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.
See Newson (2013) for the justification for this practice.

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, Imperial College London, UK.
Email: r.newson@imperial.ac.uk

References

Newson, R. B. 2013.
Bonferroni and Holm approximations for idák and Holland-Copenhaver q-values.
The Stata Journal 13(2): 379-381.
Download from
The Stata Journal website.

Newson, R. B. 2010a.
Frequentist q-values for multiple-test procedures.
The Stata Journal 10(4): 568-584.
Download from
The Stata Journal website.

Newson, R. B. 2010b.
Post-parmest peripherals: fvregen, invcise, and qqvalue.
Presented at
the 16th UK Stata User Meeting, 9-10 September, 2010.

Newson, R. and the ALSPAC Study Team. 2003a.
Multiple-test procedures and smile plots.
The Stata Journal 3(2): 109-132.
Download from
The Stata Journal website.

Newson, R. and the ALSPAC Study Team. 2003b.
Multiple test procedures and smile plots.
Presented at
the 9th UK Stata User Meeting, 19-20 May, 2003.

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, fvregen, invcise if installed