3
$\begingroup$

we are trying to establish numerical equivalence (within reasonable precision) for selected statistical models across programming languages such as SAS, R & Python or even different packages within same programming language. This will allow one to confidently use programming language/Package of their choice, within the parameters used for this testing, while ensuring accurate results. This problem becomes acute for less used advanced statistical models.

So far I haven't noticed any existing work/framework to systematically address this problem.

Our plan so far is to start with creation of simulated test cases with both usual data and edge cases. Thereafter we plan use these to compare results across software for same models.

We want to design & implement a statistically sound & generalizable approach that isn't too expensive computationally.

I would appreciate if someone can point towards any existing work or provide us helpful suggestions in this direction.

Thanks in advance.

$\endgroup$
2
  • 1
    $\begingroup$ I seem to recall that an article in The American Statistician did something similar maybe 20 - 25 years ago. There is an index of this and other ASA journals, maybe covering that period. Sorry I can't locate the article immediately. $\endgroup$
    – BruceET
    Jun 22, 2022 at 17:27
  • 2
    $\begingroup$ I was able to find following publication and its 2nd part, addressing this issue. I am still curious if anyone would have suggestions. McCullough, Bruce D. "Assessing the reliability of statistical software: Part II." The American Statistician 53.2 (1999): 149-159. $\endgroup$
    – Vineet
    Jun 22, 2022 at 18:21

1 Answer 1

4
$\begingroup$

Logically, your exact approach would depend on the audience for your paper, the types of problems they routinely compute, and the specific software to which they have access and know how to use.

Also, your paper may draw attention to seldom-used software that is freely available and very accurate. Or warn against popular software that often makes serious computational errors for particular tasks. So the useful scope of your paper may not be entirely clear until some of your results are known.

If there has not been a paper along these lines since McCullough (1999), you might begin by checking whether any important shortcomings noted there have been overcome in later releases of the software.

You begin your Question with, "[W]e are trying to establish numerical equivalence (within reasonable precision) for selected statistical models across programming languages such as SAS, R & Python or even different packages within same programming language. This will allow one to confidently use programming language/Package of their choice,..." which makes good sense.

Note: I am pretty sure that McCullough (1999) is the paper I vaguely remembered. Glad you found it.

Addenda:

(1) Do R and Minitab give the same P-values for the Welch t test, for (essentially) the same two normal samples.

The Welch two-sample t test is widely used and included in many computer programs. Unlike the pooled 2-sample t test it does not assume that the two samples come from populations with equal variances. For the same data, do two computer programs give essentially the same results.

Consider fictitious data generated in R as follows; R first:

set.seed(1234)
x1 = rnorm(10, 50, 5)
x2 = rnorm(100, 55, 2)

summary(x1);  sd(x1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  38.27   45.94   47.22   48.08   51.96   55.42 
[1] 4.978938
summary(x2);  sd(x2)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  50.64   53.32   54.39   54.73   55.74   60.10 
[1] 1.932823
t.test(x1, x2)$p.val
[1] 0.002190574

Now Minitab:

Now, I input the sample sizes, means, and standard deviations into a recent release of Minitab, which happens to be installed on my computer as I type this (you should use more decimal places of accuracy and the latest release of Minitab). Also, Minitab is one of the few programs that will directly accept such summarized data.]

Sample    N   Mean  StDev  SE Mean
1        10  48.08   4.98      1.6
2       100  54.73   1.93     0.19

Difference = μ (1) - μ (2)
Estimate for difference:  -6.65
95% CI for difference:  (-10.24, -3.06)
T-Test of difference = 0 (vs ≠): 

T-Value = -4.19 P-Value = 0.002 DF = 9

The P-value agrees with the tha P-value from R to two places. (You could also compare the T statistics and DF.)

Difference = μ (1) - μ (2)
Estimate for difference:  -6.65

T-Test of difference = 0 (vs ≠): 
 T-Value = -4.22  P-Value = 0.002  DF = 9

One might also discuss whether it is a good idea to allow the use of summarized data, which cannot be tested for normality (or other assumptions), and Minitab's frequent use of 'sample' to mean observation.

(2) Bootstrapping. There are many styles of bootstrap confidence intervals for parameters and many of them have gained popularity since 2000. Even if two programs are based on the same method, you can't expect exactly the same result because bootstraps use random re-sampling of data. Moreover as this Answer suggests, the best bootstrap CI for a parameter (here the sample variance) may depend on what is known and what is assumed. What would be your standard of judging whether two bootstrap CIs are (essentially) "the same"?

$\endgroup$
2
  • 3
    $\begingroup$ Thanks for taking out time. Our audience are regulatory agencies, statisticians and programmers, who are hesitant to use other software besides the default software (a.k.a SAS). We are brainstorming towards creation of a repository (perhaps open source it in future) via a systematic approach which can be applied on commonly used models in regulatory environment. This can become a reference for statisticians & programmers to freely use application of their choice, instead of sticking with default software. I read McCullough & other relevant papers citing it. This is a good starting point. $\endgroup$
    – Vineet
    Jun 23, 2022 at 17:28
  • $\begingroup$ Computational statistics is a much larger field now than at the turn of the Century. So to get a coherent and useful article you will have to choose topics carefully and keep each within bounds. I will add some suggestions to my Answer and hope that others also have relevant ideas. $\endgroup$
    – BruceET
    Jun 24, 2022 at 22:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.