linear-regression - 在 Python 中比较线性模型的对比(比如 Rs 对比库？)

Question

在 R 中，我可以执行以下操作来比较线性模型中的两个对比:

url <- "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extdata/spider_wolff_gorb_2013.csv"
filename <- "spider_wolff_gorb_2013.csv"
install.packages("downloader", repos="http://cran.us.r-project.org")
library(downloader)
if (!file.exists(filename)) download(url, filename)
spider <- read.csv(filename, skip=1)
head(spider, 5)
#   leg type friction
# 1  L1 pull     0.90
# 2  L1 pull     0.91
# 3  L1 pull     0.86
# 4  L1 pull     0.85
# 5  L1 pull     0.80
fit = lm(friction ~ type + leg, data=spider)
fit
# Call:
# lm(formula = friction ~ type + leg, data = spider)
# 
# Coefficients:
# (Intercept)     typepush        legL2        legL3        legL4
#      1.0539      -0.7790       0.1719       0.1605       0.2813
install.packages("contrast", repos="http://cran.us.r-project.org")
library(contrast)
l4vsl2 = contrast(fit, list(leg="L4", type="pull"), list(leg="L2",type="pull"))
l4vsl2
# lm model parameter contrast
# 
#   Contrast       S.E.      Lower     Upper    t  df Pr(>|t|)
#  0.1094167 0.04462392 0.02157158 0.1972618 2.45 277   0.0148

我已经找到了如何在 Python 中完成上述大部分工作:

import pandas as pd
df = pd.read_table("https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extdata/spider_wolff_gorb_2013.csv", sep=",", skiprows=1)
df.head(2)

import statsmodels.formula.api as sm

model1 = sm.ols(formula='friction ~ type + leg', data=df)
fitted1 = model1.fit()
print(fitted1.summary())

现在剩下的就是找到腿对 L4 与腿对 L2 对比的 t 统计量。这在 Python 中可行吗？

Answer 1

statsmodels 仍然缺少一些预定义的对比，但是可以使用模型结果类的 t_test 和 wald_test 或 f_test 方法来测试线性(或仿射)限制。限制由数组或使用参数名称的字符串给出。

关于如何指定对比/限制的细节应该在文档中

例如

>>> tt = fitted1.t_test("leg[T.L4] - leg[T.L2]")
>>> print(tt.summary())
                             Test for Constraints                             
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
c0             0.1094      0.045      2.452      0.015       0.022       0.197
==============================================================================

结果是 t_test 返回的实例中的属性或方法。例如 conf_int 可以通过以下方式获得

>>> tt.conf_int()
array([[ 0.02157158,  0.19726175]])

t_test 被矢量化并将每个限制或对比视为单独的假设。 wald_test 将限制列表视为联合假设:

>>> tt = fitted1.t_test(["leg[T.L3] - leg[T.L2], leg[T.L4] - leg[T.L2]"])
>>> print(tt.summary())
                             Test for Constraints                             
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
c0            -0.0114      0.043     -0.265      0.792      -0.096       0.074
c1             0.1094      0.045      2.452      0.015       0.022       0.197
==============================================================================


>>> tt = fitted1.wald_test(["leg[T.L3] - leg[T.L2], leg[T.L4] - leg[T.L2]"])
>>> print(tt.summary())
<F test: F=array([[ 8.10128575]]), p=0.00038081249480917173, df_denom=277, df_num=2>

此外:如果 cov_type 被指定为 fit 的参数，这也适用于稳健的协方差矩阵。