Java中的anova.lm()是否有等效函数？

2022-12-22 11:32:17

我将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),其中拟合是lm的结果,在Java中,我使用org.apache.commons.math3中的TestUtils.oneWayAnovaPValue.

使用R我得到一个0.797的pvalue,而使用Java我得到一个0.817的pvalue,所以这不是正确的方法,但我找不到如何正确地做到这一点.在Java中是否有相当于R的anova.lm？

完整代码如下.

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

> 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

解决方法:

在比较两个回归的情况下,你非常误解ANOVA.这不是oneWayAnova意义上的ANOVA. R中的onewayAnova相当于函数aov.另一方面,anova功能实现了大量的模型比较测试,而anova的名称至少令人困惑……

如果比较两个回归模型,则需要对平方和进行F检验.您在代码中执行的操作是单向ANOVA,以查看两组回归参数是否存在显着差异.这不是你想要做的,但这正是你的JAVA代码所做的.

要计算正确的F测试,您需要执行以下操作：

>通过将剩余平方和(RSS)除以*度(df)来计算最大模型的MSE(在R表中：0.77053 / 6
>通过减去两个模型的RSS来计算MSE差异(结果是R表中的“Sum of Sq.”),减去两个模型的df(结果是R表中的“Df”),并除以这些数字.
>将2除以1,得到F值
>使用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的anova()表中看到的完全相同,而df2是R表中的Res.Df列,差异应该是R表中的Df,以及MSEdiff应该与Sq的总和相同.除以R表中的Df.

免责声明：我是一个糟糕的JAVA程序员,因此上面的代码比实际代码更具概念性.请查找拼写错误或愚蠢的错误,并查看我在此处使用的FDistribution类的文档：

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

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

码农公寓

相关文章