与R不同的Statsmodels-Poisson-glm

Question

我正试图根据R中提供的一些代码来拟合一些模型（空间交互模型）。我已经能够在python框架中使用statsmodels来获得一些代码，但其中一些根本不匹配。我相信我为R和Python编写的代码应该给出相同的结果。有人觉得有什么不同吗？或者说，是否存在一些根本性的差异，这些差异可能会把事情抛在脑后？R代码是与教程中给出的数字相匹配的原始代码（在这里可以找到：http://www.bartlett.ucl.ac.uk/casa/pdf/paper181）。

R样本代码：

library(mosaic)
Data = fetchData('http://dl.dropbox.com/u/8649795/AT_Austria.csv')
Model = glm(Data~Origin+Destination+Dij+offset(log(Offset)), family=poisson(link="log"), data = Data)
cor = cor(Data$Data, Model$fitted, method = "pearson", use = "complete")
rsquared = cor * cor
rsquared

R输出：

> Model = glm(Data~Origin+Destination+Dij+offset(log(Offset)), family=poisson(link="log"), data = Data)
Warning messages:
1: glm.fit: fitted rates numerically 0 occurred 
2: glm.fit: fitted rates numerically 0 occurred 
> cor = cor(Data$Data, Model$fitted, method = "pearson", use = "complete")
> rsquared = cor * cor
> rsquared
[1] 0.9753279

Python代码：

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm
from scipy.stats.stats import pearsonr

Data= pd.DataFrame(pd.read_csv('http://dl.dropbox.com/u/8649795/AT_Austria.csv'))
Model = smf.glm('Data~Origin+Destination+Dij', data=Data, offset=np.log(Data['Offset']), family=sm.families.Poisson(link=sm.families.links.log)).fit()

cor = pearsonr(doubleConstrained.fittedvalues, Data["Data"])[0]

print "R-squared for doubly-constrained model is: " + str(cor*cor)

Python输出：

R-squared for doubly-constrained model is: 0.104758481123

Answer 1

看起来GLM在statsmodels中有收敛问题。也可能在R中，但R只给出这些警告。

Warning messages:
1: glm.fit: fitted rates numerically 0 occurred 
2: glm.fit: fitted rates numerically 0 occurred 

这可能意味着在Logit/Probit上下文中实现完美分离。我得考虑一下泊松模型。

R正在做一个更好的，如果微妙的工作，告诉你，你的装置可能有问题。例如，如果你看statsmodels中的拟合可能性，它是-1.12e27。这应该是一个线索，就在那里，有些东西是关闭的。

直接使用Poisson模型（在可能的情况下，我总是喜欢最大似然法而不是GLM），我可以复制R结果（但我得到了收敛警告）。很明显，同样，默认的牛顿-拉夫森解算器失败了，所以我使用BFG。

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm
from scipy.stats.stats import pearsonr

data= pd.DataFrame(pd.read_csv('http://dl.dropbox.com/u/8649795/AT_Austria.csv'))

mod = smf.poisson('Data~Origin+Destination+Dij', data=data, offset=np.log(data['Offset'])).fit(method='bfgs')

print mod.mle_retvals['converged']