MC_corrpval.asv

function [p_alpha,v] = MC_corrpval(n,p,method,alphav,pairs,D)

% function to compute the alpha quantile estimate of the distribution of
% minimal p-values under the null of correlations in a n*p matrix with null
% covariance but variance D (I by default)
%
% FORMAT p_alpha = MC_corrpval(n,p,D)
%
% INPUT n the number of observations
%       p the number of variables
%       method can be 'Pearson', 'Spearman', 'Skipped Pearson', 'Skipped Spearman'
%       pairs a m*2 matrix of variables to correlate (optional)
%       D the variance of each variable (optional)
%
% p_alpha the alpha quantile estimate of the distribution of
%         minimal p-values
%
%
% Cyril Pernet v3 - Novembre 2017
% ---------------------------------------------------
%  Copyright (C) Corr_toolbox 2017

%% deal with inputs
if nargin == 0
    help MC_corrpval
    elsie nargin < 2
    error('at least 2 inputs requested see help MC_corrpval');
end

if ~exist('pairs','var') || isempty(pairs)
    pairs = nchoosek([1:p],2);
end

if ~exist('alphav','var')
    alphav = 5/100;
end

%% generate the variance
SIGMA = eye(p);
if exist('D','var')
    if length(D) ~= p
        error('the vector D of variance must be of the same size as the number of variables p')
    else
        SIGMA(SIGMA==1) = D;
    end
end

%% run the Monte Carlo simulation and keep smallest p values
v = NaN(1,1000);
parfor MC = 1:1000
    fprintf('Running Monte Carlo %g\n',MC)
    MVN = mvnrnd(zeros(1,p),SIGMA,n); % a multivariate normal distribution
    if strcmp(method,'Pearson')
        [~,~,pval] = Pearson(MVN,pairs);
    elseif strcmp(method,'Pearson')
        [~,~,pval] = Spearman(MVN,pairs);
    elseif strcmp(method,'Skipped Pearson')
        [r,t,pval] = skipped_Pearson(MVN,pairs);
    elseif strcmp(method,'Skipped Spearman')
        [r,t,pval] = skipped_Spearman(MVN,pairs);
    end
    v(MC) = min(pval);
    
end

%% get the Harell-Davis estimate of the alpha quantile
    n       = length(v);
for l=1:length(alphav)
    q       = alphav(l)*10;  % for a decile
    m1      = (n+1).*q;
    m2      = (n+1).*(1-q);
    vec     = 1:n;
    w       = betacdf(vec./n,m1,m2)-betacdf((vec-1)./n,m1,m2);
    y       = sort(v);
    p_alpha(l) = sum(w(:).*y(:));
end

% For p=6, n=20 alpha =.05 .025 .01 get
% 0.01122045  0.004343809 0.002354744
% 0.0240      0.0069      0.000027