with open('un_explore.pkl') as picklefile:
    undata = pickle.load(picklefile)

"""
Explanation: Open previously saved pickle file
End of explanation
"""


for col in sorted(undata.columns):
    print undata[col].count(), col

"""
Explanation: Print column names and counts for non-NaN values for each column
End of explanation
"""


contr_use = ['condom_use_at_last_highrisk_sex_1524_years_old_men_percentage',
    'condom_use_at_last_highrisk_sex_1524_years_old_women_percentage',
    'condom_use_population_ages_1524_female__of_females_ages_1524',
    'condom_use_population_ages_1524_male__of_males_ages_1524',
    'condom_use_to_overall_contraceptive_use_among_currently_married_women_1549_years_old_percentage',
    'condom_use_with_non_regular_partner__adults1549_female',
    'condom_use_with_non_regular_partner__adults1549_male',
    'contraceptive_prevalence__of_women_ages_1549']

for col in contr_use:
    print undata[col].count(), col

"""
Explanation: Print counts specifically for contraceptive use variables
End of explanation
"""


def model_lr_data(dep, ind, data):
    formula = dep + '~' + '+'.join(ind)
    y, X = dmatrices(formula, data=data, return_type='dataframe' )
    print X.shape
    print '\n****PREDICTING\n' + dep
    model = sm.OLS(y, X)
    results = model.fit()
    print results.summary()

hiv = 'hiv_incidence_rate_1549_years_old_percentage_midpoint'
pop = 'population_total'
oda = 'net_oda_mill'
gni = 'gni_per_capita_atlas_method_current_us'
prev = 'prevalence_of_hiv_total__of_population_ages_1549'
contr = 'condom_use_at_last_highrisk_sex_1524_years_old_women_percentage'
pvty = 'population_below_national_poverty_line_total_percentage'
abr = 'adolescent_birth_rate_per_1000_women'

"""
Explanation: Create funtion for linear regression modeling
The model_lr_data function creates matrices and a linear regression model and outputs a summary of the results. The inputs are:

dep: dependent variable
ind: a list of features / independent variables
data: which dataframe to use
End of explanation
"""


features1 = [
            'net_oda_mill', 
            pop,
            gni,
            prev,
            contr,
            pvty,
            abr,
            'year', 
            'mdgregions',
            'islldc',
            'isldc2014',
            'isdeveloped'
           ]

model_lr_data(hiv, features1, undata)

"""
Explanation: Linear model 1
Use a number of features which may be important to predict HIV incidence rate worldwide.
End of explanation
"""


features2 = [
            'net_oda_mill', 
            gni,
            prev,
            abr,
            'year', 
            'mdgregions',
            'isdeveloped'
           ]

model_lr_data(hiv, features2, undata)

"""
Explanation: The first model has a high R-squared, but is only using 15 observations, since many of the selected features include rows with missing values.

Linear model 2
Limit the features so that more observations will be included. Predict HIV incidence rate worldwide.
End of explanation
"""


ssafrica = undata[undata['mdgregions'] == 'Sub-Saharan Africa']

ssafrica.shape

features3 = [
            'net_oda_mill', 
            gni,
            prev,
            'year'
           ]

model_lr_data(hiv, features3, ssafrica)

"""
Explanation: Contraceptive data is too sparse to make accurate predictions. Prevelance of HIV appears significant, GNI appears to be significant, adolescent birth rate is close to significant. ODA is not significant.

Linear model 3
Sub-Saharan Africa is the area of greatest concern. Create a subset of sub-Saharan Africa data for building models. Use limited features to predict HIV incidence rate in sub-Saharan Africa.
End of explanation
"""


keys = ['countryname', 'iso3code', 'year', 'mdgregions', 'islldc', 'isldc2014', 'ismdgcountry']

df1 = undata[keys + [hiv]]
df1 = df1.set_index(keys)

df1 = df1.sort()
df2 = df1.pct_change()
df2 = df2.reset_index()
df2.columns = keys + ['hiv_pctchange']
df2.tail()

undata.shape

undata2 = pd.merge(undata, df2, how='outer', on=keys)
undata2.shape

"""
Explanation: ODA received and GNI do not appear to be significant. Preveleance of HIV and year are the most significant.

Create new column - HIV percent change
Create new column to use as dependent variable, hiv_pctchange, which indicates the percent change in the HIV incidence rate from the previous year.
End of explanation
"""


ssafrica2 = undata2[undata2['mdgregions'] == 'Sub-Saharan Africa']

features4 = [
            'net_oda_mill', 
            gni,
            prev,
            abr,
            contr,
            'year', 
            'isdeveloped'
           ]

model_lr_data('hiv_pctchange', features4, undata2)

"""
Explanation: Linear model 4
Predict on new variable, percent change of HIV incidence rate, in sub-Saharan Africa.
End of explanation
"""


undata3 = copy.deepcopy(undata2[np.isfinite(undata2['hiv_pctchange'])])
undata3.shape

"""
Explanation: R-squared is very low at 0.17, no features appear to be significant, and the model is trained on very few observations.

Classification models
Convert to a classification problem and investigate feature importance.
Is HIV incidence rate decreasing from previous year?
End of explanation
"""


isdecrease = []

for change in undata3['hiv_pctchange']:
    if change >= 0:
        isdecrease.append(0)
    if change < 0:
        isdecrease.append(1)
        
len(isdecrease)

undata3['isdecrease'] = isdecrease

"""
Explanation: Create new variable isdecrease.
End of explanation
"""


rs = 6

models = {'logistic': LogisticRegression(),
          'naive bayes': GaussianNB(),
          'SVM': SVC(random_state=rs, probability=True),
          'gradient boosting': GradientBoostingClassifier(n_estimators=250, random_state=rs),
          'decision tree': DecisionTreeClassifier(random_state=rs),
          'random forest': RandomForestClassifier(n_estimators=250, random_state=rs),
          'extra trees': ExtraTreesClassifier(n_estimators=250, random_state=rs)
         }

def process_data(dep, ind, data):
    formula = dep + '~' + '+'.join(ind)
    y, X = dmatrices(formula, data=data, return_type='dataframe')
    X = X.iloc[:,1:]
    y = y.iloc[:,0]
    X_tr, X_ts, y_tr, y_ts = train_test_split(X, y, test_size = 0.25, random_state=rs)
    return X_tr, X_ts, y_tr, y_ts

def get_scores(model_dict):
    for mname, m in model_dict.iteritems():
        print "*** %s" % mname
        m.fit(X_tr, y_tr)
        preds = m.predict(X_ts)
        proba = m.predict_proba(X_ts)
        print 'accuracy: %f' % accuracy_score(y_ts, preds)
        print 'precision: %f' % precision_score(y_ts, preds)
        print 'recall: %f' % recall_score(y_ts, preds)
        print 'f1 score: %f' % f1_score(y_ts, preds)
        print '\n'
        all_preds[mname] = preds
        all_proba[mname] = proba

def get_crossval_scores(X, y, model_dict):
    print 'CROSS VALIDATION SCORES'
    for mname, m in model_dict.iteritems():
        print '\n*** %s' % mname
        acc = np.mean(cross_val_score(m, X, y, scoring='accuracy'))
        pre = np.mean(cross_val_score(m, X, y, scoring='precision'))
        rec = np.mean(cross_val_score(m, X, y, scoring='recall'))
        f1 = np.mean(cross_val_score(m, X, y, scoring='f1'))
        print 'cv score: %f' % np.mean(cross_val_score(m, X, y))
        print 'accuracy: %f' % acc
        print 'precision: %f' % pre
        print 'recall: %f' % rec
        print 'f1 score: %f' % f1

ssafrica3 = undata3[undata3['mdgregions'] == 'Sub-Saharan Africa']

features5 = [
            'net_oda_mill', 
            gni,
            prev,
            'isdeveloped'
           ]

X_tr, X_ts, y_tr, y_ts = process_data('isdecrease', features5, ssafrica3)

X_tr.shape, X_ts.shape

all_preds = {}
all_proba = {}

get_scores(models)

get_crossval_scores(X_tr, y_tr, models)

"""
Explanation: Classification
Create functions to run multiple classification models at once and output metrics to assess accuracy.
End of explanation
"""


trees = ExtraTreesClassifier(random_state=rs)

trees.fit(X_tr, y_tr)
importances = trees.feature_importances_

indices = np.argsort(importances)[::-1]

print("Feature ranking:")

for f in range(len(features5) - 1):
    print("%d. feature: %s (%f)" % (f + 1, features5[indices[f]], importances[indices[f]]))

forest = RandomForestClassifier(random_state=rs)

forest.fit(X_tr, y_tr)
importances = forest.feature_importances_

indices = np.argsort(importances)[::-1]

print("Feature ranking:")

for f in range(len(features5) - 1):
    print("%d. feature: %s (%f)" % (f + 1, features5[indices[f]], importances[indices[f]]))

gradient = GradientBoostingClassifier(random_state=rs)

gradient.fit(X_tr, y_tr)
importances = gradient.feature_importances_

indices = np.argsort(importances)[::-1]

print("Feature ranking:")

for f in range(len(features5) - 1):
    print("%d. feature: %s (%f)" % (f + 1, features5[indices[f]], importances[indices[f]]))

"""
Explanation: Investigate feature importance
Feature importance for best models, Extra Trees, Random Forest, and Gradient Boosting.
End of explanation
"""


