<jupyter_start><jupyter_text>Phase 2: Conducting EDA and model construction using the master_reading.csv file that contains all necessary features and target variable<jupyter_code>import sys
print(sys.base_prefix)
#pip install --user scikit-learn
#!{sys.executable} -m pip install xgboost
#!{sys.executable} -m pip install lightgbm<jupyter_output><empty_output><jupyter_text>Importing all necessary libraries<jupyter_code>import pandas
pandas.__version__
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import pandas as pd
from sklearn.model_selection import cross_val_score, RepeatedStratifiedKFold, GridSearchCV, learning_curve, cross_val_predict
from sklearn.linear_model import LinearRegression, LassoCV, Ridge, LassoLarsCV,ElasticNetCV
from sklearn.ensemble import RandomForestRegressor, AdaBoostRegressor, ExtraTreesRegressor, GradientBoostingRegressor

from sklearn.pipeline import Pipeline
import warnings
import xgboost as xgb
import lightgbm as lgb
from scipy.stats import skew
from scipy import stats
from scipy.stats.stats import pearsonr
from scipy.stats import norm
from collections import Counter
from sklearn.preprocessing import StandardScaler, Normalizer, RobustScaler, OneHotEncoder
from sklearn.model_selection import train_test_split
warnings.filterwarnings('ignore')
%matplotlib inline
cd /Users/dansa/Documents/GitHub/Phase1/Data/MASTER<jupyter_output>C:\Users\dansa\Documents\GitHub\Phase1\Data\MASTER<jupyter_text>Loading the data and reformatting the school id column<jupyter_code>master_reading_new = pandas.read_csv("master_reading.csv")
master_reading_new['NCESSCH'] = master_reading_new['NCESSCH'].apply(lambda x: '{0:0>12}'.format(x))
master_reading_new.head()
master_reading_new.shape<jupyter_output><empty_output><jupyter_text>Inspecting the data file<jupyter_code>master_reading_new.columns<jupyter_output><empty_output><jupyter_text>Create a data frame with only the needed columns for further analysis<jupyter_code>reading=pd.DataFrame(master_reading_new, columns=[ 'NCESSCH','SCH_TYPE_x', 
       'TITLEI_STATUS','TEACHERS', 'FARMS_COUNT', 'Total_enroll_students',
       'SCH_FTETEACH_TOT', 'SCH_FTETEACH_CERT', 'SCH_FTETEACH_NOTCERT',
       'FTE_teachers_count', 'SalaryforTeachers', 'Total_SAT_ACT_students',
       'SCH_IBENR_IND_new', 'Total_IB_students', 'SCH_APENR_IND_new',
       'SCH_APCOURSES', 'SCH_APOTHENR_IND_new', 'Total_AP_other_students',
       'Total_students_tookAP', 'Income_Poverty_ratio','ALL_RLA00PCTPROF_1718_new'])
reading.head()
reading.rename(columns={'NCESSCH':'School_ID', 'SCH_TYPE_x':'School_type','FARMS_COUNT':'No.FARMS_students',
                       'SCH_FTETEACH_TOT':'FTE_teachcount','SCH_FTETEACH_CERT':'Certified_FTE_teachers','SCH_FTETEACH_NOTCERT':
                       'Noncertified_FTE_teachers','Total_SAT_ACT_students':'Students_participate_SAT_ACT','SCH_IBENR_IND_new':'IB_Indicator','SCH_APENR_IND_new':'AP_Indicator',
                        'SCH_APCOURSES':'No.ofAP_courses_offer','SCH_APOTHENR_IND_new':'Students_enroll_inotherAP?','ALL_RLA00PCTPROF_1718_new':'Percent_Reading_Proficient'}, inplace=True)
reading.describe()
print(reading.info())<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 14741 entries, 0 to 14740
Data columns (total 21 columns):
School_ID                       14741 non-null object
School_type                     14741 non-null int64
TITLEI_STATUS                   14741 non-null object
TEACHERS                        14741 non-null float64
No.FARMS_students               14741 non-null float64
Total_enroll_students           14741 non-null float64
FTE_teachcount                  14741 non-null float64
Certified_FTE_teachers          14741 non-null float64
Noncertified_FTE_teachers       14741 non-null float64
FTE_teachers_count              14741 non-null float64
SalaryforTeachers               14741 non-null float64
Students_participate_SAT_ACT    14741 non-null float64
IB_Indicator                    14741 non-null float64
Total_IB_students               14741 non-null float64
AP_Indicator                    14741 non-null float64
No.ofAP_courses_offer           14741 non-null float64
Students_enroll[...]<jupyter_text>We have 14,741 entries and no null values in any column. There are 21 columns, but we can drop the school_id and we'll want to split off the Percent_Reading_Proficient.The object type features should be strings.Let's take a quick look at some of the data.<jupyter_code>reading.hist(bins=50, figsize=(20,15))
plt.show()<jupyter_output><empty_output><jupyter_text>We can see that some features have most of their instances at or near zero and relatively few instances at higher values, in some cases much higher. Other features cluster close to zero and have long tails. We also see the percent_reading_proficient is almost normally distributed.<jupyter_code>#reading['TITLEI_STATUS'].value_counts()
sns.set_style('darkgrid')
_plt = sns.countplot(x='TITLEI_STATUS', data=reading)
_plt.set_title('School Title I status')
_plt.set_xticklabels(['Title I  schoolwide school','Not a Title I school','Title I schoolwide eligible- Title I targeted assistance program','Title I schoolwide eligible school-No program','Title I targeted assistance eligible school– No program','Title I targeted assistance school'])
_plt.set_xticklabels(_plt.get_xticklabels(), rotation=45, horizontalalignment='right')
plt.savefig('/Users/dansa/Documents/Title1_dist.png', dpi=300, bbox_inches='tight')
plt.show()
# Plot Histogram
sns.distplot(reading['Percent_Reading_Proficient'] , fit=norm);

# Get the fitted parameters used by the function
(mu, sigma) = norm.fit(reading['Percent_Reading_Proficient'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
            loc='best')
plt.ylabel('Frequency')
plt.title('Proficiency distribution')

fig = plt.figure()
res = stats.probplot(reading['Percent_Reading_Proficient'], plot=plt)
plt.show()

print("Skewness: %f" % reading['Percent_Reading_Proficient'].skew())
print("Kurtosis: %f" % reading['Percent_Reading_Proficient'].kurt())
sns.set_style('darkgrid')
Type_plt = sns.countplot(x='School_type', data=reading)
Type_plt.set_title('School_types')
Type_plt.set_xticklabels(Type_plt.get_xticklabels(), rotation=45, horizontalalignment='right')
Type_plt.set_xticklabels(["1-Regular School", "2-Special Education School", "3-Career and Technical School", "4-Alternative Education School"])
plt.savefig('/Users/dansa/Documents/Type_dist.png', dpi=300, bbox_inches='tight')
plt.show()
reading['School_type'].value_counts()<jupyter_output><empty_output><jupyter_text>There are very few alternative, charter and magnet high schools in the datasetGiven that the number of alt, charter, magnet and sped schools are so small in the dataset we might want to exclude that columns in our model tuning later<jupyter_code>reading['Pct_certified_teachers']=(reading['Certified_FTE_teachers']/reading['FTE_teachcount']*100) # lets find the percent of certified teachers
reading['Pct_noncertified_teachers']=(reading['Noncertified_FTE_teachers']/reading['FTE_teachcount']*100) # lets find the percent of noncertified teachers
reading['Salary_perFTE_teacher'] = reading['SalaryforTeachers']/reading['FTE_teachers_count'] # Lets find the salary per FTE in each school
reading['IPR_estimate'] = reading['Income_Poverty_ratio'] #Income poverty ratio is reported as a percent 
#reading['Percent_Reading_Proficient'] = reading['Percent_Reading_Proficient']/100
reading_clean=reading.drop(['School_ID','Certified_FTE_teachers', 'Noncertified_FTE_teachers','FTE_teachcount','FTE_teachers_count','SalaryforTeachers','Income_Poverty_ratio' ], axis=1)
reading_clean.info()
reading_clean.dtypes
reading_clean.describe()
sns.heatmap(reading_clean.isnull(),yticklabels=False,cbar=True,cmap='viridis')
reading_clean.shape<jupyter_output><empty_output><jupyter_text>3. Looking for Correlations and VisualizingWe should calculate data correlations and plot a scatter matrix.For training the ML models, we'll want to separate the Percent_Reading_Proficient from the rest of the data. But for investigating correlations, we'll want to include the target.<jupyter_code>reading_clean=reading_clean[['School_type', 'TITLEI_STATUS', 'TEACHERS', 'No.FARMS_students',
       'Total_enroll_students', 'Students_participate_SAT_ACT', 'IB_Indicator',
       'Total_IB_students', 'AP_Indicator', 'No.ofAP_courses_offer',
       'Students_enroll_inotherAP?', 'Total_AP_other_students',
       'Total_students_tookAP',
       'Pct_certified_teachers', 'Pct_noncertified_teachers',
       'Salary_perFTE_teacher', 'IPR_estimate','Percent_Reading_Proficient']]
correlation_matrix = reading_clean.corr()
correlation_matrix['Percent_Reading_Proficient'].sort_values(ascending=False)<jupyter_output><empty_output><jupyter_text>It seems like a few features (IPR_estimate, No.ofAP_courses_offer, Total_students_tookAP) have a weak to moderate positive correlation to the target (Percent_Reading_Proficient), and a couple are somewhat negatively correlated (School_type).IPR_estimate is the Neighborhood Income Poverty Ratio.No.ofAP_courses_offer is the count of AP courses offered at the school.Total_students_tookAP is the number of students who took an AP course.School_type is refers to whether the school is a "1-Regular School, 2-Special Education School, 3-Career and Technical School and 4-Alternative Education School"We can look at a heatmap of the correlations of all numeric features to visualize which features are correlated.<jupyter_code># correlation matrix heatmap
plt.figure(figsize=(28,15))
corr_heatmap = sns.heatmap(correlation_matrix, annot=True, linewidths=0.2, center=0, cmap="RdYlGn")
corr_heatmap.set_title('Correlation Heatmap')
plt.savefig('/Users/dansa/Documents/corr_heatmap.png', dpi=300, bbox_inches='tight')
#test
corr_pairs = {}
feats = correlation_matrix.columns
for x in feats:
    for y in feats:
        if x != y and np.abs(correlation_matrix[x][y]) >= 0.7:  # which pairs are strongely correlated?
            if (y, x) not in corr_pairs.keys():
                corr_pairs[(x, y)] = correlation_matrix[x][y]
corr_pairs
weaker_label = []
for pair in corr_pairs:
    if np.abs(correlation_matrix[pair[0]]['Percent_Reading_Proficient']) < np.abs(correlation_matrix[pair[1]]['Percent_Reading_Proficient']):
        weaker_label.append(pair[0])
    else:
        weaker_label.append(pair[1])
poss_redundant_feats = set(weaker_label)
poss_redundant_feats
cov_matrix = reading_clean.cov()
plt.figure(figsize=(28,15))
covar_heatmap = sns.heatmap(data=cov_matrix, cmap='coolwarm', linewidth=0.2) 
covar_heatmap.set_title('Covariance Heatmap')
plt.savefig('/Users/dansa/Documents/covar_heatmap.png', dpi=300, bbox_inches='tight')
attrs = ['IPR_estimate','No.ofAP_courses_offer','Total_students_tookAP','Percent_Reading_Proficient']
sns.set(style='ticks', color_codes=True)
_ = sns.pairplot(data=correlation_matrix[attrs], height=3, aspect=1, kind='scatter', plot_kws={'alpha':0.9})
#from pandas.plotting import scatter_matrix
#attributes = ["Percent_Reading_Proficient", "IPR_estimate", "No.ofAP_courses_offer","Total_students_tookAP","Total_AP_other_students"]
#scatter_matrix(reading_clean[attributes], figsize=(12, 8))
sns.jointplot(x="IPR_estimate", y="Percent_Reading_Proficient", data=reading_clean)<jupyter_output><empty_output><jupyter_text>ML prep Separate labelsLet's separate out the target from the predicting features.<jupyter_code>from sklearn.model_selection import train_test_split
train, test = train_test_split(reading_clean, test_size=0.2, random_state=0)
df_X_train = train.drop(['Percent_Reading_Proficient'], axis=1)
df_y_train = train['Percent_Reading_Proficient'].copy()

df_X_test = test.drop(['Percent_Reading_Proficient'], axis=1)
df_y_test = test['Percent_Reading_Proficient'].copy()
df_X_train.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
Int64Index: 11792 entries, 3420 to 2732
Data columns (total 17 columns):
School_type                     11792 non-null int64
TITLEI_STATUS                   11792 non-null object
TEACHERS                        11792 non-null float64
No.FARMS_students               11792 non-null float64
Total_enroll_students           11792 non-null float64
Students_participate_SAT_ACT    11792 non-null float64
IB_Indicator                    11792 non-null float64
Total_IB_students               11792 non-null float64
AP_Indicator                    11792 non-null float64
No.ofAP_courses_offer           11792 non-null float64
Students_enroll_inotherAP?      11792 non-null float64
Total_AP_other_students         11792 non-null float64
Total_students_tookAP           11792 non-null float64
Pct_certified_teachers          11792 non-null float64
Pct_noncertified_teachers       11792 non-null float64
Salary_perFTE_teacher           11792 non-null float64
IPR_estimate[...]<jupyter_text>Transform Categorical FeaturesSince these categorical features don't appear to have an inherent ordering, let's try encoding them as one-hot vectors for better ML performance.<jupyter_code>onehot_encoder = OneHotEncoder()
#df_X_test.head()
onehot_encoder.fit(df_X_train)
train_data_onehot = onehot_encoder.transform(df_X_train)
test_data_onehot = onehot_encoder.transform(df_X_test)
train_data_onehot = pd.get_dummies(df_X_train, columns=['TITLEI_STATUS'], prefix=['TITLEI_STATUS'])
train_data_onehot.head()

test_data_onehot = pd.get_dummies(df_X_test, columns=['TITLEI_STATUS'], prefix=['TITLEI_STATUS'])

train_data_onehot.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
Int64Index: 11792 entries, 3420 to 2732
Data columns (total 22 columns):
School_type                     11792 non-null int64
TEACHERS                        11792 non-null float64
No.FARMS_students               11792 non-null float64
Total_enroll_students           11792 non-null float64
Students_participate_SAT_ACT    11792 non-null float64
IB_Indicator                    11792 non-null float64
Total_IB_students               11792 non-null float64
AP_Indicator                    11792 non-null float64
No.ofAP_courses_offer           11792 non-null float64
Students_enroll_inotherAP?      11792 non-null float64
Total_AP_other_students         11792 non-null float64
Total_students_tookAP           11792 non-null float64
Pct_certified_teachers          11792 non-null float64
Pct_noncertified_teachers       11792 non-null float64
Salary_perFTE_teacher           11792 non-null float64
IPR_estimate                    11792 non-null float64
TITLEI_STAT[...]<jupyter_text>Scale FeaturesWe can check out the statistics for each feature, do they need to be normalized?<jupyter_code>train_data_onehot.describe()<jupyter_output><empty_output><jupyter_text>We'll probably want to scale these features using normalization or standardization.<jupyter_code>sc= StandardScaler()
train_scaler = sc.fit(train_data_onehot)
#test_scaler = sc.fit(test_data_onehot)
train_data_standardized = train_scaler.transform(train_data_onehot)
test_data_standardized = train_scaler.transform(test_data_onehot)
print(train_data_standardized)
print(train_data_standardized.mean(axis=0))
print(train_data_standardized.std(axis=0))
train_data_std = pd.DataFrame(train_data_standardized, columns=train_data_onehot.columns)
test_data_std = pd.DataFrame(test_data_standardized, columns=test_data_onehot.columns)
#train_data_std.describe()<jupyter_output><empty_output><jupyter_text>That should work better, the standard deviation for each feature is 1 and the mean is ~0.  Using a multiple linear regression model<jupyter_code>from sklearn.linear_model import LinearRegression
regressor = LinearRegression()
regressor.fit(train_data_std, df_y_train)
y_pred = regressor.predict(test_data_std)
from termcolor import colored as cl
print(cl('R-Squared :', attrs = ['bold']), regressor.score(test_data_std, df_y_test))
# Method 2 Using Statsmodels
import statsmodels.api as sm
from scipy import stats
y = list(df_y_train)
X2 = sm.add_constant(train_data_std) 
# Define the model
est = sm.OLS(y, X2)
# Fit the model
est2 = est.fit()
print("summary()\n",est2.summary())
print("pvalues\n",est2.pvalues)
print("tvalues\n",est2.tvalues)
print("rsquared\n",est2.rsquared)
print("rsquared_adj\n",est2.rsquared_adj)
sns.distplot(y_pred, hist = False, color = 'r', label = 'Predicted Values')
sns.distplot(df_y_test, hist = False, color = 'b', label = 'Actual Values')
plt.title('Actual vs Predicted Values', fontsize = 16)
plt.xlabel('Values', fontsize = 12)
plt.ylabel('Frequency', fontsize = 12)
plt.legend(loc = 'upper left', fontsize = 13)
plt.savefig('ap.png')
df = pandas.DataFrame({'Actual':  df_y_test, 'Predicted': y_pred})
df
from sklearn import metrics
print('Mean Absolute Error:', metrics.mean_absolute_error( df_y_test, y_pred))
print('Mean Squared Error:', metrics.mean_squared_error( df_y_test, y_pred))
print('Root Mean Squared Error:', np.sqrt(metrics.mean_squared_error( df_y_test, y_pred)))
import matplotlib.pyplot as plt
from sklearn import linear_model
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
regr = linear_model.LinearRegression()
regr.fit(train_data_std, df_y_train)
y_pred = np.round(regr.predict(test_data_std))
meansqr=[]
meanabs=[]
r2=[]
rmse=[]
meansqr.append(mean_squared_error(df_y_test, y_pred))
meanabs.append(abs(df_y_test-y_pred).mean())
rmse.append(np.sqrt(mean_squared_error(df_y_test, y_pred)))
r2.append(r2_score(df_y_test, y_pred))
print("Mean squared error: %.2f"% mean_squared_error(df_y_test, y_pred))
print("Mean absolute error: %.2f"% mean_absolute_error(df_y_test, y_pred))
print('Root Mean Squared Error:%.2f'% np.sqrt(mean_squared_error( df_y_test, y_pred)))
print("r2 score: %.2f"% r2_score(df_y_test, y_pred))
#plotting y_pred and y_test
t = np.arange(0,len(y_pred) , 1)
plt.plot(t,y_pred,t,df_y_test)
# RandomForestRegressor
from sklearn.ensemble import RandomForestRegressor
# Instantiate model with 1000 decision trees
rf = RandomForestRegressor(n_estimators = 1000, random_state = 42)
# Train the model on training data
rf.fit(train_data_std, df_y_train);
y_pred =np.round(rf.predict(test_data_std))
meansqr.append(mean_squared_error(df_y_test, y_pred))
meanabs.append(abs(df_y_test-y_pred).mean())
rmse.append(np.sqrt(mean_squared_error(df_y_test, y_pred)))
r2.append(r2_score(df_y_test, y_pred))
print("Mean squared error: %.2f"% mean_squared_error(df_y_test, y_pred))
print("Mean absolute error: %.2f"% mean_absolute_error(df_y_test, y_pred))
print('Root Mean Squared Error:%.2f'% np.sqrt(mean_squared_error( df_y_test, y_pred)))
print("r2 score: %.2f"% r2_score(df_y_test, y_pred))
#plotting y_pred and y_test
t = np.arange(0,len(y_pred) , 1)
plt.plot(t,y_pred,t,df_y_test)
#svm regressor
from sklearn.svm import SVR
regressor=SVR(kernel="linear",epsilon=1.0,degree=3)
regressor.fit(train_data_std, df_y_train)
y_pred=regressor.predict(test_data_std)
meansqr.append(mean_squared_error(df_y_test, y_pred))
meanabs.append(abs(df_y_test-y_pred).mean())
rmse.append(np.sqrt(mean_squared_error(df_y_test, y_pred)))
r2.append(r2_score(df_y_test, y_pred))
print("Mean squared error: %.2f"% mean_squared_error(df_y_test, y_pred))
print("Mean absolute error: %.2f"% mean_absolute_error(df_y_test, y_pred))
print('Root Mean Squared Error:%.2f'% np.sqrt(mean_squared_error( df_y_test, y_pred)))
print("r2 score: %.2f"% r2_score(df_y_test, y_pred))
#plotting y_pred and y_test
t = np.arange(0,len(y_pred) , 1)
plt.plot(t,y_pred,t,df_y_test)
#knearest neighbhors
from sklearn import neighbors
n_neighbors=5
knn=neighbors.KNeighborsRegressor(n_neighbors,weights='uniform')
knn.fit(train_data_std, df_y_train)
y_pred=knn.predict(test_data_std)
meansqr.append(mean_squared_error(df_y_test, y_pred))
meanabs.append(abs(df_y_test-y_pred).mean())
rmse.append(np.sqrt(mean_squared_error(df_y_test, y_pred)))
r2.append(r2_score(df_y_test, y_pred))
print("Mean squared error: %.2f"% mean_squared_error(df_y_test, y_pred))
print("Mean absolute error: %.2f"% mean_absolute_error(df_y_test, y_pred))
print('Root Mean Squared Error:%.2f'% np.sqrt(mean_squared_error( df_y_test, y_pred)))
print("r2 score: %.2f"% r2_score(df_y_test, y_pred))
#plotting y_pred and y_test
t = np.arange(0,len(y_pred) , 1)
plt.plot(t,y_pred,t,df_y_test)
objects=('LinearReg','RandForReg','SVMregg','Knn')
plt.bar(np.arange(len(meansqr)),meansqr)
plt.xticks(np.arange(len(meansqr)), objects)
plt.title('Mean Square Error')
plt.show()
plt.bar(np.arange(len(meanabs)),meanabs)
plt.xticks(np.arange(len(meanabs)), objects)
plt.title('Mean Absolute Error')
plt.show()
plt.bar(np.arange(len(rmse)),rmse)
plt.xticks(np.arange(len(rmse)), objects)
plt.title('Root mean squared error')
plt.show()
plt.bar(np.arange(len(r2)),r2)
plt.xticks(np.arange(len(r2)), objects)
plt.title('r2_score')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10727969348659004, 'fraction_numerical': 0.026867816091954023, 'mean_word_length': 6.334387074113102, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 62, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '350067', 'license': 'CC0-1.0', 'n_tokens_mistral': 7857, 'n_tokens_neox': 7104, 'n_words': 1466, 'path': 'Phase2/Notebooks/.ipynb_checkpoints/Reading_data_EDA_Modeling_4.11-checkpoint.ipynb', 'repo_name': 'dansari12/Predicting-high-school-assessment-proficiency', 'size': '20880'}