DVC-ML-Experiments/Crop_Pred_ Nigeria.py

#!/usr/bin/env python
# coding: utf-8

# # Crop  Prediction for Nigeria Agriculture using Essmbel Machine learning
# 
# Crop yield which is too uncertain to estimate is an esential factor for determining farmer's income level. The difficult to estimate yiled of agricultural crop is duetto its reliance on weather conditions (rain, temperature, etc), pesticides and information about history of crop yield even. These all makes possibel making decisions related to agricultural risk management and future predictions if obtained and articulate exactly. Dueto its yet not innovated sector and its natural behaviour obtaing concise and  organized agricultural data even today is very challenging. One important source of information in this sector thta provide a panel data is Food and agricultural organization(FAO).From FAO Database data related to crop yiled, rainfall and pestcide extracted fot Nigeria.  To make the climatic factor complet from abiotic factor perspectiv, we  also considere nigeria's temperature to incorporate with crop yiled prediction.

# # Table of contents 
# 
# 1. Data Importing and Prearation
# 
# 2. Exploratory Analysis
# 
# 3. Model Developmemnt and Selection
# 
# 6. Hybriding algorithm  with pipline
# 
# 7. Drawing conclusion

# # 1. Data Importing and preparation

# As the outpuet indicates there are about 56717 rows and 12 columns for the crop yiled dataset. only column like Area, Itme and value columns are importnat  to our analsysis. To extract yield data for nigeria we can use  Area column where  the country can be extracted wit its name Nigeria uing the function `(df_yield.loc[df_yield['Area'] == 'Nigeria'])`. It is also possible to exttract using Area Code column and `(df_yield.loc[df_yield['Area Code'] == '159'])` Thanks to William Green from the Omdena Zowasel chalenege for sharng Nigerian agriculture dataset from FAO (https://files.slack.com/files-pri/T03NM44BRB9-F03TAHA5FBM/download/crop_yield-20220811t144327z-001.zip?origin_team=T03NM44BRB9)

# In[1]:


Nig_yield = pd.read_csv('C:\\Users\\moga\\Desktop\\Files\\Crop_dataset\\crop_yield\\FAOSTAT_data_en_8-9-2022_yield.csv')
Nig_yield.head()


# In[2]:


# Crop Nierogen dataset
Nig_Nitro = pd.read_csv('C:\\Users\\moga\\Desktop\\Files\\Crop_dataset\\crop_yield\\FAOSTAT_data_en_8-9-2022_nitrogen.csv')
Nig_Nitro.head()


# In[3]:


# import  Rainfal dataset
Nig_rain = pd.read_csv('C:\\Users\\moga\\Desktop\\Files\\Crop_dataset\\crop_yield\\1901-2021_NGA_rain.csv')
Nig_rain .head()


# In[4]:


#import temprature dataset
Nig_temp = pd.read_csv('C:\\Users\\moga\\Desktop\\Files\\Crop_dataset\\crop_yield\\observed-average-annual-mean-temperature-of-nigeria-for-1901-2021.csv')
Nig_temp.head()


# ## 1.2 Checking for Missing value

# In[5]:


#check for missikng value
Nig_yield.info()


# In[6]:


Nig_yield.describe()


# In[7]:


#check missing data for nitrogen 
Nig_Nitro.info()


# In[8]:


#  descriptive statistics for niterogen data
Nig_Nitro.describe()


# In[9]:


Nig_rain.info()


# In[10]:


#Dropunwanted column for raainfall
Nig_rain.describe()


# In[11]:


Nig_temp.info()


# In[12]:


Nig_temp.describe()


# ## 1.2 Droping unwanted Columns

# In[13]:


Nig_yield.columns


# In[14]:


# drop unwanted columns for yiled data
Nig_yield = Nig_yield.drop(['Domain Code','Area Code (FAO)','Area','Element Code',
                            'Item Code (FAO)','Year Code','Unit','Flag','Flag Description'], axis=1)
Nig_yield.head()


# In[15]:


#Drop unwanted column for niterogen data
Nig_Nitro=Nig_Nitro.drop(['Domain Code','Domain','Area Code (FAO)','Area','Element Code','Element',
                          'Item Code','Item','Year Code','Unit','Flag','Flag Description'],axis=1)
Nig_Nitro.head()


#  Though the datasetsfor the four features are cleaned and checked for missing value, they are in separate dataframe and mereging into one datafrmae is done next.
# 
# ## 1.3  Combining DataSets
#  Three possible ways or functions to combine data in panda includes
#  - **Merege functio,** `Mereg()`:for combining data on common columns or indices
#  
#  - **Join Functio**`.join()`:for combining data on common columns or indices
#  
#  - **concatnate function**, `Concat()`: for combining DataFrames across rows or columns

# In[16]:


yile_nitr = pd.merge(Nig_yield,Nig_Nitro, on=["Year"] )
yile_nitr.head()


# In[17]:


yile_nitr_rain = pd.merge(yile_nitr,Nig_rain, on =["Year"] )
yile_nitr_rain.head()


# Since column name for temprature dataset is`Category` it must be renamed as `Year`

# In[18]:


Nig_temp=Nig_temp.rename(index=str, columns={"Category": 'Year'})
Nig_temp.head()


# In[19]:


# defining the combibed dataframe
df_com = pd.merge(yile_nitr_rain,Nig_temp, on =["Year"] )
df_com


# In[20]:


df_com.columns


# In[21]:


# still droping of some column seems essential 
df_com = df_com.drop(['Domain','Element'], axis =1)
df_com 


#  - let's also rename some columns for convinece

# In[22]:


df_com = df_com.rename(index=str, columns={"Value_x": 'Yiled/ha',"Value_y":'Niro/ha',
                                           "Annual Mean":'Mean_Temp',"5-yr smooth":'smooth_Temp'})
df_com.head()


# In[23]:


df_com.describe()


# # 2. EXploratory Analysis

# In[24]:


from mpl_toolkits.mplot3d import Axes3D
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt # plotting
import seaborn as sns
import numpy as np # linear algebra
import os # accessing directory structure
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)from mpl_toolkits.mplot3d import Axes3D
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt # plotting
import numpy as np # linear algebra
import os # accessing directory structure
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_c


# # 2.1 Exploring Data Aggregation uing Seaborn Visualization

# In[25]:


#grouping  dataset by column Item for production information
df_com.groupby(['Item'],sort=True)['Yiled/ha'].sum().nlargest(20)


# In[26]:


#aggregation of the data is based on the various momnets and respective variance 
plt.figure(figsize=(50,30))
sns.set_context('poster',font_scale=1.6)
plt.subplot(3,1,1)
sns.barplot(x='Item',y='Value', data=Nig_yield)
plt.legend(['mean aggregate'],loc=1)
plt.xticks(rotation=90)
plt.subplot(3,1,2)
sns.barplot(x='Item',y='Value', data=Nig_yield, estimator =np.median)
plt.legend(['medain aggregate'],loc=1)
plt.xticks(rotation=90)
plt.subplot(3,1,3)
sns.barplot(x='Item',y='Value', data=Nig_yield, estimator =np.std)
plt.legend(['std aggregate'],loc=1)
plt.xticks(rotation=90)


# In[27]:


plt.figure(figsize=(50,30))
sns.set_context('paper',font_scale=4,)
plt.subplot(3,1,1)
sns.barplot(x='Item',y='Value', data=Nig_yield, estimator =np.var)
plt.legend(['var aggregate'],loc=1)
plt.xticks(rotation=90)
plt.subplot(3,1,2)
sns.barplot(x='Item',y='Value', data=Nig_yield, estimator =np.cov)
plt.legend(['Cov aggregate'],loc=1)
plt.xticks(rotation=90)


# ## 2.2  Pivoting and Ternd Analsysi

# In[28]:


get_ipython().run_line_magic('matplotlib', 'inline')
plt.figure(figsize=(30,25))
sns.set()
table = pd.pivot_table(Nig_yield, values = 'Value', index=['Year'],
                    columns=['Item'], aggfunc=np.sum)
table.plot()


# In[29]:


get_ipython().run_line_magic('matplotlib', 'inline')
plt.figure(figsize=(30,25))
sns.set()
Nig_yield.pivot_table('Value', index=['Year'],
                    columns=['Item'], aggfunc=np.sum).plot()
plt.ylabel('yield per year')


# In[30]:


sns.pairplot(Nig_yield)


# In[31]:


# use hue to categorize base on categorical vriable  sex for example
sns.pairplot(Nig_yield,hue='Item', palette='Blues')


# In[32]:


plt.figure(figsize=(50,30))
sns.set_theme(color_codes=True)
#sns.set_context('paper',font_scale=1.4)
sns.lmplot(x='Mean_Temp',y='Yiled/ha', hue='Item',data=df_com)


# In[33]:


plt.figure(figsize=(50,30))
sns.set_theme(color_codes=True)
#sns.set_context('paper',font_scale=1.4)
sns.lmplot(x='smooth_Temp',y='Yiled/ha', hue='Item',data=df_com)


# In[34]:


plt.figure(figsize=(50,30))
sns.set_theme(color_codes=True)
sns.set_context('paper',font_scale=1.4)
sns.lmplot(x='Niro/ha',y='Yiled/ha', hue='Item',data=df_com)


# In[35]:


plt.figure(figsize=(50,30))
sns.set_theme(color_codes=True)
sns.set_context('paper',font_scale=1.4)
sns.lmplot(x='Year',y='Yiled/ha', hue='Item',data=df_com)


# In[36]:


#option 2: using heatmap
plt.figure(figsize=(8,6))
sns.set_context('paper',font_scale =1.4)
yield_mx= df_com.corr()
sns.heatmap(yield_mx,annot =True, cmap='Blues')


# ## 2.4 curse of symetery
# 
# In <font color ='magneta'>Cramer’s_V</font> it seems possible to guarantee what value of the feature for known value of the response variable but is difficult for the revers. This is known as <font color='magneta'>symetry property</font> and works for <font color ='magneta'>Cramer’s V</font>.
# 
# For utilizing asymmetric measure of association between categorical features,<font color='magneta'>Theil’s U</font> is being recommnded and also referred to as the <font color='magneta'>Uncertainty Coefficient</font>.
# 
# It is based on the <font color ='magneta'>conditional entropy</font> between x and y — or in human language, given the value of x, how many possible states does y have, and how often do they occur.
# Just like<font color='magneta'>Cramer’s V</font>, the output value is on the range of [0,1], with the same interpretations as before — but unlike <font color='magneta'> Cramer’s V</font>, it is asymmetric, meaning U(x,y)$\ne$ U(y,x) (while V(x,y)=V(y,x), where V is <font color='magneta'>Cramer’s V</font>. Using<font color='magneta'>Theil’s U</font> in the simple case above will let us find out that knowing y means we know x, but not vice-versa. <font color='magneta'>Theil's U</font>, also known as the<font color ='magneta'>Uncertainty Coefficient</font>. 
# 
#  Formaly marked as U(x|y), this coefficient provides a value in the range of [0,1], where 0 means that feature y provides no information about feature x, and 1 means that feature y provides full information abpout features x's value.
# 
# Theil’s U statistic is a relative accuracy measure that compares the forecasted results with the results of forecasting with minimal historical data. It also squares the deviations to give more weight to large errors and to exaggerate errors, which can help eliminate methods with large errors
# 
# The formula for calculating Theil’s U statistic:
# ![image.png](attachment:image.png)
# 
# where $Y_t$ is the actual value of a point for a given time period t, n is the number of data points, and $\hat{Y_t}$ is predicted valye.

# In[37]:


import math
from collections import Counter
import numpy as np
import seaborn as sns
import pandas as pd
import scipy.stats as ss
import matplotlib.pyplot as plt
import sklearn.preprocessing as sp
from sklearn.tree import DecisionTreeClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report, confusion_matrix, accuracy_score
from subprocess import check_output


# In[38]:


def conditional_entropy(x,y):
    # entropy of x given y
    y_counter = Counter(y)
    xy_counter = Counter(list(zip(x,y)))
    total_occurrences = sum(y_counter.values())
    entropy = 0
    for xy in xy_counter.keys():
        p_xy = xy_counter[xy] / total_occurrences
        p_y = y_counter[xy[1]] / total_occurrences
        entropy += p_xy * math.log(p_y/p_xy)
    return entropy

def theil_u(x,y):
    s_xy = conditional_entropy(x,y)
    x_counter = Counter(x)
    total_occurrences = sum(x_counter.values())
    p_x = list(map(lambda n: n/total_occurrences, x_counter.values()))
    s_x = ss.entropy(p_x)
    if s_x == 0:
        return 1
    else:
        return (s_x - s_xy) / s_x


# In[39]:


theilu = pd.DataFrame(index=['Item'],columns=df_com.columns)
columns = df_com.columns
for j in range(0,len(columns)):
    u = theil_u(df_com['Item'].tolist(),df_com[columns[j]].tolist())
    theilu.loc[:,columns[j]] = u
theilu.fillna(value=np.nan,inplace=True)
plt.figure(figsize=(20,1))
sns.heatmap(theilu,annot=True,fmt='.2f')
plt.show()


# # 3.Prediction Model Development and Selection
# 
# - Here,we will explore both classification and Regression algorithms of the supervised machine learning. For the classification model dvelopement those categorical columns (Item) used as target variables whereas the regration models development constructed using anaual yiled as target variable.

# In[40]:


# Column Varaible is categorical
df_com.info()


# ####  3.1  Baseline Setting
# 
# While forecasting there are available accuracy measures with there own prons and cons but generally gropuded as (i) scale dependent errors, (ii) percentage errors and (iii) scaled errors[<font color='blue'>vorgelegt von,2019</font>] the first is for <font color ='Magneat'>mean absolute error</font> (<font color='pureple'>MAE</font>) and <font color ='Magneat'>root mean square error</font> (<font color='pureple'>RMSE</font>) which they typically used to compare dataset having same units and minimizaing them respectively leads to median and mean of the distribution.  
# 
# <font color ='Mgneat'>mean absolute percentage error </font> (<font color='pureple'> MAPE</font>) and <font color='Mgneat'>symmetric mean absolute percentage error</font> (<font color='pureple'>SMAP</font>) are the Percentage errors tyeps which aallows to measure time seriesat different scale while <font color='Magneat'>mean percentage error</font>(<font color='pureple'>MPE</font>) in this category  help as a bias indicator as negative and positive errors offset each other.
# 
# Scaled errors as an alternative to percentage-based errors for comparing forecasts on datasets having different scales or units. It incudes the <font color='Magneat'>mean absolute scaled error</font> (<font color='pureple'>MASE</font>).
# 
# Let $X_u$ stands for uneconded features and the target varial y which is demand in our case has shown no change at all with regard to encoding as it is contineuous numeric.

# In[41]:


import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.metrics import mean_absolute_error, make_scorer 
from sklearn.metrics import mean_squared_error
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.impute import SimpleImputer
from sklearn.model_selection import KFold 
from sklearn.model_selection import train_test_split
from sklearn.model_selection import cross_val_score
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.linear_model import BayesianRidge
from scipy.special import comb
from xgboost import XGBRegressor

np.random.seed(12345)

get_ipython().run_line_magic('matplotlib', 'inline')
get_ipython().run_line_magic('config', "InlineBackend.figure_format = 'svg'")


# In[42]:


Xu = df_com.drop(columns=['Yiled/ha'],axis=1)
yu = df_com['Yiled/ha']


# In[43]:


Xu_train, Xu_test, yu_train, yu_test = train_test_split(Xu, yu, test_size=0.2, random_state=42)
mean_absolute_error(yu_train, 
                    np.full(yu_train.shape[0], yu_train.mean()))


# In[44]:


mean_squared_error(yu_train, 
                    np.full(yu_train.shape[0], yu_train.mean()))


# Hence, the encoding method is expected to  make ready the data for our algorthim  that will give <font color ='pureple'>MAE</font> below 38021.27671310749. Since the prediction model will not be perfect a kind of nois which is random must be acknowledged but with normally distributed having a standard deviation of $\sigma^2$=1 and therefore whatever our model is perefect we have to expect <font color ='pureple'>MAE</font>=0.7868967904399754

# In[45]:


mean_absolute_error(np.random.randn(1000), 
                    np.zeros(1000))


# ##  3.2 Crop Prediction Using KNN

# In[46]:


from mpl_toolkits.mplot3d import Axes3D
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt # plotting
import numpy as np # linear algebra
import os # accessing directory structure
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)from mpl_toolkits.mplot3d import Axes3D
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt # plotting
import numpy as np # linear algebra
import os # accessing directory structure
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)


# In[47]:


#define functions for plotting data.
# Distribution graphs (histogram/bar graph) of column data
def plotPerColumnDistribution(df, nGraphShown, nGraphPerRow):
    nunique = df.nunique()
    df = df[[col for col in df if nunique[col] > 1 and nunique[col] < 50]] # For displaying purposes, pick columns that have between 1 and 50 unique values
    nRow, nCol = df.shape
    columnNames = list(df)
    nGraphRow = (nCol + nGraphPerRow - 1) / nGraphPerRow
    plt.figure(num = None, figsize = (6 * nGraphPerRow, 8 * nGraphRow), dpi = 80, facecolor = 'w', edgecolor = 'k')
    for i in range(min(nCol, nGraphShown)):
        plt.subplot(nGraphRow, nGraphPerRow, i + 1)
        columnDf = df.iloc[:, i]
        if (not np.issubdtype(type(columnDf.iloc[0]), np.number)):
            valueCounts = columnDf.value_counts()
            valueCounts.plot.bar()
        else:
            columnDf.hist()
        plt.ylabel('counts')
        plt.xticks(rotation = 90)
        plt.title(f'{columnNames[i]} (column {i})')
    plt.tight_layout(pad = 1.0, w_pad = 1.0, h_pad = 1.0)
    plt.show()


# In[48]:


# Correlation matrix
def plotCorrelationMatrix(df, graphWidth):
    filename = df.dataframeName
    df = df.dropna('columns') # drop columns with NaN
    df = df[[col for col in df if df[col].nunique() > 1]] # keep columns where there are more than 1 unique values
    if df.shape[1] < 2:
        print(f'No correlation plots shown: The number of non-NaN or constant columns ({df.shape[1]}) is less than 2')
        return
    corr = df.corr()
    plt.figure(num=None, figsize=(graphWidth, graphWidth), dpi=80, facecolor='w', edgecolor='k')
    corrMat = plt.matshow(corr, fignum = 1)
    plt.xticks(range(len(corr.columns)), corr.columns, rotation=90)
    plt.yticks(range(len(corr.columns)), corr.columns)
    plt.gca().xaxis.tick_bottom()
    plt.colorbar(corrMat)
    plt.title(f'Correlation Matrix for {filename}', fontsize=15)
    plt.show()


# In[49]:


# Scatter and density plots
def plotScatterMatrix(df, plotSize, textSize):
    df = df.select_dtypes(include =[np.number]) # keep only numerical columns
    # Remove rows and columns that would lead to df being singular
    df = df.dropna('columns')
    df = df[[col for col in df if df[col].nunique() > 1]] # keep columns where there are more than 1 unique values
    columnNames = list(df)
    if len(columnNames) > 10: # reduce the number of columns for matrix inversion of kernel density plots
        columnNames = columnNames[:10]
    df = df[columnNames]
    ax = pd.plotting.scatter_matrix(df, alpha=0.75, figsize=[plotSize, plotSize], diagonal='kde')
    corrs = df.corr().values
    for i, j in zip(*plt.np.triu_indices_from(ax, k = 1)):
        ax[i, j].annotate('Corr. coef = %.3f' % corrs[i, j], (0.8, 0.2), xycoords='axes fraction', ha='center', va='center', size=textSize)
    plt.suptitle('Scatter and Density Plot')
    plt.show()


# In[50]:


correlation_data=df_com.select_dtypes(include=[np.number]).corr()

mask = np.zeros_like(correlation_data, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True

f, ax = plt.subplots(figsize=(11, 9))

# Generate a custom diverging colormap
cmap = sns.palette="vlag"

# Draw the heatmap with the mask and correct aspect ratio
sns.heatmap(correlation_data, mask=mask, cmap=cmap, vmax=.3, center=0,
            square=True, linewidths=.5, cbar_kws={"shrink": .5});


# In[51]:


import numpy as np
import pandas as pd
from sklearn.neighbors import KNeighborsClassifier
from matplotlib import pyplot as plt
from scipy.interpolate import make_interp_spline
#from pandas_profiling import ProfileReport
from pylab import rcParams
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import accuracy_score as acc
from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay
import seaborn as sns

import warnings
warnings.filterwarnings('ignore')


# For the classification model developemnet we prefer to drop yield column of the combined dataframe as it is a result than a factor to influnce itself. However, we run it for checking on classification purpose and it outperforms the KNN model that ignors yield column we are not sure why it is so if any one can to justify. 

# In[52]:


df= df_com.drop(['Yiled/ha'],axis=1)
df.head() 


# In[53]:


column_names = ['Year','Rain','Mean_Temp','smooth_Temp','Item']
df = df.reindex(columns=column_names)
df


# In[54]:


col = list(df.columns)
classes = df["Item"].unique()


# In[55]:


xdata = df.iloc[:, 0:4].values
ydata = df.iloc[:, 4].values
plt.figure(figsize = (16, 9))
sns.countplot(classes, palette = 'rocket')
plt.xticks(rotation=90)


# **Plotting the various values of the input parameter for a particular output value**

# In[56]:


all_col = df.columns[:-1]

for col in all_col:
  plt.figure(figsize = (16, 9))
  sns.barplot(x = 'Item', y = col, data = df, palette = 'rocket')
  plt.xlabel('Item', fontsize = 12)
  plt.ylabel(col, fontsize = 12)
  plt.xticks(rotation=90)
  plt.title(f'{col} vs Crop')
  plt.show()


# In[57]:


plt.figure(figsize = (10, 17))
sns.pairplot(df, hue = 'Item', palette = 'rocket')
plt.show()


# **Split & scale the data into train and test dataset**
# 
# The data set is splited into test and train dataset based on the 80 - 20 rule and it is scaled for better performance for training.

# In[58]:


from sklearn.linear_model import LinearRegression, SGDRegressor
from sklearn.ensemble import RandomForestRegressor
from xgboost import XGBClassifier, XGBRegressor
from sklearn.model_selection import train_test_split
from imblearn.datasets import make_imbalance
from category_encoders.target_encoder import TargetEncoder
import statsmodels.api as sm
xtrain, xtest, ytrain, ytest = train_test_split(xdata, ydata, test_size=0.2, random_state=884)
x_st = StandardScaler()
xtrain = x_st.fit_transform(xtrain)
xtest = x_st.fit_transform(xtest)


# **Requied only for optimising the performance of the KNN algorithm**
# 
# The training model is being optimised using trial and error method to determine the best random_state in test_train_split function and the nearest neighbor in the KNN algorithm.

# In[59]:


acc_list = []
err_rate = []

neighbors = np.linspace(1, 50, 50)
neighbors = neighbors.astype(int)

for K in neighbors:
  classifier = KNeighborsClassifier(n_neighbors = K)
  classifier.fit(xtrain, ytrain)
  y_pred = classifier.predict(xtest)

  accuracy = round(acc(ytest, y_pred)*100, 3)

  acc_list.append(accuracy)
  err_rate.append(np.mean(y_pred != ytest))

xy = make_interp_spline(neighbors, acc_list)
xz = make_interp_spline(neighbors, err_rate)
x = np.linspace(1, 50, 1000)
y = xy(x)
z = xz(x)

plt.figure(figsize = (13, 7))
plt.subplot(2, 1, 1)
sns.lineplot(x, y, linewidth = 2, color = '#5C284F')
plt.xlabel('K value')
plt.ylabel('Accuracy (%)')
plt.title('Accuracy vs K', fontweight = 'bold')
plt.xlim(min(neighbors), max(neighbors))

plt.subplot(2, 1, 2)
sns.lineplot(x, z, linewidth = 2, color = '#D96856')
plt.xlabel('K value')
plt.ylabel('Loss')
plt.title('Loss vs K', fontweight = 'bold')
plt.xlim(min(neighbors), max(neighbors))

plt.tight_layout()
plt.show()

K_opt = acc_list.index(max(acc_list))
print('\nOptimal value of K = ', K_opt)


# **K Nearest Neighbor Algorithm**
# 
# The classification is done based on the k nearest neighbour algorithm where the euclidian distance is calculated for each input output combination and the classification is done depending on the minimum euclidian distance.

# In[60]:


classifier = KNeighborsClassifier(n_neighbors=K_opt+1)
classifier.fit(xtrain, ytrain)
y_pred = classifier.predict(xtest)

accuracy = acc(ytest, y_pred)*100
print('Accuracy of the training Model : ', round(accuracy, 3), '%')


# **Thi is Too poor accuracy and needs further explortion, i don't know where to start actually.**

# **Display the Performance of the trained model**
# 
# A confusion matrix is used to determine the performance of a trained model by mapping the training accuracy for each input and output combination. Ideally it should be a diagonal matrix but due to training uncertainity , the matrix obtained is not a diagonal matrix.

# In[61]:


cm = confusion_matrix(ytest, y_pred, normalize = 'pred')

fig, ax = plt.subplots(figsize=(25,15))
sns.heatmap(cm, annot = True)
plt.xlabel('Predicted Crop', fontsize = 12)
plt.ylabel('Actual Crop', fontsize = 12)
plt.title('Confusion Matrix', fontweight = 'bold', fontsize = 15)

plt.xticks(rotation=90)
plt.yticks(rotation=0)

ax.xaxis.set_ticklabels(classes)
ax.yaxis.set_ticklabels(classes)
plt.show()


# # 3.3 Crop Prediction model Regresion approach 
# 
# Another important supervised learing  alogorithm is the Regression approach by using anaual yield as target variable which is numerical than Items that are categorical.
# In our classification based prediction to the dataset,perfomance obtained is too poor even after incorporating the yield column of orginal dataframe into the analysis (**accuracy become 3.4%**). 
# 
# let's turn into another supervised alogorithm,regression,now since peritty of having with numerical target variable that is yield. We lend oursveles here to further data processing since we do have a categorical variable, Item column.
# 
# Unlike label encoding and One-Hot encoding that respectively experinced value implication and dimenionality problems.

# In[62]:


import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
import seaborn as sns
from scipy import stats
from sklearn.preprocessing import *
from sklearn. metrics import *
from sklearn.linear_model import LinearRegression, SGDRegressor
from sklearn.ensemble import RandomForestRegressor
from xgboost import XGBClassifier, XGBRegressor
from sklearn.model_selection import train_test_split
from imblearn.datasets import make_imbalance
from category_encoders.target_encoder import TargetEncoder
import statsmodels.api as sm


# In[63]:


df.info()


# For a detail to Label and One Hot encoding read https://www.analyticsvidhya.com/blog/2020/03/one-hot-encoding-vs-label-encoding-using-scikit-learn/
# 
# One-Hot encoding is the process of creating dummy variables.

# - One-Hot Encoding results in a Dummy Variable Trap as the outcome of one variable can easily be predicted with the help of the remaining variables.
# 
# - Dummy Variable Trap is a scenario in which variables are highly correlated to each other.
# 
# - The Dummy Variable Trap leads to the problem known as **multicollinearity**. 
# 
# - Multicollinearity occurs where there is a dependency between the independent features.
# 
# - It is a serious issue in machine learning models like `Linear Regression` and `Logistic Regression`.
# 
# So, in order to overcome the problem of multicollinearity, one of the dummy variables has to be dropped. One of the common ways to check for multicollinearity is the Variance Inflation Factor (VIF):
# 
# - VIF=1, Very Less Multicollinearity
# 
# - VIF<5, Moderate Multicollinearity
# 
# - VIF>5, Extreme Multicollinearity (This is what we have to avoid)
# Compute the VIF scores:

# In order not traped by target leakage due to the underlings probablity based target encoding we deploye Target encoder with prior smoothing by Vinícius Trevisan (https://towardsdatascience.com/dealing-with-categorical-variables-by-using-target-encoder-a0f1733a4c69) and https://towardsdatascience.com/categorical-feature-encoding-547707acf4e5

# ### 3.3.1 One Hot Encoding
# 
# one-Hot encoding some times refered as dummies

# In[64]:


from pandas import DataFrame,get_dummies
df_com.columns = df_com.columns.to_series().apply(lambda x: x.strip())
data_to_encode = df_com[['Item']]
df_com_dumy = get_dummies(df_com, prefix={k:"dmy_%s"%k for k in data_to_encode},
                                 columns = list(data_to_encode))
df_com_dumy.head (5)


# In[65]:


# Compare sizes
print('Original size:', df_com.shape)
print('One-hot encoded size:', df_com_dumy.shape)


# It is observed that the dataset now changed to a shape of 2790x54 which was 2790x7 in the orginal dataset.

# # 3.3.2 Evaluate Onehot encoding

# Now lets evaluate one-hot encoding using
# - <font color='magneta'>Linear regression</font>(<font color='pureple'>lr</font>)
# 
# - <font color='Magneta'>Suport vector Regression</font>,
# 
# - (<font color='pureple'>SVR</font>) and
# 
# - <font color='Magnea'>Random Forest</font>(<font color='pureple'>RF</font>) algorithm

# Let $X_{oh}$ stands for feature variable affter one-hot encding and $X_t$ for target encoding and the response variable remains constant as it is free from encoding. 

# In[66]:


Xoh = df_com_dumy.drop(columns=['Yiled/ha'],axis=1)
y = df_com_dumy['Yiled/ha']


# In[67]:


# convvert dataframe to array to make the ML ease
#array  = df_com_dumy.values
#Xoh = array[:,1:54]
#y = array[:,1]
#Xoh.shape


# In[68]:


Xoh_train, Xoh_test, y_train, y_test = train_test_split(Xoh, y, test_size=0.3, random_state=42)
print ("training dataset", Xoh_train,y_train)
print("testing datasdet", Xoh_test,y_test)


# In[69]:


#import libraries 
from pandas import Series,DataFrame
from sklearn.preprocessing import OneHotEncoder
from sklearn.linear_model import LinearRegression
#train the model
lr_oh = LinearRegression()
lr_oh.fit(Xoh_train,y_train)
#predict ton the test data
pred_oh =lr_oh.predict(Xoh_test)
# error determination
errors_lr_oh =abs(pred_oh -y_test)
mse =np.mean((pred_oh -y_test)**2)
#print out mean absolute error(mae)
print('Mean Absolute Error(MAE):',round(np.mean(errors_lr_oh),2),'yield/ha')
print('Mean Square Error (MSE):',round(np.mean(mse),2),'yiled/ha')


# An improvement about 25% on MAE (9845.51/38021.27671310749) and 1.8% (57039460.81/3162915966.08407)on MSE obtained
# lets try another meterics for prediction model, i.e., coefficient of determination <font color='Magneat'> $R^2$</font> and <font color='Magneat'>adjusted-$R^2$</font>

# In[70]:


#model score
lr_oh.score(Xoh_test,y_test)


# Determining the performance metrics  for a model is ofcourse to perform accuracy evaluation  using <font color='blue'> mean average percentage errore</font> (<font color='purple'>MAPE</font> subtracted from 100% 

# In[71]:


# calculate MAPE
mape_lr_oh=100*(errors_lr_oh/y_test)
#calculate and display accuracy 
accuracy =100-np.mean(mape_lr_oh)
print('Accuracy:',round(accuracy,2),'%')


# In[72]:


_plot=plt.scatter(pred_oh,(pred_oh-y_test),c='b')
plt.hlines(y=0,xmin=-0.4,xmax=3)
plt.title('Residual plot')


# In[73]:


# get importance 
from matplotlib import pyplot
importance = lr_oh.coef_
# summarize model importance

for i, v in enumerate(importance):
    print('Feature:%0d, score:%.5f'%(i,v))
#plot feature importance 
pyplot.bar([x for x in range(len(importance))],importance )
plt.xlabel('Features')
plt.ylabel('features importance')
pyplot.show()


# In[80]:


predictors = df_com_dumy.drop(['Yiled/ha'],axis=1).columns
coef = Series(lr_oh.coef_,predictors).sort_values()
coef.plot(kind='bar',title='Model Coefficients')


# In[ ]:


def one_hot_encoder_one(data,Item,keep_first=True):
    oh = OneHotEncoder()
    oh_df = pd.DataFrame(oh.fit_transform(df_com[[Item]]).toarray())
    oh_df.columns = oh.get_feature_names()
    for col in oh_df.columns:
        oh_df.rename({col:f'{Item}_'+col.split('_')[1]},axis=1,inplace=True)
    new_data = pd.concat([df_com,oh_df],axis=1)
    new_data.drop(Item,axis=1,inplace=True)
    if keep_first == False:
        new_data=new_data.iloc[:,1:]
    return new_data


# combine the data and potentially remove one column for auto correlation

# In[ ]:


def one_hot_encoder_two(data,Fetaure,keep_first=True):

    one_hot_cols = pd.get_dummies(data[feature])
    
    for col in one_hot_cols.columns:
        one_hot_cols.rename({col:f'{Item}_'+col},axis=1,inplace=True)
    
    new_data = pd.concat([df_com,one_hot_cols],axis=1)
    new_data.drop(Item,axis=1,inplace=True)
    
    if keep_first == False:
        new_data=new_data.iloc[:,1:]
    
    return new_data


# In[ ]:


df_com.head()


# In[ ]:


one_hot_encoder_one(df_com,'Item')


# In[ ]:


stats = df_com['Yiled/ha'].groupby(df_com['Item']).agg(['count', 'mean'])


# In[ ]:


smoothing_factor = 1.0 # The f of the smoothing factor equation 
min_samples_leaf = 1 # The k of the smoothing factor equation
prior = df_com['Yiled/ha'].mean()
smoove = 1 / (1 + np.exp(-(stats['count'] - min_samples_leaf) / smoothing_factor))
smoothing = prior * (1 - smoove) + stats['mean'] * smoove
encoded = pd.Series(smoothing, name = 'genre_encoded_complete')


# In[ ]:


from category_encoders import TargetEncoder
encoder = TargetEncoder()
df_com['Item_encoded'] = encoder.fit_transform(df_com['Item'], df_com['Yiled/ha'])


# In[ ]:


df_com['Item_encoded'] 


# The above target encoding is based on binary classification whiel the dataframe used here is a multiclass problem and must be encodedd using `category_encoders.TargetEncoder` object in this scenario:

# In[ ]:


from category_encoders import TargetEncoder
targets = df_com['Yiled/ha'].unique()
for t in targets:
    target_aux = df_com['Yiled/ha'].apply(lambda x: 1 if x == t else 0)
    encoder = TargetEncoder()
    df_com['Item_encoded_sklearn_target_' + str(t)] = encoder.fit_transform(df_com['Item'], target_aux)


# In[ ]:


te_df=df_com.copy()
for col in te_df.select_dtypes(include='O').columns:
    te=TargetEncoder()
    te_df[col]=te.fit_transform(te_df[col],te_df['Yiled/ha'])


# In[ ]:


te_df.head()


# Just like label and ordinal encoding, we lose the name of the actual values for each particular feature

# **Building** `target_encoding function`
# 
# - Inputs are data, a feature, and the target feature
# 
# - A new data frame is created containing each unique value of a feature from the data which is then grouped by it’s mean target value
# 
# - An empty dictionary is then created, filled with this data, and mapped to each unique value of the particular feature
# 
# - The return value is the data frame with new target encodings in place

# In[ ]:


def target_encoding(data, column, target):
    
    grouped = data[[column,target]].groupby(column,as_index=False).mean()
    empty_dict = {}
    for i in range(len(grouped)):
        empty_dict[grouped.iloc[i,0]]=grouped.iloc[i,1]
    data[column]=data[column].map(lambda x: empty_dict[x])
    
    return data


# - Let’s see the same result with the new function

# In[ ]:


te_df=df_com.copy()
for col in te_df.select_dtypes(include='O').columns:
    target_encoding(te_df,col,'Yiled/ha')


# - ordinal data

# In[ ]:


df_com.head()


# - New data

# In[ ]:


te_df.head()


# In[ ]:


def reg_model(data):
    X = df_com.drop('Yiled/ha',axis=1)
    y = df_com['Yiled/ha']
    X_tr,X_te,y_tr,y_te = train_test_split(X,y,random_state=14,test_size=0.25)
    linreg = LinearRegression()
    linreg.fit(X_tr,y_tr)
    p = linreg.predict(X_te)
    print(f'R-squared: {r2_score(y_te,p)}')
    print('-'*20)
    print(f'Error: {mean_absolute_error(y_te,p)}')
    coefs = pd.DataFrame(linreg.coef_).T
    coefs.columns=X.columns
    print('-'*20)
    return coefs


# In[ ]:


reg_model(te_df)


# ![image-2.png](attachment:image-2.png)

# **2. Crop recomendation Essenble learning**
# credti to `Prakash Kumar.K`

# In[ ]:


import pandas as pd
import pandas_profiling as pp
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from jinja2.utils import escape
import warnings
import os
import plotly.graph_objects as go
import plotly.io as pio
import pickle
from sklearn.utils import resample
# Metrics
from sklearn.metrics import accuracy_score, classification_report, confusion_matrix, auc, roc_curve

# Validation
from sklearn.model_selection import train_test_split, cross_val_score, KFold
from sklearn.pipeline import Pipeline, make_pipeline

# Tuning
from sklearn.model_selection import GridSearchCV

# Feature Extraction
from sklearn.feature_selection import RFE

# Preprocessing
from sklearn.preprocessing import MinMaxScaler, StandardScaler, Normalizer, Binarizer, LabelEncoder

# Models
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.linear_model import LogisticRegression
from sklearn.naive_bayes import GaussianNB
from sklearn.svm import SVC
from sklearn.neighbors import KNeighborsClassifier
from sklearn.tree import DecisionTreeClassifier

# Ensembles
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import BaggingClassifier
from sklearn.ensemble import AdaBoostClassifier
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.ensemble import ExtraTreesClassifier

warnings.filterwarnings('ignore')


# In[ ]:


sns.set_style("whitegrid", {'axes.grid' : False})
pio.templates.default = "plotly_white"

# Analyze Data
def explore_data(df):
    print("Number of Instances and Attributes:", df.shape)
    print('\n')
    print('Dataset columns:',df.columns)
    print('\n')
    print('Data types of each columns: ', df.info())


# In[ ]:


# Checking for duplicates
def checking_removing_duplicates(df):
    count_dups = df.duplicated().sum()
    print("Number of Duplicates: ", count_dups)
    if count_dups >= 1:
        df.drop_duplicates(inplace=True)
        print('Duplicate values removed!')
    else:
        print('No Duplicate values')


# In[ ]:


# Split training and validation set
def read_in_and_split_data(data, target):
    X = data.drop(target, axis=1)
    y = data[target]
    X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.2, random_state=0)
    return X_train, X_test, y_train, y_test


# In[ ]:


# appending Algorithms
def GetModel():
    Models = []
    Models.append(('LR'   , LogisticRegression()))
    Models.append(('LDA'  , LinearDiscriminantAnalysis()))
    Models.append(('KNN'  , KNeighborsClassifier()))
    Models.append(('CART' , DecisionTreeClassifier()))
    Models.append(('NB'   , GaussianNB()))
    Models.append(('SVM'  , SVC(probability=True)))
    return Models


# In[ ]:


#appending Algorithms
def ensemblemodels():
    ensembles = []
    ensembles.append(('AB'   , AdaBoostClassifier()))
    ensembles.append(('GBM'  , GradientBoostingClassifier()))
    ensembles.append(('RF'   , RandomForestClassifier()))
    ensembles.append(( 'Bagging' , BaggingClassifier()))
    ensembles.append(('ET', ExtraTreesClassifier()))
    return ensembles


# In[ ]:


# Normalized Models
def NormalizedModel(nameOfScaler):
    
    if nameOfScaler == 'standard':
        scaler = StandardScaler()
    elif nameOfScaler =='minmax':
        scaler = MinMaxScaler()
    elif nameOfScaler == 'normalizer':
        scaler = Normalizer()
    elif nameOfScaler == 'binarizer':
        scaler = Binarizer()

    pipelines = []
    pipelines.append((nameOfScaler+'LR'  , Pipeline([('Scaler', scaler),('LR'  , LogisticRegression())])))
    pipelines.append((nameOfScaler+'LDA' , Pipeline([('Scaler', scaler),('LDA' , LinearDiscriminantAnalysis())])))
    pipelines.append((nameOfScaler+'KNN' , Pipeline([('Scaler', scaler),('KNN' , KNeighborsClassifier())])))
    pipelines.append((nameOfScaler+'CART', Pipeline([('Scaler', scaler),('CART', DecisionTreeClassifier())])))
    pipelines.append((nameOfScaler+'NB'  , Pipeline([('Scaler', scaler),('NB'  , GaussianNB())])))
    pipelines.append((nameOfScaler+'SVM' , Pipeline([('Scaler', scaler),('SVM' , SVC())])))
    pipelines.append((nameOfScaler+'AB'  , Pipeline([('Scaler', scaler),('AB'  , AdaBoostClassifier())])  ))
    pipelines.append((nameOfScaler+'GBM' , Pipeline([('Scaler', scaler),('GMB' , GradientBoostingClassifier())])  ))
    pipelines.append((nameOfScaler+'RF'  , Pipeline([('Scaler', scaler),('RF'  , RandomForestClassifier())])  ))
    pipelines.append((nameOfScaler+'ET'  , Pipeline([('Scaler', scaler),('ET'  , ExtraTreesClassifier())])  ))

    return pipelines


# In[ ]:


# Train model
def fit_model(X_train, y_train,models):
    # Test options and evaluation metric
    num_folds = 10
    scoring = 'accuracy'

    results = []
    names = []
    for name, model in models:
        kfold = KFold(n_splits=num_folds, shuffle=True, random_state=0)
        cv_results = cross_val_score(model, X_train, y_train, cv=kfold, scoring=scoring)
        results.append(cv_results)
        names.append(name)
        msg = "%s: %f (%f)" % (name, cv_results.mean(), cv_results.std())
        print(msg)
        
    return names, results


# In[ ]:


# Save trained model
def save_model(model,filename):
    pickle.dump(model, open(filename, 'wb'))


# In[ ]:


# Performance Measure
def classification_metrics(model, conf_matrix):
    print(f"Training Accuracy Score: {model.score(X_train, y_train) * 100:.1f}%")
    print(f"Validation Accuracy Score: {model.score(X_test, y_test) * 100:.1f}%")
    fig,ax = plt.subplots(figsize=(8,6))
    sns.heatmap(pd.DataFrame(conf_matrix), annot = True, cmap = 'YlGnBu',fmt = 'g')
    ax.xaxis.set_label_position('top')
    plt.tight_layout()
    plt.title('Confusion Matrix', fontsize=20, y=1.1)
    plt.ylabel('Actual label', fontsize=15)
    plt.xlabel('Predicted label', fontsize=15)
    plt.show()
    print(classification_report(y_test, y_pred))


# In[ ]:


# ROC_AUC
def roc_auc(y_test, y_pred):
    fpr, tpr, thresholds = roc_curve(y_test, y_pred)
    plt.figure(figsize=(8,6))
    print(f"roc_auc score: {auc(fpr, tpr)*100:.1f}%")
    plt.plot(fpr, tpr, color='orange', label='ROC')
    plt.plot([0, 1], [0, 1], color='darkblue', linestyle='--')
    plt.xlabel('False Positive Rate',fontsize=12)
    plt.ylabel('True Positive Rate', fontsize=12)
    plt.title('Receiver Operating Characteristic (ROC) Curve', fontsize=20)
    plt.legend()
    plt.show()


# In[ ]:


#DataAnalysis
explore_data(df)


# In[ ]:


#Removing_Duplicates
checking_removing_duplicates(df)


# In[ ]:


df.isna().sum()


# In[ ]:


#Removing Outliers outliers
Q1 = df.quantile(0.25)
Q3 = df.quantile(0.75)
IQR = Q3 - Q1

df_out = df[~((df < (Q1 - 1.5 * IQR)) |(df > (Q3 + 1.5 * IQR))).any(axis=1)]


# In[ ]:


#Training the model
target ='Item'
X_train, X_test, y_train, y_test = read_in_and_split_data(df, target)

models = GetModel()
names,results = fit_model(X_train, y_train,models)


# In[ ]:


#Preprocessing Technique
ScaledModel = NormalizedModel('minmax')
name,results = fit_model(X_train, y_train, ScaledModel)


# In[ ]:


#Saving the model
pipeline = make_pipeline(MinMaxScaler(),  GaussianNB())
model = pipeline.fit(X_train, y_train)
y_pred = model.predict(X_test)
conf_matrix = confusion_matrix(y_test,y_pred)
classification_metrics(pipeline, conf_matrix)
# save model
save_model(model, 'model.pkl')


# In[ ]: