<jupyter_start><jupyter_text>MIT DATATHHON CHALLENGE COVID<jupyter_code>import pandas as pd
import numpy  as np
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
pd.set_option('display.max_columns', None)
final_dataset  = pd.read_csv('final_dataset.csv')<jupyter_output><empty_output><jupyter_text>Regression goalOur regression goal is to predict the number of cases on any given day.y = the vector of cases counts seen on days 1 through n. Thus y = [y_1, y_2, y_3,…,y_n]. y_i is the number of bicyclists on day i.X = the matrix of predictors a.k.a. regressors a.k.a explanatory variables a.k.a. regression variables. The size of matrix X is a (n x m) since there are n independent observations (rows) in the data set and each row contains values of m explanatory variables.λ = the vector of event rates. The vector λ is a primary characteristic of count based data sets. λ is a vector of size (n x 1). It contains n rates [λ_0, λ_1, λ_2,…,λ_n], corresponding to the n observed counts in the counts vector y. The rate λ_i for observation ‘i’ is assumed to drive the actual observed count y_i in the counts vector y. The λ column is not present in the input data.Instead, λ vector is a deduced variable that is calculated by the regression model during the training phase.For the bicyclist counts data, each one of the λ_i values is defined as the number of cases in ‘unit’ time on day i. Unit time can be 1 second, 1 hour, 1 day, 1 week — whatever unit time interval we want to measure the rate over. This rate λ_i is assumed to drive the observed count of bicyclists y_i on day i. The training algorithm of the Negative Binomial regression model will fit the observed counts y to the regression matrix X. Recollect that the Negative Binomial regression model does not make the variance = mean assumption that the Poisson regression model does. NB model requires us to define a new parameter α which it uses to express the variance in terms of the mean as follows- When p = 1: - When p=2: Then The p=2 case is referred to as the NB2 model.<jupyter_code>final_dataset.head()
final_dataset.tail()
#final_dataset.info()
final_dataset.date.unique()
final_dataset['date']= pd.to_datetime(final_dataset['date'], errors='coerce') 
#final_dataset['date'] = final_dataset['date'].astype('datetime64[ns]')
final_dataset.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 65 entries, 0 to 64
Data columns (total 26 columns):
 #   Column                  Non-Null Count  Dtype         
---  ------                  --------------  -----         
 0   countyfips              65 non-null     int64         
 1   date                    62 non-null     datetime64[ns]
 2   super_spreader          65 non-null     int64         
 3   total_businesses        65 non-null     int64         
 4   county                  65 non-null     object        
 5   state                   65 non-null     object        
 6   cases                   65 non-null     int64         
 7   deaths                  65 non-null     int64         
 8   county_name             65 non-null     object        
 9   population              65 non-null     int64         
 10  households              65 non-null     int64         
 11  age                     65 non-null     float64       
 12  income                  65 non-null     int64       [...]<jupyter_text>STEP1: fitting the Poisson regression mode<jupyter_code>final_dataset['MONTH'] = final_dataset['date'].dt.month
final_dataset['DAY_OF_WEEK'] = final_dataset['date'].dt.dayofweek
final_dataset['DAY'] = final_dataset['date'].dt.day
final_dataset['YEAR'] = final_dataset['date'].dt.year
final_dataset.columns
final_dataset.head(2)
final_dataset.head(2)
mask = np.random.rand(len(final_dataset)) < 0.8
final_dataset_train = final_dataset[mask]
final_dataset_test = final_dataset[~mask]
print('Training data set length='+str(len(final_dataset_train)))
print('Testing data set length='+str(len(final_dataset_test)))
final_dataset_test.head(2)
from patsy import dmatrices, dmatrix, demo_data
expr = """daily_cases ~ DAY + DAY_OF_WEEK + MONTH + population +  households + age + income +  females + over65 + white +  black  + asian + hispanic + 

amerindian +  other +  poverty +  no_hsdiploma + super_spreader_density """
y_train, X_train = dmatrices(expr, final_dataset_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, final_dataset_test, return_type='dataframe')

poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()
print(poisson_training_results.summary())
##  λ vector  
print(poisson_training_results.mu)
print(len(poisson_training_results.mu))
print(type(poisson_training_results.mu))<jupyter_output>[  23.79620265   13.69147513   42.9582616    16.55406062 1827.06986428
 1294.61475184 2426.29128654 2910.4755498  2070.5485471   277.84822158
  592.65054918 1136.76899212 1343.73223699   95.94505867   70.70768599
  110.34725534  148.67984393  249.24030811  150.07984795   82.91921615
  160.97401983  117.10676401 1104.64366785  831.47587756 1640.26942038
 1845.14087652 1293.47015744   70.7247974   145.13707872  150.9825926
  111.15553127   28.68748719   37.04337465   27.26913815  889.12607452
  677.73549136 1267.71945286 1423.1031027  1008.31587841   77.68222444
   56.20141914  123.62910071   61.48725571   49.8366921    35.77431824
   86.88400715   56.50498251    7.89025086   19.13071406   32.77044849
   40.20858659]
51
<class 'numpy.ndarray'><jupyter_text>STEP 2:  Auxiliary OLS regression model on the data set and use the fitted model to get the value of α.<jupyter_code>import statsmodels.formula.api as smf
poisson_training_results_mu = np.array(poisson_training_results.mu).ravel().tolist()

poisson_training_results_mu  = pd.Series(poisson_training_results_mu)
#poisson_training_results_mu
#type(poisson_training_results_mu)
## Add the λ vector as a new column 

final_dataset_train['daily_cases_LAMBDA'] = poisson_training_results_mu
final_dataset_train.shape
final_dataset_train.head(2)
final_dataset_train['AUX_OLS_DEP'] = final_dataset_train.apply(lambda x: ((x['daily_cases'] - x['daily_cases_LAMBDA'])**2 - x['daily_cases']) / x['daily_cases_LAMBDA'],
                                                               axis=1)
final_dataset_train.head()
#final_dataset_train
ols_expr = """AUX_OLS_DEP ~ daily_cases_LAMBDA - 1"""
aux_olsr_results = smf.ols(ols_expr, final_dataset_train).fit()
# The determination of α
print(aux_olsr_results.params)
aux_olsr_results.tvalues<jupyter_output><empty_output><jupyter_text>STEP 3: Provide the value of alpha found in STEP 2 into the statsmodels.genmod, and train the NB2 model on the training data set.<jupyter_code>nb2_training_results = sm.GLM(y_train, X_train,family=sm.families.NegativeBinomial(alpha=aux_olsr_results.params[0])).fit()
print(nb2_training_results.summary())<jupyter_output>Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:            daily_cases   No. Observations:                   51
Model:                            GLM   Df Residuals:                       34
Model Family:        NegativeBinomial   Df Model:                           16
Link Function:                    log   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -316.31
Date:                Fri, 15 May 2020   Deviance:                       9.6167
Time:                        17:46:55   Pearson chi2:                     7.75
No. Iterations:                     9                                         
Covariance Type:            nonrobust                                         
==========================================================================================
                             coef    std[...]<jupyter_text>STEP 4: Predictions using our trained NB2 model.<jupyter_code>nb2_predictions = nb2_training_results.get_prediction(X_test)
predictions_summary_frame = nb2_predictions.summary_frame()
print(predictions_summary_frame)
predicted_counts= predictions_summary_frame['mean']
actual_counts = y_test['daily_cases']
fig = plt.figure()
fig.suptitle('Predicted_cases versus actual daily_cases COVID-19')
predicted, = plt.plot(X_test.index, predicted_counts, 'go-', label='Predicted_cases counts')
actual, = plt.plot(X_test.index, actual_counts, 'ro-', label='Actual counts')
plt.legend(handles=[predicted, actual])
fig.set_size_inches(15,10)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1013277428371768, 'fraction_numerical': 0.08176100628930817, 'mean_word_length': 3.1184652278177456, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '364075', 'license': 'MIT', 'n_tokens_mistral': 3110, 'n_tokens_neox': 2611, 'n_words': 746, 'path': 'Negative_Binomial_Regression.ipynb', 'repo_name': 'garbamoussa/Negative-Binomial-Regression-', 'size': '8586'}