Statsmodels的泊松回归模型与R中不同

Question

Statsmodels的泊松回归模型与R中不同

3

我正在尝试根据提供的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

- user3311076

1

我不确定对于非线性模型来说R^2是否有意义。你可以使用解释偏差与空模型偏差之比来计算类似的指标。或者，你可以探索GLMs的其他替代指标，而不是使用R^2。 - Gavin Simpson

感谢您的回复。根据我所了解的，通常这些模型使用多种不同的度量来检查它们的拟合程度，并且不仅仅依赖于R ^ 2。话虽如此，这两个值相差如此之大肯定有原因。 "Dij"变量有时可以被解释为log（“Dij”），在这种情况下，我已经能够使用相同的代码在R和python + statsmodels中拟合所有感兴趣的模型。 - user3311076

1

你有检查过拟合是否相同吗？比较两个软件中的拟合值。进行一些模型检验。对于困难的数据，有时算法可能效果不佳，需要一些调整等。你似乎从拟合模型转到了检查一些可疑的摘要统计量，而没有检查实际的拟合是否相同。首先，只需查看系数（并确保您知道哪些刻度），Python 中的拟合值是响应或链接函数的比例尺吗？R 的比例尺是响应比例尺。 - Gavin Simpson

因此，我认为两个不同框架的glm函数背后必定存在一些不同的机制，但我搜索了很久也没有找到答案。 - user3311076

1

它们应该代表零值，因此意味着它们实际上几乎为零。这解决了问题...看起来Python处理这些值的方式与R不同，因此替换稍大于1.0^-25的数字将给出预期的模型结果。 - user3311076

显示剩余4条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- jseabold · Accepted Answer

看起来在statsmodels中GLM存在收敛问题。也许在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模型（如果可能，我总是更喜欢最大似然估计），我可以复制R的结果（但我会收到收敛警告）。同样有意思的是，使用默认的牛顿-拉普森求解器失败了，所以我使用bfgs求解器。

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']