<jupyter_start><jupyter_text>Winter Resort Visitation Prediction - Resort 1*Capstone Data Science Project / Flatiron Data Science Bootcamp*By: [<NAME>](https://github.com/lperconti) OverviewIn this project I look at three winter resorts and utilize weather and event data to predict visitation for each resort. This notebook looks at Resort 1.  The Stakeholder / Business ProblemFor a full, detailed stakeholder and business problem introduction, see the [Read.Me](https://github.com/lperconti/Visitation/blob/main/README.md) file.   Exploratory Data Analysis Resort 1 Background: **Resort 1** is close to a large metropolitan area (approx 455,000 people). The resort's lift capacity is 9,980 riders/hour. The resort generously gave four seasons' worth of daily visit totals in an excel document for use in this project. Data was given on the condition that the specific resort not be identified publically. Therefore the name and any location-based information associated with this dataset are withheld from this repository. The daily visits in the given data were further broken down into season pass holder visits, day ticket visits, day visits, and night vists. A couple notes on getting the data for this particular ski area: * `is_school_out` indicates if the public school district did not have class on a day that normally would have school (weekday). These included Winter Break, MLK Day, Presidents' Day, and Spring Break. Weekends (Sundays and Mondays) surrounding school closures such as Thanksgiving, Presidents Day, were not counted as a `SCHOOL_OUT` day. I was able to source these dates from the local school district's website. There could be other schools in the surrounding area such as colleges, private schools with school out days that this data does not capture. * `is_holiday` indicates if the date is a federally (USA) recognized holiday. Examples include: Thanksgiving, Christams Eve, Christmas Day, New Years Eve, New Years Day, MLK Day, and Presidents' Day. * Mountain Weather Data was gathered from a NOAA weather station on the mountain where the resort is located.* Town weather Data was gathered from the metropolitan area's airport.  Data Type Keys - WeatherAll weather measurements are in Standard (degrees fahrenheit, inches)- TMAX - Maximum temperature- WESD - Water equivalent of snow on the ground- TAVG - Average Temperature.- TMIN - Minimum temperature- PRCP - Precipitation- TOBS - Temperature at the time of observation- SNWD - Snow depth- AWND - Average wind speed- PRCP - Total Daily Precipitation- SNOW - Total Daily Snow- WSF2 - Fastest 2-minute wind speed (mph)- WSF5 - Fastest 5-minute wind speed (mph)- T - 'Town'- M - 'Mountain' After cleaning Visitation stats: **Daily Visits:**count     484.000000mean     2421.783058std      1709.932584min        48.0000025%      1037.50000050%      2015.50000075%      3503.500000max      7144.000000<jupyter_code>#imports

#libraries for data manipulation and visualization
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
%matplotlib inline

#libraries for modeling
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import scale
from sklearn.preprocessing import OneHotEncoder
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.compose import ColumnTransformer
from sklearn.tree import DecisionTreeRegressor
from sklearn.model_selection import GridSearchCV
from sklearn.ensemble import RandomForestRegressor
from sklearn.ensemble import GradientBoostingRegressor
import xgboost as xg

#libraries for evaluation
from sklearn.inspection import permutation_importance
from sklearn.model_selection import cross_validate
from sklearn.model_selection import cross_val_score
#initial look at the dataset
df_1 = pd.read_csv('Data/Resort_1.csv', thousands=',')
df_1.info()
#See if there's any duplicates: 
df_1.duplicated().value_counts()<jupyter_output><empty_output><jupyter_text>To get day of week and month categories, we will need to make the `Date` column a datetime object. Day of week and month categories will be used in our models to see if they have any predictive power.<jupyter_code>#Date dealings - make the 'Date' column a datetime object, create day of week and month category to do some EDA and use as variables
#make the 'Date' Column a datetime type
df_1['Date'] = pd.to_datetime(df_1['Date'])

#The day_of_week_num column created with Monday=0, Sunday=6. Weekends = 5 & 6, Friday = 4
df_1['day_of_week_num'] = df_1['Date'].dt.dayofweek

#creating Month Column: 1 = January, 2 = February, ect.
df_1['month'] = df_1['Date'].dt.month

#no need for day of week column in a string type anymore
df_1.drop(['Day of week'], axis=1, inplace = True)

df_1['month'].value_counts()
#visiualize date category results: 
df_1.describe()
#renaming `Total VIsits Alpine` into something more succinct:
df_1.rename(columns = {'Total VIsits Alpine':'TOTAL_VISITS'}, inplace = True)<jupyter_output><empty_output><jupyter_text>Next we want to look at how many Null Values are in our dataset. Null Values may result in errors when modeling, and we will need to decide how to deal with missing data:<jupyter_code>df_1.isnull().sum()<jupyter_output><empty_output><jupyter_text>It appears the missing values are in `Night Ticket`, `SP Day`, and `SP Night`, which are all visitation targets. We will likely drop these columns and focus on only prediciting `TOTAL_VISITS` for now.   Correlation Matrix Let's start exploring correlation, taking out `TOTAL_VISITS` (our target), SP, Day and Night related ticket totals to see if any of our variables are highly correlated. Also taking out `Date` as the specific date is not a predictor.<jupyter_code>#removing predictor (visit) and date columns for correlation testing
df_1_corr = df_1.drop(columns=['Date', 'TOTAL_VISITS', 'SP Day', 'SP Night', 'Night Ticket'])

# Create a df with the target as the first column,
# then compute the correlation matrix
heatmap_data = df_1_corr
corr = heatmap_data.corr()

# Set up figure and axes
fig, ax = plt.subplots(figsize=(20, 15))

# Plot a heatmap of the correlation matrix, with both
# numbers and colors indicating the correlations
sns.heatmap(data=corr, mask=np.triu(np.ones_like(corr, dtype=bool)), ax=ax, annot=True, 
            cbar_kws={"label": "Correlation", "orientation": "horizontal", "pad": .2, "extend": "both"})

ax.set_title("Heatmap of Correlation Between Attributes (Including Target)");
#identifying the highest correlated variables, correlation >.75

df = df_1_corr.corr().abs().stack().reset_index().sort_values(0, ascending=False)
df['pairs'] = list(zip(df.level_0, df.level_1))
df.set_index(['pairs'], inplace = True)
df.drop(columns=['level_1', 'level_0'], inplace = True)

# cc for correlation coefficient
df.columns = ['cc']
df.drop_duplicates(inplace=True)
df[(df.cc>.75) & (df.cc<1)]<jupyter_output><empty_output><jupyter_text>Multicollinairity decisionBased on the information above it makes sense to drop some data, especially the weather data: **Weather:** - `T_WSF5` (fastest 5 min wind speed in the metropolotan(town) area)- `M_WESD` (water equivalent snow on the ground on the mountain)- `M_TMIN` + `T_TMIN` (temperature min on the mountain and in town)- `M_TMAX` + `T_TMAX` (temperature max on the mountain and in town)- `T_TOBS` + `M_TOBS` (temperature at time of observation on mtn and in town)- `T_WSF2` (fastest 2 min wind speend, keeping the average wind speed for town)**Ticketing:**To keep things simple I think for now I am going to remove the ticketing and visit informatoin that breaks down season pass vs ticket and night vs day ticketing. I may want to use these as a prediction if there's time or the business need becomes obvious. Unfortunately I don't have pre-sale data to work with. - `SP Night`- `SP Day`- `Day Tickets`- `Night Ticket`- `Total Tickets`- `Total SP`<jupyter_code>df_1_preprocessed = df_1.drop(columns=['T - WSF5', 'M_WESD', 'M_TMIN', 'T - TMIN', 'M_TMAX', 'T - TMAX', 'M_TOBS', 'T - WSF2', 'SP Night', 'SP Day', 'Day Tickets', 'Night Ticket', 'Total Ticket', 'Total SP'])
df_1_preprocessed.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 567 entries, 0 to 566
Data columns (total 14 columns):
 #   Column           Non-Null Count  Dtype         
---  ------           --------------  -----         
 0   Date             567 non-null    datetime64[ns]
 1   is_school_out    567 non-null    int64         
 2   is_holiday       567 non-null    int64         
 3   TOTAL_VISITS     567 non-null    float64       
 4   T - AWND         567 non-null    float64       
 5   T - PRCP         567 non-null    float64       
 6   T - SNOW         567 non-null    float64       
 7   T - SNWD         567 non-null    float64       
 8   T - TAVG         567 non-null    int64         
 9   M_PRCP           567 non-null    float64       
 10  M_SNWD           567 non-null    int64         
 11  M_TAVG           567 non-null    int64         
 12  day_of_week_num  567 non-null    int64         
 13  month            567 non-null    int64         
dtypes: datetime64[ns](1), float64(6), int64(7)[...]<jupyter_text>0 visit daysTaking a look at rows with 0 visits (these are the days the resort is closed)<jupyter_code>df_1_preprocessed[df_1_preprocessed.TOTAL_VISITS == 0]
#dropping the 0 visit days: 
df_1_preprocessed = df_1_preprocessed[df_1_preprocessed.TOTAL_VISITS != 0]
df_1_preprocessed['TOTAL_VISITS'].value_counts()
#seeing what's left: 
df_1_preprocessed['TOTAL_VISITS'].describe()<jupyter_output><empty_output><jupyter_text>Indexing with Date<jupyter_code>#df_1_preprocessed.reset_index(inplace=True)
df_1_preprocessed.set_index('Date', inplace=True)
df_1_preprocessed.head(2)<jupyter_output><empty_output><jupyter_text>Visualizations<jupyter_code>#visualizing the visitations over time
df_1_preprocessed['TOTAL_VISITS'].plot(figsize = (40, 15));<jupyter_output><empty_output><jupyter_text>Visualizations: BoxplotThis basic Boxplots, or sometimes called a box-and-whisker plot, is helpful to visualize distributions of our visitation data. The box shows the quartiles of the data and whiskers extend to show the rest of the distribution. It's helpful to see how the data is distributed.<jupyter_code>#visualizing daily visits by month: 
medians = df_1.groupby(['month'])['TOTAL_VISITS'].median()
print(medians)

import matplotlib as mpl

fig, ax = plt.subplots(figsize=(15, 12))
sns.boxplot(x='month', y='TOTAL_VISITS', data=df_1_preprocessed, width=.7)

medians = df_1_preprocessed.groupby(['month'])['TOTAL_VISITS'].median()
print(medians)

ax.set_xlabel("")

ax.set_ylabel("Total Daily Visits", fontsize=24)

ax.set_title("Resort One Daily Visits by Month")

plt.yticks(fontsize=20)

plt.xticks([0, 1, 2, 3, 4, 5], ['Jan', 'Feb', 'Mar', 'Apr', 'Nov', 'Dec'], fontsize=20)

ax.yaxis.set_major_formatter(mpl.ticker.StrMethodFormatter('{x:,.0f}'))

plt.savefig('Visuals/resort1visits.png', dpi=300);<jupyter_output>month
1     2482.5
2     2602.0
3     1221.5
4      247.5
11       0.0
12    1687.0
Name: TOTAL_VISITS, dtype: float64
month
1     2571.0
2     2602.0
3     1337.0
4      449.0
11     193.0
12    2525.0
Name: TOTAL_VISITS, dtype: float64<jupyter_text>Let's look at the relationship between visitors and Precipitation on the Mountain:<jupyter_code>fig, ax = plt.subplots()

ax.scatter(df_1_preprocessed['M_PRCP'], df_1_preprocessed["TOTAL_VISITS"], alpha=0.5)
#visualization of target - Visits
fig, ax = plt.subplots(figsize=(10, 5))
ax.hist(df_1_preprocessed['TOTAL_VISITS'], bins=100);
#Relationship between mountain precipitation and visitation
plt.scatter(y=df_1_preprocessed['M_PRCP'], x=df_1_preprocessed['TOTAL_VISITS'], alpha=.4)
plt.savefig('Visuals/resort1precip.png', dpi=300);
#relationship between snowdepth and visitation
plt.scatter(y=df_1_preprocessed['M_SNWD'], x=df_1_preprocessed['TOTAL_VISITS'], alpha=.4)
plt.savefig('Visuals/resort1msnowdepth.png', dpi=300);<jupyter_output><empty_output><jupyter_text>At first glance, calendar events seem to have more impact on visitation:<jupyter_code>Holiday_averages=df_1_preprocessed.groupby(df_1_preprocessed['is_holiday'])['TOTAL_VISITS'].mean()
print(Holiday_averages)

Schoolout_averages=df_1_preprocessed.groupby(df_1_preprocessed['is_school_out'])['TOTAL_VISITS'].mean()
print(Schoolout_averages)

fig, ax = plt.subplots()
sns.stripplot(x='M_PRCP', y='TOTAL_VISITS', data=df_1_preprocessed)
fig, ax = plt.subplots()
sns.swarmplot(x='month', y='TOTAL_VISITS', data=df_1_preprocessed, size=3, palette='twilight', edgecolor='none')

ax.set_xlabel("")
ax.set_title("Resort 1 Daily vists by month")
ax.set_ylabel("Daily Visits", fontsize=15)
plt.yticks(fontsize=14)
plt.xticks([0, 1, 2, 3, 4, 5], ['Jan', 'Feb', 'Mar', 'Apr', 'Nov', 'Dec'], fontsize=14)

ax.yaxis.set_major_formatter(mpl.ticker.StrMethodFormatter('{x:,.0f}'))
plt.savefig('Visuals/resort1swarmplotmonth.png', dpi=300);<jupyter_output><empty_output><jupyter_text>Saving Cleaned Resort 3 Data to CSV for reference:<jupyter_code>df_1_preprocessed.to_csv("Data/DF_1_clean.csv")<jupyter_output><empty_output><jupyter_text>First Simple Model**Linear Regression**We will be using a regression analysis to estimate the relationships between the dependent (Visits) and independent variables in our dataset. A simple Linear Regression model makes sense to build ourbaseline model with. The R-squared score tells us how much of the total variability in the data is explained by the model and is very useful when comparing different models. This will help us establish and determine if there is a realationship between our independent variables and the visitation this ski resort experiences. We can also use regression analysis to determine which specific features are good predictors for the target and which are not. Regression analysis is fairly interpretable, and being able to explain these models to a non-technical audience at ski resorts wil be key to the models being used for prediction. The final reason for using Regression Analysis is the ability to generate predictions, which is what this analysis is all about: seeing if we can predict visitation for ski resorts. To set up modeling, we will first split our data into X (independent) and Y(dependent) dataframes. Then we will split the data into train and test splits. Since we are working with a fairly small dataset, we will be using cross validation on our training data to score the models without having to utlize the test data until a best model has been determined.<jupyter_code>#Split into test and train X and Y
X = df_1_preprocessed.drop(columns=['TOTAL_VISITS'])
y = df_1_preprocessed['TOTAL_VISITS']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
#visualizing our X to ensure the dataframe split correctly
X.head()
# creating categorical lists and Transformer: 
cat = ['is_school_out', 'is_holiday', 'day_of_week_num', 'month']
non_cat = ['T - AWND', 'T - PRCP', 'T - SNOW', 'T - SNWD', 'T - TAVG', 'M_PRCP', 'M_SNWD', 'M_TAVG']

Transformer = ColumnTransformer([
	('ohe', OneHotEncoder(), cat),
	('scaler', StandardScaler(), non_cat),
	], remainder ='passthrough')<jupyter_output><empty_output><jupyter_text>Baseline Linear Regression model with One Hot Encoding and Standard Scaler in a pipeline<jupyter_code>#setting up the pipe and training the model: 

LRPipe = Pipeline([('trans', Transformer),
                ('linreg', LinearRegression())])

LRPipe.fit(X_train, y_train)

#cross validation and scoring
cv_LR = cross_validate(LRPipe, X_train, y_train, return_train_score=True, cv=5)
cv_score_LR = cross_val_score(LRPipe, X_train, y_train, cv = 5)

print(f'Cross Validation Scores, Linear Regression: {cv_LR}')

print(f'mean CV Score, Linear Regression, Resort 3: {cv_score_LR.mean()}')<jupyter_output>Cross Validation Scores, Linear Regression: {'fit_time': array([0.04224586, 0.03379774, 0.04090691, 0.01570988, 0.01338315]), 'score_time': array([0.0224452 , 0.032547  , 0.008744  , 0.00675607, 0.00983882]), 'test_score': array([0.47984437, 0.53637444, 0.73232436, 0.66131302, 0.48757651]), 'train_score': array([0.68906124, 0.67728235, 0.63319602, 0.64723468, 0.68273521])}
mean CV Score, Linear Regression, Resort 3: 0.5794865402447538<jupyter_text>With Encoding and scaling, the Basic Linear Regression CV score is 0.579There is definitly some room for improvement, but we're doing OK so far.  Tree Regressorwhy tree next?<jupyter_code>Transformer = ColumnTransformer([
	('ohe', OneHotEncoder(), cat),
	# ('scaler', StandardScaler(), non_cat), #no need to scale on Trees
	], remainder ='passthrough')

#Pipeline second Step: treeregressor
dt_clf2 = DecisionTreeRegressor(random_state = 42)

#pipeline
pipeline_dt = Pipeline([('transformer', Transformer),
						('dtr', dt_clf2)])

#set up grid search paramater grid
dt_param_grid = {
    'dtr__min_samples_split': [2, 4, 6, 8],
    'dtr__max_depth': [None, 2, 3, 4, 5, 6],
	'dtr__splitter': ['random', 'best']
}

#set up gridsearch:
gridsearch_dt = GridSearchCV(estimator=pipeline_dt,
						  param_grid=dt_param_grid,
						  cv=5)

#fit gridsearch
gridsearch_dt.fit(X_train, y_train)

print("Best Parameter Combination for Tree Regressor Found During Grid Search:")
print(gridsearch_dt.best_params_)
Transformer = ColumnTransformer([
	('ohe', OneHotEncoder(), cat),
	# ('scaler', StandardScaler(), non_cat), #no need to scale on Trees
	], remainder ='passthrough')

#Pipeline second Step: treeregressor
dt_clf3 = DecisionTreeRegressor(max_depth=5, min_samples_split=4, splitter='best', random_state = 42)

#pipeline
pipeline_dt_best = Pipeline([('transformer', Transformer),
						('dtr', dt_clf3)])


pipeline_dt_best.fit(X_train, y_train)

cv = cross_validate(pipeline_dt_best, X_train, y_train, return_train_score=True, cv=5)
cv_score = cross_val_score(pipeline_dt_best, X_train, y_train, cv = 5)
print(f'Cross Validation Scores bDecision treeest Decision Tree Regressor: {cv}')

print(f'mean CV Score: {cv_score.mean()}')<jupyter_output>Cross Validation Scores bDecision treeest Decision Tree Regressor: {'fit_time': array([0.03358006, 0.02057004, 0.02237988, 0.01736331, 0.02448511]), 'score_time': array([0.00809598, 0.00996399, 0.01013017, 0.02309561, 0.01484895]), 'test_score': array([0.44124527, 0.54965897, 0.45015904, 0.688071  , 0.41617338]), 'train_score': array([0.73431034, 0.76736684, 0.7726827 , 0.7539595 , 0.7848853 ])}
mean CV Score: 0.509061532768853<jupyter_text>Decision Tree model with optomized hyperparamaters performed worse than Linear Regression (LR Mean CV score = .579) at a mean CV score of 0.509. Trying Random Forests:Our next model will be Random Forests. While Random Forests tend to take more time to run, they tend to do well with 'weak learners', which we may have here. Random forests combine several decision trees which individually are weak performers and combine them together for a stronger learner.<jupyter_code>Transformer = ColumnTransformer([
	('ohe', OneHotEncoder(), cat),
	# ('scaler', StandardScaler(), non_cat), #no need to scale on Trees
	], remainder ='passthrough')

#Pipeline second Step: treeregressor
rfr = RandomForestRegressor(random_state = 42)

#pipeline
pipeline_rfr = Pipeline([('transformer', Transformer),
						('dtr', rfr)])

#set up grid search paramater grid
rfr_param_grid = {
    'dtr__n_estimators': [2, 4, 6, 8, 10, 12, 14, 16],
    'dtr__max_depth': [None, 2, 3, 4, 5, 6, 8, 10],
	'dtr__max_features': ['auto', 'sqrt', 'log2'],
	'dtr__criterion': ['mse', 'mae']
}

#set up gridsearch:
gridsearch_rfr = GridSearchCV(estimator=pipeline_rfr,
						  param_grid=rfr_param_grid,
						  cv=5)

#fit gridsearch
gridsearch_rfr.fit(X_train, y_train)

print("Best Parameter Combination Found During Grid Search:")
print(gridsearch_rfr.best_params_)

#modeling with best paramaters found in Grid Search: 
Transformer = ColumnTransformer([
	('ohe', OneHotEncoder(), cat),
	# ('scaler', StandardScaler(), non_cat), #no need to scale on Trees
	], remainder ='passthrough')

#Pipeline second Step: treeregressor
rfr_best = RandomForestRegressor(criterion = 'mse', max_depth = None, max_features = 'sqrt', n_estimators = 10, random_state = 42)

#pipeline
pipeline_rfr_best = Pipeline([('transformer', Transformer),
						('dtr', rfr_best)])


pipeline_rfr_best.fit(X_train, y_train)

cv = cross_validate(pipeline_rfr_best, X_train, y_train, return_train_score=True, cv=5)
cv_score = cross_val_score(pipeline_rfr_best, X_train, y_train, cv = 5)
print(f'Cross Validation Scores best Random Forest Regressor: {cv}')

print(f'mean CV Score: {cv_score.mean()}')<jupyter_output>Cross Validation Scores best Random Forest Regressor: {'fit_time': array([0.02399516, 0.05153203, 0.03333497, 0.02297592, 0.02309179]), 'score_time': array([0.00611687, 0.00928688, 0.00577426, 0.00757027, 0.0045619 ]), 'test_score': array([0.55986679, 0.67336618, 0.74058046, 0.6732444 , 0.54516542]), 'train_score': array([0.92655425, 0.91821211, 0.9111    , 0.90889369, 0.91310719])}
mean CV Score: 0.6384446490641958<jupyter_text>As we can see, the Random Forest Model performed better than the decision Tree. Random Forest regressor mean CV score = .6638, which is an improvement on both the Linear Regression and decision tree Model.   Gradient Boosting: Gradient boosting is another model that combines weak learners and individual trees to create a stronger learner. Gradient boosting builds one tree at a time, and each new tree helps correct errors made by the previous tree. This is also a computationally expensive model and takes a bit of time to run, and we will have to watch out for overfitting of the model.<jupyter_code>Transformer = ColumnTransformer([
	('ohe', OneHotEncoder(), cat),
	# ('scaler', StandardScaler(), non_cat), #no need to scale on Trees
	], remainder ='passthrough')

GBR = GradientBoostingRegressor(random_state = 42)

pipeline_gbr = Pipeline([('transformer', Transformer),
						('gbr', GBR)])

#set up grid search paramater grid
gbr_param_grid = {
    'gbr__n_estimators': [50, 100, 200, 300],
    'gbr__criterion': ['friedman_mse', 'mse', 'mae'],
	'gbr__min_samples_split': [2, 3, 4, 5]

}

#set up gridsearch:
gridsearch_gbr = GridSearchCV(estimator=pipeline_gbr,
						  param_grid=gbr_param_grid,
						  cv=5)

#fit gridsearch
gridsearch_gbr.fit(X_train, y_train)

print("Best Parameter Combination Found During GBR Grid Search:")
print(gridsearch_gbr.best_params_)

Transformer = ColumnTransformer([
	('ohe', OneHotEncoder(), cat),
	# ('scaler', StandardScaler(), non_cat), #no need to scale on Trees
	], remainder ='passthrough')

GBR_best = GradientBoostingRegressor(criterion='mse', n_estimators=200, min_samples_split= 5, random_state = 42)

#pipeline
pipeline_gbr_best = Pipeline([('transformer', Transformer),
						('gbr', GBR_best)])

pipeline_gbr_best.fit(X_train, y_train)


cv = cross_validate(pipeline_gbr_best, X_train, y_train, return_train_score=True, cv=5)
cv_score = cross_val_score(pipeline_gbr_best, X_train, y_train, cv = 5)
print(f'Cross Validation Score Gradient Boosting Model: {cv}')

print(f'mean CV Score: {cv_score.mean()}')<jupyter_output>Cross Validation Score Gradient Boosting Model: {'fit_time': array([0.13171697, 0.12799406, 0.16082287, 0.14935923, 0.14530873]), 'score_time': array([0.00510597, 0.00635886, 0.00764894, 0.00540376, 0.00605416]), 'test_score': array([0.62728481, 0.68505887, 0.6747616 , 0.7939104 , 0.6750317 ]), 'train_score': array([0.96385551, 0.96685321, 0.96533828, 0.95868339, 0.96596248])}
mean CV Score: 0.6912094773793598<jupyter_text>So far Gradient Booster Model is performing the best with a mean CV score of .691. we see some overfitting with better train scores than test scores.  XGBoostFinally to try and get even better mean CV scores and model, Let's see if we can improve with X Boost. XGBoost is a powerful regressor that works off the assumption that the base learners are bad. it uses the error residuals of the previous model to fit the next model, and the final prediction is a weighted sum of all the tree predictions. This model is good for minimizing variance/overfitting.<jupyter_code>Transformer = ColumnTransformer([
	('ohe', OneHotEncoder(), cat),
	('scaler', StandardScaler(), non_cat),
	], remainder ='passthrough')

xgb_r = xg.XGBRegressor(seed = 123)

#pipeline
xboostPipe = Pipeline([('transformer', Transformer),
						('xgb_r', xgb_r)])

# objective ='reg:squarederror', n_estimators = 10, 
boost_param_grid = {
    'xgb_r__eta': [.2, .3, .4, .5, .6, .7, .8],
	'xgb_r__max_depth': [3, 4, 5, 6, 7, 8, 9, 10]
}

gridsearch_boost = GridSearchCV(estimator=xboostPipe,
						  param_grid=boost_param_grid,
						  cv=5)

#fit gridsearch
gridsearch_boost.fit(X_train, y_train)

print("Best Parameter Combination Found During XBoost Grid Search:")
print(gridsearch_boost.best_params_)
Transformer = ColumnTransformer([
	('ohe', OneHotEncoder(), cat),
	# ('scaler', StandardScaler(), non_cat),
	], remainder ='passthrough')

#Pipeline second Step: 
xgb_r_best = xg.XGBRegressor(eta = 0.2, max_depth = 3, seed = 123)

xboost_pipe_best = Pipeline([('transformer', Transformer),
						('xgb_r_2', xgb_r_best)])

xboost_pipe_best.fit(X_train, y_train)

cv = cross_validate(xboost_pipe_best, X_train, y_train, return_train_score=True, cv=5)
cv_score = cross_val_score(xboost_pipe_best, X_train, y_train, cv = 5)
print(f'Cross Validation Scores X Boost: {cv}')

print(f'mean CV Score X Boost: {cv_score.mean()}')<jupyter_output>Cross Validation Scores X Boost: {'fit_time': array([0.07625198, 0.04629803, 0.04954314, 0.05134416, 0.04434633]), 'score_time': array([0.00893188, 0.009408  , 0.00683093, 0.00692201, 0.00708604]), 'test_score': array([0.63686313, 0.7072925 , 0.71510326, 0.7883116 , 0.65902666]), 'train_score': array([0.95631625, 0.96415222, 0.96199409, 0.95670808, 0.96122295])}
mean CV Score X Boost: 0.7013194264453662<jupyter_text>X boost for resort 1 is preforming slightly better than Gradient Booster Regressor at CV score of .701  EvaulationSince XGBoost model performed the best and looks like it's the least overfit, we will evaluate the model on our test data:<jupyter_code>#boost Test Score
xboost_pipe_best.score(X_test, y_test)<jupyter_output><empty_output><jupyter_text>XGBoost gets a test CV score of .782, this means we can account for 78% of the variance in the Resort 1 data with our Gradient Booster Regressor model. Let's take a look at our predictions:<jupyter_code>y_pred_xboost = xboost_pipe_best.predict(X_test)

df = pd.DataFrame({'Real Visits Resort 1':y_test, 'Predicted XBoost Visits':y_pred_xboost, 
                    'Difference': abs(y_test - y_pred_xboost)})
print('Mean Difference, predicted - actual')
print(df['Difference'].mean())
df.head()
ohe = xboost_pipe_best.named_steps["transformer"].transformers_[0][1].get_feature_names()
numbers = non_cat

columns = [*ohe, *numbers]

for name, importance in zip(columns, xboost_pipe_best['xgb_r_2'].feature_importances_):
    print(name, importance)
feature_importance = xboost_pipe_best['xgb_r_2'].feature_importances_
sorted_idx = np.argsort(feature_importance)
pos = np.arange(sorted_idx.shape[0])
fig = plt.figure(figsize=(10, 10))
plt.subplot(1, 2, 1)
plt.barh(pos, feature_importance[sorted_idx], align="center")
plt.yticks(pos, np.array(columns)[sorted_idx])
plt.title("Feature Importance")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10795931016881176, 'fraction_numerical': 0.056550114850329965, 'mean_word_length': 5.331486611265005, 'pattern_counts': {'":': 4, '<': 77, '<?xml version=': 0, '>': 78, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '322070', 'license': 'Unlicense', 'n_tokens_mistral': 10298, 'n_tokens_neox': 8993, 'n_words': 2852, 'path': 'Resort_1_EDA_Modeling_notebook.ipynb', 'repo_name': 'lperconti/Visitation', 'size': '27436'}