Subscribe Bookmark RSS Feed

Noncentrality Parameter Calculation

mewing

Community Trekker

Joined:

Jun 23, 2011

I'm curious if anyone knows and can explain to me how JMP calculates its noncentrality parameter for custom designs.  I know the general form for calculating it, but I can't quite get my math to add up to what JMP is using.  The reason for this is that the d.f. for error used in calculating power for split plot designs is wrong by default in JMP 10.  While I can adjust the d.f. manually it would be nice if I could write a script to do the calculations correctly from the get-go but the design evaluation platform isn't scriptable in JMP 10 (not that I can figure out anyway - even trying a << Get Scriptable Object returned an empty set) so I have to be able to do the math 'by hand' in JSL.

Thanks in advance!

2 REPLIES
bradleyjones

Staff

Joined:

Mar 30, 2012

There is a difference between JMP9 and JMP10 for the power calculations. In JMP9 the power was calculated for a Signal to Noise ratio where the regression coefficient was the signal and the error standard deviation was the noise. In JMP 10 the Signal for a 2-level factor is the change in the response when going from the low value of the factor to the high value of the factor. So now the signal is 2 times the regression coefficient for a 1 degree of freedom effect.

Here is the JSL you need to compute the power calculations in JMP 9 for a 1 degree of freedom effect

Let X be the X matrix for the model.

Then, the variance of the coefficients is

varBetaHat = vecdiag(inverse(X`*X));

alpha = 0.05; // Default significance level

denominatorDF = 11; // df for a 16 run orthogonal design in 4 factors

Fcrit = F Quantile(1-alpha,1,denominatorDF); // critical value of the F distribution

SNRatio = 1;  // default SNRatio

noncentrality = (SNRatio^2)/varBetaHat[2]; // noncentrality for the 1st main effect assuming the 1st effect is the intercept

power = 1 - F Distribution(Fcrit,1,denominatorDF,noncentrality);

Here is the JSL you need to compute the power calculations in JMP 10 for a 1 degree of freedom effect using the same X matrix.

varContrast = 4*vecdiag(inverse(X`*X)); // since the contrast is twice the coefficient, its variance is 4 times larger.

alpha = 0.05; // Default significance level

denominatorDF = 11; // df for a 16 run orthogonal design in 4 factors

Fcrit = F Quantile(1-alpha,1,denominatorDF); // critical value of the F distribution

SNRatio = 2; // default SNRatio is twice as big to accommodate and give the same default value for the power.

noncentrality = (SNRatio^2)/varContrast[2]; // noncentrality for the 1st main effect assuming the 1st effect is the intercept

power = 1 - F Distribution(Fcrit,1,denominatorDF,noncentrality);

mewing

Community Trekker

Joined:

Jun 23, 2011

Dr. Jones,

Thank you for your reply.  Would you mind explaining how this works in the presence of split-plot structure?  I used the custom designer to create a simple 5 factor (2 HTC, 3 ETC factors) and obtained a design with 6 WP and 8 replicates within each WP to make a perfectly replicated design.  I fit the model using REML and then saved the coding table to get the X matrix.  However, I cannot solve X'X with the coding used for the 6 whole plots because they are (intentionally) confounded with the HTC factors.  If I do as you've specified for the ETC factors I get the same ncp as JMP does.  If I collapse the HTC factors down to just the 6 independent settings to treat them as if there were only 6 runs (which is, in effect, what is going on) I get a nice even 6 for the ncp.  However, JMP gives me 5.333 and 4.741 (the exact values are not important).  So, there's something else happening and I'm not sure what it is or how to account for it.

Again, I appreciate you taking the time to provide the above code to get the ncp.