Java中的anova.lm()是否有一个等价的函数？

Question

我正在比较R中的两个线性模型和Anova，我想在Java中做同样的事情。为了简化它，我从https://stats.stackexchange.com/questions/48854/why-am-i-getting-different-intercept-values-in-r-and-java-for-simple-linear-regr获取了示例代码，并在下面做了一些修改。模型为test_trait ~ geno_A + geno_B和test_trait ~ geno_A + geno_B + geno_A:geno_B。在R和Java中实现的模型的系数是相同的。在R中，我使用anova(fit, fit2)，其中fits是lm的结果，而在Java中，我使用来自org.apache.commons.math3的TestUtils.oneWayAnovaPValue。

用R得到一个p值的0.797，而用0.817得到一个pvalue，所以这不是正确的方法，但是我无法找到正确的方法。在Java中有相当于R的anova.lm吗？

完整的代码如下。

R

test_trait <- c( -0.48812477 , 0.33458213, -0.52754476, -0.79863471, -0.68544309, -0.12970239,  0.02355622, -0.31890850,0.34725819 , 0.08108851)
geno_A <- c(1, 0, 1, 2, 0, 0, 1, 0, 1, 0)
geno_B <- c(0, 0, 0, 1, 1, 0, 0, 0, 0, 0) 
fit <- lm(test_trait ~ geno_A+geno_B)
fit2 <- lm(test_trait ~ geno_A + geno_B + geno_A:geno_B)

它给出了系数

> fit
Call:
lm(formula = test_trait ~ geno_A + geno_B)

Coefficients:
(Intercept)       geno_A       geno_B  
   -0.03233     -0.10479     -0.60492  

> fit2
Call:
lm(formula = test_trait ~ geno_A + geno_B + geno_A:geno_B)

Coefficients:
  (Intercept)         geno_A         geno_B  geno_A:geno_B  
    -0.008235      -0.152979      -0.677208       0.096383  

和阿诺瓦

> anova(fit, fit2) # 0.797 
Analysis of Variance Table

Model 1: test_trait ~ geno_A + geno_B
Model 2: test_trait ~ geno_A + geno_B + geno_A:geno_B
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1      7 0.77982                           
2      6 0.77053  1 0.0092897 0.0723  0.797

Java

    double [] y =  {-0.48812477,  0.33458213,  
            -0.52754476, -0.79863471,
            -0.68544309, -0.12970239,
             0.02355622, -0.31890850,
             0.34725819,  0.08108851};
double [][] x = {{1,0}, {0,0},
                 {1,0}, {2,1},
                 {0,1}, {0,0},
                 {1,0}, {0,0},
                 {1,0}, {0,0}};
double [][] xb = {{1,0,0}, {0,0,0},
                  {1,0,0}, {2,1,2},
                  {0,1,0}, {0,0,0},
                  {1,0,0}, {0,0,0},
                  {1,0,0}, {0,0,0}};

OLSMultipleLinearRegression regr = new OLSMultipleLinearRegression();
regr.newSampleData(y, x);
double[] beta = regr.estimateRegressionParameters();   

System.out.printf("First model: y = int + genoA + genoB\n");
System.out.printf("Intercept: %.3f\t", beta[0]);
System.out.printf("beta1: %.3f\t", beta[1]);
System.out.printf("beta2: %.3f\n\n", beta[2]);

regr.newSampleData(y, xb);
double[] betab = regr.estimateRegressionParameters();   

System.out.printf("Second model: y = int + genoA + genoB + genoA:genoB\n");
System.out.printf("Intercept: %.3f\t", betab[0]);
System.out.printf("beta1: %.3f\t", betab[1]);
System.out.printf("beta2: %.3f\t", betab[2]);
System.out.printf("beta2: %.3f\n", betab[3]);

它给出与R中相同的系数。

First model: y = int + genoA + genoB
Intercept: -0.032   beta1: -0.105   beta2: -0.605

Second model: y = int + genoA + genoB + genoA:genoB
Intercept: -0.008   beta1: -0.153   beta2: -0.677   beta2: 0.096

但是Anova给出了一个不同的结果

List classes = new ArrayList();
classes.add(beta);
classes.add(betab);
double pvalue = TestUtils.oneWayAnovaPValue(classes);
double fvalue = TestUtils.oneWayAnovaFValue(classes);
System.out.println(pvalue); 
System.out.println(fvalue); 

0.8165390406874127
0.05979444576790511

Answer 1

在比较两种回归的情况下，你非常误解ANOVA。这不是oneWayAnova意义上的变异。R中的onewayAnova等价于函数aov。另一方面，函数anova实现了许多用于比较模型的测试，至少可以说，anova的名称令人困惑.

如果你比较两个回归模型，你想对平方和做一个F检验。在您的代码中，您所做的是一个单向的方差，以查看这两组回归参数是否存在显著差异。这不是您想要做的，但这正是您的JAVA代码所做的。

为了计算正确的F测试，需要执行以下操作：

1. 用残差平方和除以自由度(df)计算最大模型的均方误差(R表: 0.77053 /6)
2. 通过减去两个模型的RSS来计算MSEdifference (结果是“Sq的和”。在R表中，为两个模型减除df (结果是R表中的"Df“)，并对这些数字进行除法。
3. 除以2乘以1，得到F值。
4. 用3中的F值计算p值，对于df，分子中的df差和分母中最大模型的df。

据我所知，类OLSMultipleLinearRegression没有任何方便的方法来提取自由度数，所以这在Java中并不是直接的。您必须手动计算df，然后使用类FDistribution计算p值。

例：

OLSMultipleLinearRegression regr = new OLSMultipleLinearRegression();
regr.newSampleData(y, x);
double SSR1 = regr.calculateResidualSumOfSquares();
double df1 = y.length - (x[0].length + 1); 
    //df = n - number of coefficients, including intercept

regr.newSampleData(y, xb);
double SSR2 = regr.calculateResidualSumOfSquares();
double df2 = y.length - (xb[0].length + 1);

double MSE = SSR2/df2; // EDIT: You need the biggest model here!
double MSEdiff = Math.abs ((SSR2 - SSR1) / (df2 - df1));
double dfdiff = Math.abs(df2 - df1);

double Fval = MSEdiff / MSE;

FDistribution Fdist = new FDistribution(dfdiff, df2);
double pval = 1 - Fdist.cumulativeProbability(Fval);

现在F值和p值都应该正好是R. df1和df2的anova()表中的列Res.Df，R表中的差异应该是Df，MSEdiff应该与Sum of Sq.除以Df和R表。

免责声明:我是一个糟糕的JAVA程序员，所以上面的代码比实际代码更具有概念性。请查找输入错误或愚蠢的错误，并检查我在这里使用的FDistribution类的文档：

https://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/distribution/FDistribution.html#cumulativeProbability%28double%29

现在你知道为什么统计学家使用R而不是Java ;-)

编辑:上面代码中使用的FDistribution是类

org.apache.commons.math3.distribution.FDistribution

FDistribution在JSci中也有：

JSci.maths.statistics.FDistribution

如果使用该代码，则代码的最后一部分如下：

FDistribution Fdist = new FDistribution(dfdiff, df2);
double pval = 1 - Fdist.cumulative(Fval);

根据具体的实现，累积概率可能略有不同。唉，我不知道有什么区别和/或哪个更值得信任。

Answer 2

问题是，您正在比较的方法并不相同。

R中的anova()实际上执行了一个似然比检验，通过添加一个新变量更多的信息在这里的答案来检查您的第二个模型是否有显著的改进

另一方面，java中的oneWayAnovaPValue()只是运行t测试，以检查组之间的均值差异是否显著。在这种情况下，您所做的是比较第一组系数的平均值是否与第二组显著不同，这是不相关的。

据我所知，java中没有现成的函数可以很容易地执行似然比测试。但是你可以很容易地创建一个。在R中，您可以执行以下操作

anova(fit, fit2,test="Chisq")
#p: 0.788

#or manually:
df.diff = fit$df.residual - fit2$df.residual
vals <- (sum(residuals(fit)^2) - sum(residuals(fit2)^2))/sum(residuals(fit2)^2) * fit2$df.residual 
pchisq(vals, df.diff, lower.tail = FALSE)
#p: 0.7879634

因此，您可以在java中采用相同的方法。在google上进行一个简短的搜索，我就可以在java 这里中实现pchisq (请注意，lower.tail=FALSE命令与1-pchisq(lower.tail=TRUE)相同，所以我们并不真正需要这个选项)。

这允许我们执行以下操作

public void regressionRun(){
OLSMultipleLinearRegression regr = new OLSMultipleLinearRegression();
OLSMultipleLinearRegression regr2 = new OLSMultipleLinearRegression();

double[] y = new double[] { -0.48812477, 0.33458213, -0.52754476,
        -0.79863471, -0.68544309, -0.12970239, 0.02355622, -0.31890850,
        0.34725819, 0.08108851 };
double[][] x = new double[10][];
double[][] x2 = new double[10][];
x[0] = new double[] { 1, 0 };
x[1] = new double[] { 0, 0 };
x[2] = new double[] { 1, 0 };
x[3] = new double[] { 2, 1 };
x[4] = new double[] { 0, 1 };
x[5] = new double[] { 0, 0 };
x[6] = new double[] { 1, 0 };
x[7] = new double[] { 0, 0 };
x[8] = new double[] { 1, 0 };
x[9] = new double[] { 0, 0 };
//
x2[0] = new double[] { 1, 0, 0 };
x2[1] = new double[] { 0, 0, 0 };
x2[2] = new double[] { 1, 0, 0 };
x2[3] = new double[] { 2, 1, 2 };
x2[4] = new double[] { 0, 1, 0 };
x2[5] = new double[] { 0, 0, 0 };
x2[6] = new double[] { 1, 0, 0 };
x2[7] = new double[] { 0, 0, 0 };
x2[8] = new double[] { 1, 0, 0 };
x2[9] = new double[] { 0, 0, 0 };

regr.newSampleData(y, x);
double[] b = regr.estimateResiduals();

regr2.newSampleData(y, x2);
double[] b2 = regr2.estimateResiduals();

//calculate sum of squares
double sumsq_b = 0;
double sumsq_b2 = 0;
for (double res : b){
    sumsq_b += res**2;
}
for (double res : b2){
    sumsq_b2 += res**2;
}
//calculate degrees of freedom
int df_b = y.length-(x[0].length+1);
int df_b2 = y.length-(x2[0].length+1);

double vals = (sumsq_b-sumsq_b2)/sumsq_b2*df_b2;

double pvalue = 1-pchisq(vals,df_b-df_b2);
System.out.println(pvalue);
}
//0.7879633810167291