Example 61.6 Adaptive Adjustments and ODS Graphics

An experiment was performed using Affymetrix^®gene chips on the CD4 lymphocyte white blood cells of patients with and without a hereditary allergy (atopy) and possibly with asthma. The Asthma-Atopy microarray data set and analysis are discussed in Gibson and Wolfinger (2004): a one-way ANOVA model of the log2mas5 variable ( $\log _2(\mbox{MAS 5.0}$ summary statistics) is fit against a classification variable trt describing different asthma-atopy combinations in the patients, and the least squares means of the trt variable are computed.

For this example, a 1% random sample of least squares means having p-values exceeding 1E–6 is taken. The resulting data are recorded in the test data set, where the Probe_Set_ID variable identifies the probe and the Probt variable contains the p-values for the m = 121 tests, as follows:

data test; 
   length Probe_Set_ID $9.; 
   input Probe_Set_ID $ Probt @@; 
   datalines;
200973_s_ .963316  201059_at .462754  201563_at .000409  201733_at .000819  
201951_at .000252  202944_at .106550  203107_x_ .040396  203372_s_ .010911  
203469_s_ .987234  203641_s_ .019296  203795_s_ .002276  204055_s_ .002328  
205020_s_ .008628  205199_at .608129  205373_at .005209  205384_at .742381  
205428_s_ .870533  205653_at .621671  205686_s_ .396440  205760_s_ .000002  
206032_at .024661  206159_at .997627  206223_at .003702  206398_s_ .191682  
206623_at .010030  206852_at .000004  207072_at .000214  207371_at .000013  
207789_s_ .023623  207861_at .000002  207897_at .000007  208022_s_ .251999  
208086_s_ .000361  208406_s_ .040182  208464_at .161468  209055_s_ .529824  
209125_at .142276  209369_at .240079  209748_at .071750  209894_at .000042  
209906_at .223282  210130_s_ .192187  210199_at .101623  210477_x_ .300038  
210491_at .000078  210531_at .000784  210734_x_ .202931  210755_at .009644  
210782_x_ .000011  211320_s_ .022896  211329_x_ .486869  211362_s_ .881798  
211369_at .000030  211399_at .000008  211572_s_ .269788  211647_x_ .001301  
213072_at .005019  213143_at .008711  213238_at .004824  213391_at .316133  
213468_at .000172  213636_at .097133  213823_at .001678  213854_at .001921  
213976_at .000299  214006_s_ .014616  214063_s_ .000361  214407_x_ .609880  
214445_at .000009  214570_x_ .000002  214648_at .001255  214684_at .288156  
214991_s_ .006695  215012_at .000499  215117_at .000136  215201_at .045235  
215304_at .000816  215342_s_ .973786  215392_at .112937  215557_at .000007  
215608_at .006204  215935_at .000027  215980_s_ .037382  216010_x_ .000354  
216051_x_ .000003  216086_at .002310  216092_s_ .000056  216511_s_ .294776  
216733_s_ .004823  216747_at .002902  216874_at .000117  216969_s_ .001614  
217133_x_ .056851  217198_x_ .169196  217557_s_ .002966  217738_at .000005  
218601_at .023817  218818_at .027554  219302_s_ .000039  219441_s_ .000172  
219574_at .193737  219612_s_ .000075  219697_at .046476  219700_at .003049  
219945_at .000066  219964_at .000684  220234_at .130064  220473_s_ .000017  
220575_at .030223  220633_s_ .058460  220925_at .252465  221256_s_ .721731  
221314_at .002307  221589_s_ .001810  221995_s_ .350859  222071_s_ .000062  
222113_s_ .000023  222208_s_ .100961  222303_at .049265  37226_at  .000749  
60474_at  .000423  
run;

The following statements adjust the p-values in the test data set by using the adaptive adjustments (ADAPTIVEHOLM, ADAPTIVEHOCHBERG, ADAPTIVEFDR, and PFDR), which require an estimate of the number of true null hypotheses ( $\hat{m}_0$ ) or proportion of true null hypotheses ( $\hat{\pi }_0$ ). This example illustrates some of the features and graphics for computing and evaluating these estimates. The NOPVALUE option is specified to suppress the display of the “p-Values” table.

ods graphics on;
proc multtest inpvalues(Probt)=test plots=all seed=518498000
              aholm ahoc afdr pfdr(positive) nopvalue;
   id Probe_Set_ID;
run;
ods graphics off;

Output 61.6.1 lists the requested p-value adjustments, along with the selected value of the “Lambda” tuning parameter and the seed (specified with the SEED= option) used in the bootstrap method of estimating the number of true null hypotheses. The “Lambda Values” table lists the estimated number of true nulls for each value of $\lambda$ , where you can see that the minimum MSE (0.002315) occurs at $\lambda =0.4$ . Output 61.6.2 shows that the SPLINE method failed due to a large slope at $\lambda =0.95$ , so the bootstrap method is used and the MSE plot is displayed.

Output 61.6.1: p and Lambda Values

The Multtest Procedure

P-Value Adjustment Information
P-Value Adjustment	Adaptive Holm
P-Value Adjustment	Adaptive Hochberg
P-Value Adjustment	Adaptive FDR
P-Value Adjustment	pFDR Q-Value
Lambda	0.4
Seed	518498000

Lambda Values
Lambda	MSE	NTrueNull Observed	NTrueNull Spline
0	0.657880	121.000000	67.318707
0.050000	0.030212	43.157895	59.812885
0.100000	0.024897	41.111111	52.636271
0.150000	0.014904	36.470588	46.033846
0.200000	0.008580	32.500000	40.172642
0.250000	0.006476	30.666667	35.157768
0.300000	0.002719	25.714286	31.046105
0.350000	0.002471	24.615385	27.861153
0.400000	0.002378	23.333333	25.595089
0.450000	0.003285	25.454545	24.217908
0.500000	0.003036	24.000000	23.687690
0.550000	0.003567	24.444444	23.965745
0.600000	0.005813	27.500000	25.016579
0.650000	0.004118	22.857143	26.809774
0.700000	0.006647	26.666667	29.321876
0.750000	0.006260	24.000000	32.512203
0.800000	0.013242	30.000000	36.315191
0.850000	0.037624	40.000000	40.618909
0.900000	0.046906	40.000000	45.274355
0.950000	0.332183	80.000000	50.117369

Output 61.6.2: Tuning Parameter Plots

Output 61.6.3 also shows that the bootstrap estimate is used for the PFDR adjustment. The other adjustments have different default methods for estimating the number of true nulls.

Output 61.6.3: Adjustments and Their Default Estimation Method

Estimated Number of True Null Hypotheses
P-Value Adjustment	Method	Estimate	Proportion
Adaptive Holm	Decreased Slope	26	0.21488
Adaptive Hochberg	Decreased Slope	26	0.21488
Adaptive FDR	Lowest Slope	43	0.35537
Positive FDR	Bootstrap	23.3333	0.19284

Output 61.6.4 displays the estimated number of true nulls $\hat{m}_0$ against a uniform probability plot of the unadjusted p-values (if the p-values are distributed uniformly, the points on the plot will all lie on a straight line). According to Schweder and Spjøtvoll (1982) and Hochberg and Benjamini (1990), the points on the left side of the plot should be approximately linear with slope $\frac{1}{m_0+1}$ , so you can use this plot to evaluate whether your estimate of $\hat{m}_0$ seems reasonable.

Output 61.6.4: p-Value Distribution

The NTRUENULL= option provides several methods for estimating the number of true null hypotheses; the following table displays each method and its estimate for this example:

NTRUENULL=	Estimate
BOOTSTRAP	23.3
DECREASESLOPE	26
KSTEST	35
LEASTSQUARES	28
LOWESTSLOPE	43
MEANDIFF	42
SPLINE	50.1

Another method of estimating the number of true null hypotheses fits a finite mixture model (mixing a uniform with a beta) to the distribution of the unadjusted p-values (Allison et al., 2002). Osborne (2006) provides the following PROC NLMIXED statements to fit this model:

proc nlmixed data=test;
   parameters pi0=0.5 a=.1 b=.1;
   pi1= 1-pi0;
   bounds 0 <= pi0 <= 1;
   loglikelihood= log(pi0+pi1*pdf('beta',Probt,a,b));
   model Probt ~ general(loglikelihood);
run;

You might have to change the initial parameter values in the PARAMETERS statement to achieve convergence; see Chapter 64: The NLMIXED Procedure, for more information. This mixture model estimates $\hat{\pi }_0=0$ , meaning that the distribution of p-values is completely specified by a single beta distribution. If the estimate were, say, $\hat{\pi }_0=0.10$ , you could then specify it as follows:

proc multtest inpvalues(Probt)=test ptruenull=0.10
              aholm ahoc afdr pfdr(positive) nopvalue;
   id Probe_Set_ID;
run;

A plot of the unadjusted and adjusted p-values for each test is also produced. Due to the large number of tests and adjustments, the plot is not very informative and is not displayed here.

The top two plots in Output 61.6.5 show how the adjusted values compare with each other and the unadjusted p-values. The PFDR and AFDR adjustments are eventually smaller than the unadjusted p-values since they control the false discovery rate. The adaptive Holm and Hochberg adjustments are almost identical, so the adaptive Holm values are mostly obscured in all four plots. The plot of the Proportion Significant versus the Adjusted p-values tells you how many of the tests are significant for each cutoff, while the plot of the number of false positives (FPN) versus the Proportion Significant tells you how many false positives you can expect for that cutoff.

Output 61.6.5: Adjustment Diagnostics

A Manhattan plot displays –log $_{10}$ of the adjusted p-values, so the most significant tests stand out at the top of the plot. The default plot is not displayed here. The following statements create a Manhattan plot of the adaptive FDR p-values, with the most significant tests labeled with their observation number. The ID values are displayed on the X axis, and the VREF= option specifies the significance level. This plot is typically created with many more p-values, and special ODS Graphics options such as the LABELMAX= option might be required to display the graph. If memory usage is an issue, you might want to store your p-values and use the SGPLOT procedure to create a similar graph.

ods graphics on / labelmax=1000;
proc multtest inpvalues(Probt)=test afdr nopvalue
   plots=Manhattan(label=obs vref=0.0001);
   id Probe_Set_ID;
run;
ods graphics off;

Output 61.6.6: Manhattan Plot

If you have a lot of tests, the “Raw and Adjusted p-Values” and “P-Value Adjustment Diagnostics” plots can be more informative if you suppress some of the tests. In the following statements, the SIGONLY=0.001 option selects tests with adjusted p-values < 0.001 for display. Output 61.6.7 displays tests with their “significant” adjusted p-values:

ods graphics on;
proc multtest inpvalues(Probt)=test plots(sigonly=0.001)=PByTest
              aholm ahoc afdr pfdr(positive) nopvalue;
run;
ods graphics off;

Output 61.6.7: Raw and Adjusted p-Values

The MULTTEST Procedure

Example 61.6 Adaptive Adjustments and ODS Graphics