将ANOVA pvalue函数从R迁移到Python

Question

你好

有没有一种有效的方法可以在Python中找到4路ANVOA模型的p值

类似这样的东西在for循环中的R中可以进行一系列模拟

pValues[k] <- anova(lm(Yield ~ Water + Row + Column, data=y))$"Pr(>F)"[1]

我试过researchpy，后来又转到了statsmodels，但我不知道如何从这里开始

pValues[k] = statsmodels.stats.anova_lm(data=y)."Pr(>F)"[1]

Answer 1

在R中：

set.seed(111)
y = data.frame(matrix(rnorm(400),100,4))
colnames(y) = c("Yield","Water","Row","Column")

anova(lm(Yield ~ Water + Row + Column, data=y))
Analysis of Variance Table

Response: Yield
          Df  Sum Sq Mean Sq F value Pr(>F)
Water      1   0.364 0.36410  0.3122 0.5776
Row        1   0.518 0.51768  0.4440 0.5068
Column     1   0.703 0.70256  0.6025 0.4395
Residuals 96 111.942 1.16606           

write.csv(y,"y_data.csv",quote=FALSE,row.names=FALSE)

在python中，您可以使用函数anova_lm从StatsInStatsModels中获取表：

import statsmodels.api as sm
from statsmodels.formula.api import ols
import pandas as pd

y = pd.read_csv("y_data.csv")

mod = ols('Yield ~ Water + Row + Column',data=y).fit()
tab = sm.stats.anova_lm(mod)

            df      sum_sq   mean_sq         F    PR(>F)
Water      1.0    0.364100  0.364100  0.312247  0.577606
Row        1.0    0.517678  0.517678  0.443954  0.506818
Column     1.0    0.702561  0.702561  0.602508  0.439531
Residual  96.0  111.941964  1.166062       NaN       NaN

然后像这样取出p值：

tab["PR(>F)"][0]
Out[8]: 0.5776056586929655