Statsmodels Poisson-glm anders als R

Ich versuche, passen einige Modelle (Spatial interaction models) laut einigen code, die in R. ich habe in der Lage zu bekommen, einige der code, um die Arbeit mit statsmodels in einem python-framework, sondern einige von Ihnen überhaupt nicht zusammenpassen. Ich glaube, dass ich den code für R und Python geben sollte identische Ergebnisse. Hat jemand irgendwelche Unterschiede feststellen? Oder gibt es einige grundlegende Unterschiede, die möglicherweise auslösen, die Dinge aus? Der R-code ist der ursprüngliche code entspricht den Nummern in ein tutorial (zu Finden hier: http://www.bartlett.ucl.ac.uk/casa/pdf/paper181).

R Beispielcode:

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 Ausgabe:

> 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-Code:

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 Ausgabe: 

R-squared for doubly-constrained model is: 0.104758481123
  • Ich bin mir nicht sicher, dass der R^2 macht Sinn für ein nicht-lineares Modell. Sie konnte berechnen, etwas ähnliches über das Verhältnis der Abweichung erklären zu null deviance, nehme ich an. Oder gibt es alternativen zu R^2 für GLMs Sie könnten, untersucht.
  • Danke für die Antwort. Von was ich gelesen habe, in der Regel sind diese Modelle verwenden verschiedene Maßnahmen, um zu überprüfen die Passform, und nicht allein verlassen sich auf R^2. Das heißt, es muss einen Grund diese beiden sind so radikal aus. Der "Dij" - variable kann manchmal interpretted als log("Dij") und in diesem Fall war ich in der Lage zu passen alle Modelle von Interesse mit dem gleichen code in beide R und python+statsmodels
  • Gut, haben Sie überprüft, dass die Passungen sind gleich? Vergleichen Sie die Einbauküche Werte in software. Einige model-checking. Mit schwierigen Daten, manchmal auch die algorithmen, die vielleicht nicht gut und kann einige Zureden etc. Sie scheinen zu haben, Weg von der Anpassung eines Modells an die überprüfung einige zweifelhafte Zusammenfassung stat ohne zu prüfen, ob die tatsächlichen passt sind die gleichen. Zum einen, schauen Sie sich nur die Koeffizienten (und stellen Sie sicher, Sie wissen, welche Skalen gibt es auf), und sind ausgestattet Werte in Python über das Ausmaß der Reaktion oder der link-Funktion? R ' s sind auf der Skala der Reaktion.
  • Ja, ich Stimme zu, es wäre nicht hilfreich, um zu springen von der Ausrüstung eines Modells eine einzige zweifelhafte Zusammenfassung stat. Ich weiß, dass die weder die Einbauküche Werte oder die Modell-Koeffizienten sind matching zwischen den R-code und Python-code. Wenn ich den genau gleichen code, aber das Modell ändern der Formel auf "Daten~Ursprung+Ziel+log(Dij)" bekomme ich eine passende R-squared-Werte, passenden ausgestattet Werte und matching-Koeffizienten zwischen den beiden frameworks. Wenn ich"Dij" als eine variable, die nicht unter das Protokoll, nichts passt, obwohl ich schon unter der R-code als richtig, denn es wurde von einem Beispiel in einem working paper
  • Ich bin daher denken, es muss einen Unterschied mechanik hinter der glm-Funktionen zwischen den beiden anderen frameworks, die ich gesucht und nicht gefunden.
  • Das deutet auf ein problem mit der Befestigung. Überprüfung der Konvergenz-Kriterien in der R-Modell, plot der angepassten Modell, zu betrachten, die passen. Die Ausgabe, die Sie zeigen, für R schlägt vor, einige Fragen - es ist die Warnung, dass Sie passend ein Modell, wo \hat{E(mu)} vorhergesagt werden 0, D. H. der Erwartung der Mittelwert der Poisson-Verteilung ausgestattet, wie 0. Das deutet auf eine mögliche Problem mit der Passform, da einige ungerade Werte in Dij wird reduziert, wenn Sie zu transformieren, um die log-Skala.
  • Sie haben mehrere 1.000000e-300 in Ihrem Dij. Sind Sie sein sollen fehlende Werte oder wirklich fast null?
  • Auch jede zehnte Data null ist. So wie Gavin wies darauf hin, es ist etwas seltsam mit den Daten. Ich weiß nicht, was R tut, aber die Ergebnisse in statsmodels Aussehen stark beeinflusst durch diese Ausreißer. Wenn ich drop diesen Beobachtungen, dann bekomme ich mit 'Dij' als regressor eine rsquared von 0.97486 mit statsmodels (mit 72 Beobachtungen).
  • Sie repräsentieren die Werte von null und damit gemeint sind wirklich fast null. Diese behoben das Problem...es wird angezeigt, dass python die Verarbeitung dieser Werte anders als R und daher das in zahlen etwas größer 1.0^-25 sagen wird geben dem Modell erwarteten Ergebnisse.

InformationsquelleAutor user3311076 | 2014-02-14
  1. 3

    Sieht es aus wie GLM Konvergenz hat hier Probleme in statsmodels. Vielleicht ist in R auch, aber R gibt nur diese Warnungen.

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

    Das konnte bedeuten, dass etwas wie eine perfekte Trennung in Logit - /Probit-Kontext. Ich würde darüber nachdenken, Sie für eine Poisson-Modell.

    R ist besser, wenn subtile, job, Ihnen zu sagen, dass etwas falsch sein kann, die in Ihrer Armatur. Wenn Sie einen Blick auf die Einbauküche Wahrscheinlichkeit in statsmodels zum Beispiel, ist es -1.12e27. Das sollte ein Hinweis darauf sein Recht gibt, dass etwas ausgeschaltet ist.

    Verwendung von Poisson-Modell direkt (ich bevorzuge immer die maximale Wahrscheinlichkeit zu GLM wenn möglich), kann ich replizieren der R Ergebnisse (aber ich bekomme eine Konvergenz Warnung). Sich vielsagend wieder das Standard-newton-raphson solver versagt, so dass ich mit 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']
    InformationsquelleAutor jseabold
Schreibe einen Kommentar