2014-02-14 3 views
3

Я пытаюсь подгонять некоторые модели (модели пространственного взаимодействия) в соответствии с некоторым кодом, который предоставляется в R. Мне удалось заставить часть кода работать с использованием statsmodels в рамках python, но некоторые из них совсем не совпадают. Я считаю, что код, который у меня есть для R и Python, должен давать одинаковые результаты. Кто-нибудь видит какие-то различия? Или есть некоторые фундаментальные различия, которые могут отбросить все это? Код R - это исходный код, который соответствует номерам, указанным в учебнике (см. Здесь: http://www.bartlett.ucl.ac.uk/casa/pdf/paper181).Statsmodels Poisson glm отличается от R

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 
+1

Я не уверен, что R^2 имеет смысл для нелинейной модели. Я мог бы вычислить что-то подобное, используя отношение отклонения объяснения к нулевому отклонению. Или есть альтернативы R^2 для GLM, которые вы могли исследовать. –

+0

Спасибо за ответ. Из того, что я прочитал, обычно эти модели используют несколько различных мер для проверки их соответствия и не только полагаются на R^2. Тем не менее, должна быть причина, по которой эти двое так радикальны. Переменная «Dij» иногда может быть интерпретирована как log («Dij»), и в этом случае я смог сопоставить все интересующие модели с использованием того же кода в моделях R и python + statsmodels – user3311076

+1

Ну, вы проверили, что подходит одинаковы? Сравните установленные значения в обоих программах. Сделайте некоторые проверки модели. Со сложными данными иногда алгоритмы могут не работать хорошо, и может потребоваться некоторое уговоры и т. Д. Кажется, вы отказались от установки модели для проверки некоторого сомнительного сводного стата, не проверяя, совпадают ли фактические соответствия. Во-первых, просто посмотрите на коэффициенты (и убедитесь, что знаете, какие масштабы там есть), и являются ли установленные значения в Python в масштабе ответа или функции ссылки? R находятся в масштабе ответа. –

ответ

3

Похоже, что 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 (но я получаю предупреждение о сближении). Точно, опять же, решатель по умолчанию newton-raphson терпит неудачу, поэтому я использую 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'] 
+1

Я подал отчет об ошибке, чтобы посмотреть, что именно происходит здесь, в statsmodels GLM. https://github.com/statsmodels/statsmodels/issues/1391 – jseabold

Смежные вопросы