Sparsity estimation with pplot
FORMAT [nullh,nullhint]=pr_pplot(p,alpha,sortf)
Input:
p - pvalues (column vector)
alpha - probability for confidence intervals
sortf - non-zero if p values are sorted ascending
Output:
nullh - estimated number of null coefficients
The solution is constrained to between 0 and length(p).
nullhint - 100(1-alpha) confidence intervals. These incorporate the
uncertainty due to the Least Squares procedure, not that due
to the selection procedure.
Federico E. Turkheimer, June 15th, 2000
$Id: pr_pplot.m,v 1.1 2005/06/05 04:17:42 matthewbrett Exp $0001 function [nullh,nullhint]=pr_pplot(p,alpha,sortf) 0002 % Sparsity estimation with pplot 0003 % FORMAT [nullh,nullhint]=pr_pplot(p,alpha,sortf) 0004 % 0005 % Input: 0006 % p - pvalues (column vector) 0007 % alpha - probability for confidence intervals 0008 % sortf - non-zero if p values are sorted ascending 0009 % 0010 % Output: 0011 % nullh - estimated number of null coefficients 0012 % The solution is constrained to between 0 and length(p). 0013 % nullhint - 100(1-alpha) confidence intervals. These incorporate the 0014 % uncertainty due to the Least Squares procedure, not that due 0015 % to the selection procedure. 0016 % 0017 % Federico E. Turkheimer, June 15th, 2000 0018 % 0019 % $Id: pr_pplot.m,v 1.1 2005/06/05 04:17:42 matthewbrett Exp $ 0020 0021 if nargin < 1 0022 error('Need p values'); 0023 end 0024 if nargin < 2 0025 alfa = 0.05; 0026 end 0027 0028 [h,k] = size(p); % resizing to column 0029 if(k>1)p=p';end; 0030 n = length(p); 0031 0032 if nargin < 3 0033 sortf = 0; 0034 end 0035 if ~sortf 0036 q = sort(p); 0037 else 0038 q = p; 0039 end 0040 0041 q = flipud(1-q); 0042 0043 % Point-Change analysis 0044 % Reference: 0045 % Stephens MA, "Tests for the Uniform Distribution", 0046 % In Goodness-of-Fit Techniques (D'Agostino RB, Stephens MA Eds.) 0047 % NewYork, Marcel Dekker, pp.331-366, 1986 0048 i = pplot_elbow(q); 0049 0050 % Weights = inverse of the variance of a Beta variable 0051 % Reference: 0052 % Quesenberry CP and Hales C, "Concentration Bands for Uniformity 0053 % Plots", J. Statist. Simul. Comput., 1980, 11:41-53. 0054 0055 j = (1:i)'; 0056 v = (j.*(i-j+1))/((i+2)*(i+1)^2); 0057 w = 1./v(1:i); 0058 0059 % Weighted Linear Regression 0060 x = (1:i)'; 0061 y = q(1:i); 0062 sxx = sum(w.*(x.^2)); 0063 c = sum(w.*x.*y)/sxx; 0064 0065 nu = length(x)-1; % Regression degrees of freedom 0066 yhat = c*x; % Predicted responses at each data point. 0067 r = y-yhat; % Residuals. 0068 rmse = sqrt(sum((r.^2).*w)/(nu*sxx)); % error 0069 tval = spm_invTcdf((1-alpha/2),nu); 0070 0071 0072 nullh = ceil(1/c)-1; 0073 temp = 1/c; 0074 diff = nullh-temp; 0075 if (c<0), 0076 nullh = 0; 0077 nullhint = 0; 0078 end 0079 0080 if (nullh>length(p)), 0081 nullh = length(p); 0082 nullhint = 0; 0083 end 0084 0085 cint = [c-tval*rmse, c+tval*rmse]; 0086 nullhint = [1./cint-1 + diff]; % Correction for rounding