Разница в результатах GLM между iPython и R

Я пытаюсь справиться с выполнением регрессионных анализов в R. Ниже приведены некоторые случайные фиктивные данные, которые я сгенерировал в R, запустить логический glm в R. Я сохранил данные в тестовом файле, прочитал это в python с ipython (ipython notebook – это потрясающе btw, только начал его использовать!), а затем попытался запустить тот же анализ с помощью python. Результаты очень похожи, но они разные. Я бы ожидал, что они будут такими же. Я сделал что-то неправильно, есть ли параметр, который мне не хватает, или разница из-за какого-то базового расчета?

Любая помощь ценится!

EDIT: Я не знаю, гарантирует ли это закрытие, но я переписываю вещи с измененным (и более чистым) кодом Бена, и результаты между python и R теперь одинаковы. Я вообще не изменил код python, а мой предыдущий код R и код Бена, в то время как разные должны быть (насколько я вижу) делать то же самое. Независимо от того, что вопрос сейчас спорен. Тем не менее, спасибо, что посмотрели.

Генерировать случайные данные и запускать glm:

set.seed(666) dat <- data.frame(a=rnorm(500), b=runif(500), c=as.factor(sample(1:5, 500, replace=TRUE))) library(plyr) dat <- mutate(dat, y0=((jitter(a)^2+(-log10(b)))/(as.numeric(c)/10))+rnorm(500), y=(y0>=mean(y0))) fit1 <- glm(y~a+b+c,data=dat,family=binomial('logit')) summary(fit1) Call: glm(formula = y ~ a + b + c, family = binomial("logit"), data = dat) Deviance Residuals: Min 1Q Median 3Q Max -1.5369 -0.8154 -0.5479 0.9314 2.3831 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.2363 0.3007 4.111 3.94e-05 *** a -0.2051 0.1062 -1.931 0.0535 . b -1.6103 0.3834 -4.200 2.67e-05 *** c2 -0.5114 0.3091 -1.654 0.0980 . c3 -1.3169 0.3147 -4.184 2.86e-05 *** c4 -2.0017 0.3342 -5.990 2.09e-09 *** c5 -2.5084 0.3772 -6.651 2.92e-11 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 631.31 on 499 degrees of freedom Residual deviance: 537.84 on 493 degrees of freedom AIC: 551.84 Number of Fisher Scoring iterations: 4 

вывести те же данные для использования в python:

 write.table(dat, file='test.txt', row.names=F, col.names=T, quote=F, sep=',' ) 

и теперь код python:

 import pandas as pd import statsmodels.api as sm import pylab as pl import numpy as np data = pd.read_csv('test.txt') data.describe() dummy_c = pd.get_dummies(data['c'], prefix='c') data = data[['y','a','b']].join(dummy_c.ix[:,'c_2':]) data_depend = data['y'] data_independ = data.iloc[:,1:] data_independ = sm.add_constant(data_independ, prepend=False) glm_binom = sm.GLM(data_depend, data_independ, family=sm.families.Binomial()) res = glm_binom.fit() print res.summary() 

который дает:

  Generalized Linear Model Regression Results ============================================================================== Dep. Variable: y No. Observations: 500 Model: GLM Df Residuals: 493 Model Family: Binomial Df Model: 6 Link Function: logit Scale: 1.0 Method: IRLS Log-Likelihood: -268.92 Date: Sun, 27 Oct 2013 Deviance: 537.84 Time: 01:26:47 Pearson chi2: 514. No. Iterations: 6 ============================================================================== coef std err t P>|t| [95.0% Conf. Int.] ------------------------------------------------------------------------------ a -0.2051 0.106 -1.931 0.054 -0.413 0.003 b -1.6103 0.383 -4.200 0.000 -2.362 -0.859 c_2 -0.5114 0.309 -1.654 0.098 -1.117 0.094 c_3 -1.3169 0.315 -4.184 0.000 -1.934 -0.700 c_4 -2.0017 0.334 -5.990 0.000 -2.657 -1.347 c_5 -2.5084 0.377 -6.651 0.000 -3.248 -1.769 const 1.2363 0.301 4.111 0.000 0.647 1.826 ============================================================================== 

One Solution collect form web for “Разница в результатах GLM между iPython и R”

Немного очистки кода:

 set.seed(101) dat <- data.frame(a=rnorm(500), b=runif(500), c=as.factor(sample(1:5, 500, replace=TRUE))) library(plyr) dat <- mutate(dat, y0=((jitter(a)^2+(-log10(b)))/(as.numeric(c)/10))+rnorm(500), y=(y0>=mean(y0))) fit1 <- glm(y~a+b+c,data=dat,family=binomial('logit')) fit2 <- update(fit1,control=glm.control(maxit=6)) all.equal(fit1,fit2) coef(fit1) ## (Intercept) ab c2 c3 c4 ## 1.22283193 -0.07544488 -1.54732712 -0.36477556 -1.46313143 -1.95008291 ## c5 ## -3.11914945 

Я согласен с комментарием @ Roland, что поможет воспроизводимый пример. Наиболее вероятным отличием является контрастное кодирование, например:

 fit3 <- update(fit1,contrasts=list(c=contr.sum)) coef(fit3) ## ## (Intercept) ab c1 c2 c3 ## -0.15659594 -0.07544488 -1.54732712 1.37942787 1.01465231 -0.08370356 ## c4 ## -0.57065503 

Если вы используете модель с непрерывными предикторами, лучше ли результаты соответствуют?

Обновление: контрастное кодирование не может быть всей историей, потому что отклонение / логарифмическая вероятность, а также коэффициенты различаются.

  • Используя pandas, как я могу эффективно подбирать большой DataFrame по группам?
  • преобразование списка Python в числовой вектор R
  • условные суммы для совокупности панд
  • Построение векторов в системе координат с R или python
  • Python scipy chisquare возвращает разные значения, чем R chisquare
  • выполнение скрипта R из python
  • Преобразование объектов python для rpy2
  • Unix: Получение Mouse-координаты над X, как Mathematica?
  • Python - лучший язык программирования в мире.