Различия в линейной регрессии в R и Python

Question

Я пытался соответствовать линейные результаты регрессии R с этим питона

Соответствующие коэффициенты для каждой независимой переменной

и ниже код:

Данные загружаются.

https://www.dropbox.com/s/oowe4irm9332s78/X.csv?dl=0

https://www.dropbox.com/s/79scp54unzlbwyk/Y.csv?dl=0

R Код:

#define pathname = " " 

X <- read.csv(file.path(pathname,"X.csv"),stringsAsFactors = F) 
Y <- read.csv(file.path(pathname,"Y.csv"),stringsAsFactors = F) 

d1 = Y$d1 

reg_data <- cbind(d1, X) 
head(reg_data) 

reg_model <- lm(d1 ~ .,data=reg_data) 

names(which(is.na(reg_model$coefficients))) 

reg_model$coefficients 

R Результат

> summary(reg_model) 

Call: 
lm(formula = d1 ~ ., data = reg_data) 

Residuals: 
ALL 60 residuals are 0: no residual degrees of freedom! 

Coefficients: (14 not defined because of singularities) 
      Estimate Std. Error t value Pr(>|t|) 
(Intercept) -67.37752   NA  NA  NA 
v1   1.30214   NA  NA  NA 
v2   -2.93118   NA  NA  NA 
v3   7.58902   NA  NA  NA 
v4   11.88570   NA  NA  NA 
v5   1.60622   NA  NA  NA 
v6   3.71528   NA  NA  NA 
v7   -9.34627   NA  NA  NA 
v8   -3.84694   NA  NA  NA 
v9   -2.51332   NA  NA  NA 
v10   4.22403   NA  NA  NA 
v11   -9.70126   NA  NA  NA 
v12    NA   NA  NA  NA 
v13   4.67276   NA  NA  NA 
v14   -6.57924   NA  NA  NA 
v15   -3.68065   NA  NA  NA 
v16   5.25168   NA  NA  NA 
v17   14.60444   NA  NA  NA 
v18   16.00679   NA  NA  NA 
v19   24.79622   NA  NA  NA 
v20   13.85774   NA  NA  NA 
v21   2.16022   NA  NA  NA 
v22   -36.65361   NA  NA  NA 
v23   2.26554   NA  NA  NA 
v24    NA   NA  NA  NA 
v25    NA   NA  NA  NA 
v26   7.00981   NA  NA  NA 
v27   0.88904   NA  NA  NA 
v28   0.34400   NA  NA  NA 
v29   -5.27597   NA  NA  NA 
v30   5.21034   NA  NA  NA 
v31   6.79640   NA  NA  NA 
v32   2.96346   NA  NA  NA 
v33   -1.52702   NA  NA  NA 
v34   -2.74632   NA  NA  NA 
v35   -2.36952   NA  NA  NA 
v36   -7.76547   NA  NA  NA 
v37   2.19630   NA  NA  NA 
v38   1.63336   NA  NA  NA 
v39   0.69485   NA  NA  NA 
v40   0.37379   NA  NA  NA 
v41   -0.09107   NA  NA  NA 
v42   2.06569   NA  NA  NA 
v43   1.57505   NA  NA  NA 
v44   2.70535   NA  NA  NA 
v45   1.17634   NA  NA  NA 
v46   -10.51141   NA  NA  NA 
v47   -1.15060   NA  NA  NA 
v48   2.87353   NA  NA  NA 
v49   3.37740   NA  NA  NA 
v50   -5.89816   NA  NA  NA 
v51   0.85851   NA  NA  NA 
v52   3.73929   NA  NA  NA 
v53   4.93265   NA  NA  NA 
v54   3.45650   NA  NA  NA 
v55   0.12382   NA  NA  NA 
v56   -0.21171   NA  NA  NA 
v57   4.37199   NA  NA  NA 
v58   3.21456   NA  NA  NA 
v59   0.09012   NA  NA  NA 
v60   -0.85414   NA  NA  NA 
v61   -3.29856   NA  NA  NA 
v62   4.38842   NA  NA  NA 
v63    NA   NA  NA  NA 
v64    NA   NA  NA  NA 
v65    NA   NA  NA  NA 
v66    NA   NA  NA  NA 
v67    NA   NA  NA  NA 
v68    NA   NA  NA  NA 
v69    NA   NA  NA  NA 
v70    NA   NA  NA  NA 
v71    NA   NA  NA  NA 
v72    NA   NA  NA  NA 
v73    NA   NA  NA  NA 

Residual standard error: NaN on 0 degrees of freedom 
Multiple R-squared:  1, Adjusted R-squared: NaN 
F-statistic: NaN on 59 and 0 DF, p-value: NA 

Python код:

Y = pd.read_csv(pathname+"Y.csv") 
X = pd.read_csv(pathname+"X.csv") 

lr = LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False) 
lr.fit(X, Y['d1']) 

(list(zip(lr.coef_, X))) 

lr.intercept_ 

Python Результат:

intercept = 29.396033164254106 

[(-2.4463986167806304, 'v1'), 
(-1.6293010275307021, 'v2'), 
(0.89089949009506508, 'v3'), 
(-3.1021251646895251, 'v4'), 
(-1.7707078771936109, 'v5'), 
(-2.0474705122225636, 'v6'), 
(-1.5537181337496202, 'v7'), 
(-1.6391241229716156, 'v8'), 
(-1.2981646048517046, 'v9'), 
(0.89221826294889328, 'v10'), 
(-0.56694104645951571, 'v11'), 
(2.042810365310288e-14, 'v12'), 
(-2.0312478672439052, 'v13'), 
(-1.5617121392788413, 'v14'), 
(0.4583365939498274, 'v15'), 
(0.8840538748922967, 'v16'), 
(-5.5952681002058871, 'v17'), 
(2.4937042448512892, 'v18'), 
(0.45806845189176543, 'v19'), 
(-1.1648810657830406, 'v20'), 
(-1.7800004329275585, 'v21'), 
(-5.0132817522704816, 'v22'), 
(3.6862778096189266, 'v23'), 
(2.7533531010703882e-14, 'v24'), 
(1.2150003225741557e-14, 'v25'), 
(0.94669823515018103, 'v26'), 
(-0.3082823207975679, 'v27'), 
(0.53619247380957358, 'v28'), 
(-1.1339902793546781, 'v29'), 
(1.9657159583080186, 'v30'), 
(-0.63200501460653324, 'v31'), 
(1.4741013580918978, 'v32'), 
(-2.4448418291953313, 'v33'), 
(-2.0787115960875036, 'v34'), 
(0.22492914212063603, 'v35'), 
(-0.75136276693004922, 'v36'), 
(1.2838658951186761, 'v37'), 
(0.5816277993227944, 'v38'), 
(-0.11270569554555088, 'v39'), 
(-0.13430982360936233, 'v40'), 
(-3.3189296496897662, 'v41'), 
(-0.452575588270415, 'v42'), 
(6.1329755709937519, 'v43'), 
(0.61559185634634817, 'v44'), 
(-1.206315459828555, 'v45'), 
(-3.7452010299772009, 'v46'), 
(-1.1472174665136678, 'v47'), 
(2.8960489381172652, 'v48'), 
(0.0090220136972478659, 'v49'), 
(-5.264918363314754, 'v50'), 
(1.2194758337662015, 'v51'), 
(2.78655271320092, 'v52'), 
(3.106513852668896, 'v53'), 
(3.5181252502607929, 'v54'), 
(-0.34426523770507278, 'v55'), 
(-0.48792823932479878, 'v56'), 
(0.12284460490031779, 'v57'), 
(1.6860388628044991, 'v58'), 
(1.2823067194737174, 'v59'), 
(2.8352263554153665, 'v60'), 
(-1.304336378501032, 'v61'), 
(0.55226132316435139, 'v62'), 
(1.5416988124754771, 'v63'), 
(-0.2605804175310813, 'v64'), 
(1.2489066081702334, 'v65'), 
(-0.44469553013696161, 'v66'), 
(-1.4102990055550157, 'v67'), 
(3.8150423259022639, 'v68'), 
(0.12039684410168072, 'v69'), 
(-1.340699466779357, 'v70'), 
(1.7066389124439212, 'v71'), 
(0.50470752944860442, 'v72'), 
(1.0024872633969766, 'v73')] 

Но это не matching.Please помощь.

Примечание: соответствует для ниже примере

http://davidcoallier.com/blog/linear-regression-from-r-to-python/

Answer 1

TL; др, если вы хотите, чтобы повторить стратегию R в Python вы, вероятно, придется осуществить это самостоятельно, так как R делает некоторые умный материал, который не широко доступен в других местах.

Для справки (так как упоминалось только в комментариях), это некорректное соответствие/ранг-дефицит, который всегда будет иметь место, когда имеется больше переменных-предикторов, чем ответов (p> n: в этом случае p = 73, n = 61), и часто, когда есть много категориальных ответов и/или экспериментальная конструкция каким-то образом ограничена. Работа с этими ситуациями, чтобы получить ответ, который вообще означает что-либо вообще, обычно требует тщательной мысли и/или передовых методов (например, оштрафованная регрессия: см. Ссылки на регрессию лассо и гребня в Wikipedia article on linear regression).

Самый наивный способ справиться с этой ситуацией, чтобы бросить все это в стандартной линейной алгебры и надеюсь, что ничего не ломается слишком плохо, что, по-видимому, что в Python пакет statsmodels делает: от pitfalls document:

• Ранжирование дефектных матриц не приведет к возникновению ошибки.
• Случаи почти идеальной мультиколлинеарности или плохо подготовленных матриц проектирования могут приводить к численным неустойчивым результатам. Пользователям необходимо вручную проверить номер ранга или условия матрицы, если это не является желаемым поведением.

Следующая лучшая вещь (разумно, когда есть небольшой степень коллинеарности) является поворот благоразумно при выполнении линейной алгебры, то есть переставить вычислительную задачу таким образом, чтобы коллинеарных части могут быть оставлены вне. Это то, что делает R; для этого авторам R-кода пришлось модифицировать стандартные подпрограммы LINPACK.

underlying code имеет этот комментарий/объяснение:

с dqrdc2 использует домовладелец преобразование для вычисления QR
гр разложения п р матрица х. ограниченная колонка
c стратегия поворота, основанная на 2-нормах приведенных столбцов
c перемещает столбцы с почти нулевой нормой к правому краю
c матрицей x. эта стратегия означает, что последовательный метод 44 может быть вычислен естественным образом.

c Я очень нервничаю относительно изменения кода linpack таким образом.
c если вы являетесь вычислительным линейным алгеброй алгебры, и вы действительно
c понимаете, как решить эту проблему, пожалуйста, не стесняйтесь
c предложите улучшения этого кода.

Этот код (и комментарий) были в базе кода R начиная с 1998 года; Мне бы очень хотелось узнать, кто изначально написал его, но у меня возникли проблемы после истории кода за пределами реорганизации кода в 1998 году. (Немного больше копания предполагает, что этот код был в кодовой базе R, по существу, с начала его истории , т. е. файл был добавлен в версию SVN 2 1997-09-18 и не был изменен до намного позже.)

Недавно Мартин Мэчлер (2016 25 октября here) добавил дополнительную информацию о ?qr, чтобы эта информация была фактически быть в документации в следующей версии R ...

Если вы знаете, как связать скомпилированный код FORTRAN с кодом Python (я этого не делаю), было бы довольно легко скомпилировать src/appl/dqrdc2.f и перевести начинку lm.fit в Python: это ядро ​​lm.fit, минус проверка ошибок и другая обработка ...

z <- .Call(C_Cdqrls, x, y, tol, FALSE) 
coef <- z$coefficients 
pivot <- z$pivot 
r1 <- seq_len(z$rank) 
dn <- colnames(x) 
nmeffects <- c(dn[pivot[r1]], rep.int("", n - z$rank)) 
r2 <- if (z$rank < p) 
    (z$rank + 1L):p 
else integer() 
if (is.matrix(y)) { ## ... 
} else { 
    coef[r2] <- NA 
    if (z$pivoted) 
     coef[pivot] <- coef 
    names(coef) <- dn 
    names(z$effects) <- nmeffects 
} 
z$coefficients <- coef 

Answer 2

В дополнение к ответу Бен Bolker в.

Основная проблема заключается в том, что статистический пакет должен делать с некорректными проблемами. Насколько я видел, существуют большие различия между пакетами в том, как обрабатываются сингулярные и почти исключительные матрицы дизайна. Единственный детерминированный способ заключается в том, что пользователь явно выбирает алгоритм выбора или наказания.

Если ранг может быть четко идентифицирован, то результат все еще детерминирован, но зависит от выбранной политики сингулярности. Падающие переменные Stata и R, Stata сбрасывает их последовательно, поскольку перечислены переменные, т. Е. Отбрасывает последнюю коллинеарную переменную. Я не знаю, какие переменные R падают. statsmodels использует симметричную обработку переменных и сводит сингулярные значения, используя обобщенный обратный pinv, основанный на сингулярном декомпозиции.Это соответствует крошечному наказанию, как в регрессии регрессии РСС или регрессии Риджа.

Результат не является детерминированным, то есть зависит от пакета линейной алгебры и может даже не совпадать на разных компьютерах, если числовой шум влияет на идентификацию ранга или коллинеарных переменных. В частности, ранжирование, показывающее пороговые значения QR и pinv/SVD, ниже которых ранг идентифицируется как уменьшенный. В statsmodels порог по умолчанию от numpy для относительного номера условия составляет около 1e-15. Таким образом, мы получаем регуляризованное решение, если сингулярные значения меньше этого. Но в некоторых случаях порог может быть слишком низким, а в «неособое» решение преобладают численные шумы, которые невозможно воспроизвести. Я думаю, что то же самое будет верно для любого ранга, показывающего QR или другого чисто численного решения проблемы коллинеарности.
(Robustness issue of statsmodel Linear regression (ols) - Python
https://github.com/statsmodels/statsmodels/issues/2628
и связанный с
statmodels in python package, How exactly duplicated features are handled?)

О ранге выявления, поворачиваясь QR:

Он не был доступен в SciPy, когда большинство statsmodels.OLS было написано, но это было доступно в scipy, как ключевое слово pivoting, в течение некоторого времени, но еще не добавлено в качестве опции для статистических моделей.OLS https://github.com/scipy/scipy/pull/44

Я скептически отношусь к нему как к решению по умолчанию, потому что, я думаю, без его проверки этот численный шум повлияет на то, какие переменные поворачиваются. Тогда больше не будет детерминированным, какие переменные будут отброшены. На мой взгляд, выбор переменных должен быть сознательным выбором пользователя и не оставлять чисто численные варианты.

(отказ от ответственности: Я statsmodels сопровождающим)

редактировать:

Вопрос использует scikit учиться на примере.

Насколько я понимаю, scikit-learn использует другую функцию LAPACK для решения проблемы линейной регрессии, но также основан на декомпозиции сингулярных значений, таких как statsmodels. Однако scikit-learn использует в настоящее время порог по умолчанию для соответствующей функции scipy, которая меньше, чем функция numpy, которую используют statsmodels.
, например. how does sklearn do Linear regression when p >n?

Итак, я ожидаю, что scikit-learn и statsmodels имеют одинаковые результаты в единственном случае, но результаты будут отличаться в некоторых близких случаях.