在 Java 中是否有 anova.lm() 的等效函数？
Is there an equivalent function for anova.lm() in Java?
我正在用 Anova 比较 R 中的两个线性模型，我想在 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)，其中拟合是 lm 的结果，在 Java 中，我使用 org.apache.commons.math3.
中的 TestUtils.oneWayAnovaPValue
使用 R 我得到一个 pvalue 0.797，而使用 Java 我得到一个 pvalue 0.817，所以这不是正确的方法，但我找不到如何正确地做。 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
在比较两个回归的情况下，您非常误解方差分析。这不是 oneWayAnova 意义上的方差分析。 onewayAnova 在 R 中等价于函数 aov。另一方面，函数 anova 实现了多种模型比较测试，而名称 anova 至少可以说是令人困惑...
如果比较两个回归模型，您想对平方和进行 F 检验。您在代码中所做的是单向方差分析，以查看两组回归参数是否存在显着差异。那不是你想要做的，但这正是你的 JAVA 代码正在做的。
为了计算正确的 F 检验，您需要执行以下操作：
通过将残差平方和 (RSS) 除以自由度 (df) 计算最大模型的 MSE（在 R table 中：0.77053 / 6
通过减去两个模型的RSS计算MSEdifference（结果在Rtable中为"Sum of Sq."），减去两个模型的df（结果为"Df" 在 R table) 中，然后将这些数字相除。
2除以1得到F值
使用 3 中的 F 值计算 p 值，对于 df，分子中的 df 差和分母中最大模型的 df。
据我所知，class OLSMultipleLinearRegression 没有任何方便的方法来提取自由度数，因此这在 Java 中并不直接。您必须手动计算 df，然后使用 class 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 的 anova() table 中看到的完全相同。df1 和 df2 是列 Res.Df在Rtable中，差值在Rtable中应该是Df，而MSEdiff应该是Sum of Sq.除以[=23] =] 来自 R table.
免责声明：我是一个糟糕的 JAVA 程序员，所以上面的代码比实际代码更具概念性。请寻找拼写错误或愚蠢的错误，并检查我在这里使用的 FDistribution class 的文档：
现在你知道为什么统计学家使用 R 而不是 Java ;-)
编辑：上面代码中使用的 FDistribution 是 class
org.apache.commons.math3.distribution.FDistribution
JSci 中还有一个 FDistribution :
JSci.maths.statistics.FDistribution
如果你使用那个，代码的最后一部分变成：
FDistribution Fdist = new FDistribution(dfdiff, df2);
double pval = 1 - Fdist.cumulative(Fval);
根据具体实施情况，累积概率可能略有不同。唉，我不知道有什么区别 and/or 哪一个更值得信赖。
问题是您比较的方法不一样。
R 中的 anova() 实际上执行 Likelihood-Ratio test to check if your second model significantly improves by adding a new variable: more info in the answer here
另一方面，oneWayAnovaPValue() in java 只是运行 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 here 中实现了 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
我正在用 Anova 比较 R 中的两个线性模型，我想在 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)，其中拟合是 lm 的结果，在 Java 中，我使用 org.apache.commons.math3.
TestUtils.oneWayAnovaPValue
使用 R 我得到一个 pvalue 0.797，而使用 Java 我得到一个 pvalue 0.817，所以这不是正确的方法，但我找不到如何正确地做。 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
在比较两个回归的情况下，您非常误解方差分析。这不是 oneWayAnova 意义上的方差分析。 onewayAnova 在 R 中等价于函数 aov。另一方面，函数 anova 实现了多种模型比较测试，而名称 anova 至少可以说是令人困惑...
如果比较两个回归模型，您想对平方和进行 F 检验。您在代码中所做的是单向方差分析，以查看两组回归参数是否存在显着差异。那不是你想要做的，但这正是你的 JAVA 代码正在做的。
为了计算正确的 F 检验，您需要执行以下操作：
通过将残差平方和 (RSS) 除以自由度 (df) 计算最大模型的 MSE（在 R table 中：0.77053 / 6
通过减去两个模型的RSS计算MSEdifference（结果在Rtable中为"Sum of Sq."），减去两个模型的df（结果为"Df" 在 R table) 中，然后将这些数字相除。
2除以1得到F值
使用 3 中的 F 值计算 p 值，对于 df，分子中的 df 差和分母中最大模型的 df。
据我所知，class OLSMultipleLinearRegression 没有任何方便的方法来提取自由度数，因此这在 Java 中并不直接。您必须手动计算 df，然后使用 class 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 的 anova() table 中看到的完全相同。df1 和 df2 是列 Res.Df在Rtable中，差值在Rtable中应该是Df，而MSEdiff应该是Sum of Sq.除以[=23] =] 来自 R table.
免责声明：我是一个糟糕的 JAVA 程序员，所以上面的代码比实际代码更具概念性。请寻找拼写错误或愚蠢的错误，并检查我在这里使用的 FDistribution class 的文档：
现在你知道为什么统计学家使用 R 而不是 Java ;-)
编辑：上面代码中使用的 FDistribution 是 class
org.apache.commons.math3.distribution.FDistribution
JSci 中还有一个 FDistribution :
JSci.maths.statistics.FDistribution
如果你使用那个，代码的最后一部分变成：
FDistribution Fdist = new FDistribution(dfdiff, df2);
double pval = 1 - Fdist.cumulative(Fval);
根据具体实施情况，累积概率可能略有不同。唉，我不知道有什么区别 and/or 哪一个更值得信赖。
问题是您比较的方法不一样。
R 中的 anova() 实际上执行 Likelihood-Ratio test to check if your second model significantly improves by adding a new variable: more info in the answer here
另一方面，oneWayAnovaPValue() in java 只是运行 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 here 中实现了 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