Subscribe Bookmark RSS Feed

"Make Parms Table" in fit distribution: how can I give this data table a handle to work with?

ursula_garczare

Community Trekker

Joined:

May 1, 2012

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) )

);

2 REPLIES
txnelson

Super User

Joined:

Jun 22, 2012

If you modify your distribution section slightly, you can get the "pt" data table handle specified:

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 )

);

pt= lfdr <<Fit Distribution(Normal Mixtures( Clusters( 2 ) , Make Parms Table ));


you will still have an issue since you are actually going to generate 2 tables.  The good thing is that the pt will be pointed at the last table, and you can get it's name, which will be of the form "Untitled NN", which will allow you to construct it's value:


pt_first = data table(word(1,pt<<get name," ") || " " || char(num(word(2,pt<<get name," "))-1));

Jim
ursula_garczare

Community Trekker

Joined:

May 1, 2012

Thanks! That last-line construct may actually be quite useful for other platforms that generate tables as well, e.g. the response screening platform!