BESTDw.p format in proc iml

[greveam]

Hi, I want to format p in the macro below so that I can report the exact p-value. I tried using the bestd32.31 format but still get a p rounded as 0. Any help would be much appreciated. Thanks!

%macro combchi(df=,chi=);
proc iml;
df=&df;
g2={&chi};
m=ncol(g2);
g=sqrt(g2);
mg2=sum(g2)/m;
r=(1+1/m)*(ssq(g)-(sum(g)**2)/m)/(m-1);
f=(mg2/df - r*(m-1)/(m+1))/(1+r);
ddf=(m-1)*(1+1/r)**2/df**(3/m);
p=1-probf(f,df,ddf);
print f df ddf;
print p[format=bestd32.31];
quit;
%mend combchi;

%combchi(df=2,chi=102.997006 108.592438 105.224443 111.206028 106.644927)

SAS Output:

The SAS System

f df ddf

52.009228

2

3641.6389

p

0

[Rick_SAS · Accepted solution]

The problem is not the format, it is the way you are computing the p-value.  When you compute the quantity 

p = 1 - probf(...);

the PROBF function computes the left tail of the F distribution and returns 1 to within numerical precision. When you want an extreme p-value, it is more accurate numerically to compute the RIGHT tail directly by using the SDF function:

p = sdf("F", f,df,ddf);

[ballardw]

You might also want to consider using the PVALUEw.d format.

---

@greveam wrote:

Hi, I want to format p in the macro below so that I can report the exact p-value. I tried using the bestd32.31 format but still get a p rounded as 0. Any help would be much appreciated. Thanks!

31 decimal places is still likely not "the exact p-value". It is the value displayed to a precision of 31 decimals. Though when fewer digits are shown it could be the result of another rounding/truncation/precision issue in the calculations involved.