I try to do a demonstration of false discovery rate procedures. For the local false discoovery rate method, I'd like to make use of the fit distribution with normal mixtures, which means I need to work with the parameters.
I know how to work with make into data table, but I thought, it would be nice to work with the "make parms table" feature of fit distribution and then I get stuck in not knowing how to get access to the values in that table, as I cannot work out how to get a handle to that table.
I am sure there is an easy way....and this is so annoying
p0 = 0.9;
delta = 0.91;// "true" effect size in the normal distributions under the alternative, power 0.8
alpha = 0.2;
nB = 40;
nS = 10000;
n0 = Round( ns * p0 );
n1 = Round( ns / 2 * (1 - p0) );
vec0 = J( n0, 1, Random t( nB - 2 ) );
vec1 = J( n1, 1, Random t( nB - 2, Sqrt( nB ) * delta / 2 ) );
vec2 = J( n1, 1, Random t( nB - 2, Sqrt( nB ) * (-delta) / 2 ) );
allvec = V Concat( vec0, vec1, vec2 );
groupvec = V Concat( J( n0, 1, 1 ), J( n1 * 2, 1, 2 ) );
nt = New Table( "Simulation",
New Column( "Simulated Truth", numeric, nominal, set values( groupvec ) ),
New Column( "t statistic value", numeric, continuous, set values( allvec ) ),
New Column( "p value", numeric, continuous, formula( 2 * (1 - t Distribution( Abs( :t statistic value ), nB - 2 )) ) ),
New Column( "Z value", numeric, continuous, formula( Normal Quantile( p value ) ) )
);
nt << Sort( By( :p value ), Order( Ascending ), replace table );
nt << New Column( "Rank Fraction", numeric, formula( Row() / ns ) );
nt << New Column( "FDR Threshold", numeric, formula( :Rank Fraction * alpha ) );
nt << New Column( "FDR sig", numeric, nominal, formula( :p value < :FDR threshold ) );
nt << New Column( "p05 sig", numeric, nominal, formula( :p value < 0.05 ) );
lfdr = nt << Distribution(
//invisible,
Automatic Recalc( 1 ),
Continuous Distribution( Column( :t statistic value ), Quantiles( 0 ), Summary Statistics( 0 ), Vertical( 0 ) ),
Continuous Distribution( Column( :p value ), Quantiles( 0 ), Summary Statistics( 0 ), Vertical( 0 ) ),
Column( :Z value ),
Quantiles( 0 ),
Summary Statistics( 0 ),
Vertical( 0 ),
Fit Distribution( Normal Mixtures( Clusters( 2 ) ), Make Parms Table )
);
lfdr << title( Concat( "Distributions with p0=", Char( p0, 5, 1 ), " +/- delta=", Char( delta, 5, 1 ), " and n=", Char( nB, 5, 1 ) ) );
pt = The required handle to the data table!!!!!!
//I assume that the lower mu will always be the alternative, as they relate to the lower ps and I hope they are sorted by mu, otherwise this coding is risky
nt << New Column( "H1 density", numeric, continuous, formula( Normal Density( :Z value, pt:Estimate[1], pt:Estimate[3] ) ) );
nt << New Column( "H0 density", numeric, continuous, formula( Normal Density( :Z value, pt:Estimate[2], pt:Estimate[4] ) ) );
nt << New Column( "H1 posterior", numeric, continuous, formula( pt:Estimate[5] * :H1 density ) );
nt << New Column( "H0 posterior", numeric, continuous, formula( pt:Estimate[6] * :H0 density ) );
nt << New Column( "local FDR",
numeric,
continuous,
formula( pt:Estimate[6] * :H0 density / (pt:Estimate[5] * :H1 density + pt:Estimate[6] * :H0 density) )
);