We are going to fit a model to count data. The first one we will use is about the total count of damage incidents to ships. You can get it here: http://data.princeton.edu/wws509/datasets/ships.dta
We are stealing this dataset from STATA, so it's in STATA format. Pandas has a reader for that format to make your life easy:
from pandas.io.stata import StataReader reader = StataReader('ships.dta') dataframe = reader.data()
Here are some details on this dataset:
The file has 34 rows corresponding to the observed combinations of type of ship, year of construction and period of operation. Each row has information on five variables as follows:
Note that there no ships of type E built in 1960-64, and that ships built in 1970-74 could not have operated in 1960-74. These combinations are omitted from the data file.
Model the damage incident counts with a Poisson Regression.
> Hint: You can look at the previous ipython notebook with the logistic GLM example to see how you can do GLM with statsmodels
Remember that you will have to create dummy variables for categorical variables, and if you have time bins (like "1960-1964"), you have the option of either treating each bin as a category (and create dummies for each bin). Also remember to add a constant (to model the intercept).
Take a look at the statsmodels summary table, the goodness of fit indicators (Deviance, Pearson's chi square statistic) and the coefficients. Is this a good model?
The months of service provides the time interval in which a ship has chances to acquire damages. It can be thought of "exposure", and this column can be used as an offset. You can check these two resources for a quick idea on exposure and using an offset:
Offset in Wikipedia
When to use an offset in CrossValidated (StackOverflow for statistics)
Try your model with months of service as the offset. Does it perform better?
Now separate your data (even though it's only 14 rows) into a training and test set (your test will only be 4 or 5 rows), and check if you predict well (you can look at mean absolute error or mean squared error using
Deviance. Compute the difference in Deviance statistics for your model and the null model. This is called the null deviance. You can do this in one of 2 ways:
Check if this difference is extreme enough that we can reject the null hypothesis. If we can't reject the null hypothesis, we cannot say that this model tells us more than that trivial, null model. To calculate the p-value (prob. of getting a deviance difference at least as extreme as this under the null hypothesis), we need to do a hypothesis test.
You can import the chi square test from scipy for this:
from scipy.stats import chisqprob
Is your model better than the null model?
Now, instead of a poisson regression, do an ordinary least squares regression with log Y. Compare the models. Are the coefficients close? Do they perform similarly?
Now, let's do this on another dataset: Smoking and Cancer.
You can get it here: http://data.princeton.edu/wws509/datasets/smoking.dta This dataset has information on lung cancer deaths by age and smoking status. It has these columns:
That population looks a lot like an offset!
Fit poisson and OLS models and compare them.