私は R で回帰分析を実行する方法を理解しようとしています。以下は、R で生成したランダムなダミー データです。R でロジスティック glm を実行します。データをテスト ファイルに保存し、それを Python に読み込みます。 ipython (ipython ノートブックは素晴らしいところで、使い始めたばかりです!)、次に python で同じ分析を実行しようとしました。結果は非常に似ていますが、異なります。私は彼らが同じであることを期待していたでしょう。私は何か間違ったことをしましたか、欠落しているパラメータがありますか、それとも基礎となる計算による違いはありますか?
どんな助けでも大歓迎です！
編集:これで閉じる必要があるかどうかはわかりませんが、ベンの編集された (そしてよりクリーンな) コードで再実行したところ、python と R の間の結果は同じになりました。私は python コードをまったく変更しませんでした。私の以前の R コードと Ben のコードは、(私が見る限り) 同じことをしているはずです。とにかく、質問のポイントは今では議論の余地があります。それにもかかわらず、ご覧いただきありがとうございます。
ランダム データを生成し、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
==============================================================================