# -*- coding: utf-8 -*-
"""Supplementary Material_1.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1AO_vTTqq0Dc-CKU3jlp4MJ1mfxPPwiaL

Impacts of COVID-19 on Public Universities in Brazil: A Machine Learning Counterfactual Analysis

R. Rossi Jr.
romeu.rossi@ufv.br
"""

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.metrics import mean_absolute_error

"""##Importing the data

Importing data from the INEP annual Higher Education Census
"""

df_12 = pd.read_csv('/content/drive/MyDrive/Dados UFs/Censo/MICRODADOS_CADASTRO_CURSOS_2012.CSV', encoding='ISO-8859-1', sep = ';')
df_13 = pd.read_csv('/content/drive/MyDrive/Dados UFs/Censo/MICRODADOS_CADASTRO_CURSOS_2013.CSV', encoding='ISO-8859-1', sep = ';')
df_14 = pd.read_csv('/content/drive/MyDrive/Dados UFs/Censo/MICRODADOS_CADASTRO_CURSOS_2014.CSV', encoding='ISO-8859-1', sep = ';')
df_15 = pd.read_csv('/content/drive/MyDrive/Dados UFs/Censo/MICRODADOS_CADASTRO_CURSOS_2015.CSV', encoding='ISO-8859-1', sep = ';')
df_16 = pd.read_csv('/content/drive/MyDrive/Dados UFs/Censo/MICRODADOS_CADASTRO_CURSOS_2016.CSV', encoding='ISO-8859-1', sep = ';')
df_17 = pd.read_csv('/content/drive/MyDrive/Dados UFs/Censo/MICRODADOS_CADASTRO_CURSOS_2017.CSV', encoding='ISO-8859-1', sep = ';')
df_18 = pd.read_csv('/content/drive/MyDrive/Dados UFs/Censo/MICRODADOS_CADASTRO_CURSOS_2018.CSV', encoding='ISO-8859-1', sep = ';')
df_19 = pd.read_csv('/content/drive/MyDrive/Dados UFs/Censo/MICRODADOS_CADASTRO_CURSOS_2019.CSV', encoding='ISO-8859-1', sep = ';')
df_20 = pd.read_csv('/content/drive/MyDrive/Dados UFs/Censo/MICRODADOS_CADASTRO_CURSOS_2020.CSV', encoding='ISO-8859-1', sep = ';')
df_21 = pd.read_csv('/content/drive/MyDrive/Dados UFs/Censo/MICRODADOS_CADASTRO_CURSOS_2021.CSV', encoding='ISO-8859-1', sep = ';')

lista_colunas = ['NU_ANO_CENSO','NO_REGIAO', 'CO_REGIAO', 'NO_UF', 'SG_UF', 'CO_UF', 'NO_MUNICIPIO',
                            'CO_MUNICIPIO', 'IN_CAPITAL', 'CO_IES', 'NO_CURSO', 'CO_CURSO', 'CO_CINE_AREA_GERAL', 'NO_CINE_AREA_GERAL',
                            'TP_GRAU_ACADEMICO', 'TP_DIMENSAO', 'TP_MODALIDADE_ENSINO', 'TP_NIVEL_ACADEMICO',
             'TP_CATEGORIA_ADMINISTRATIVA', 'TP_REDE', 'IN_GRATUITO',
             'QT_VG_NOVA', 'QT_VG_TOTAL','QT_INSC_VG_NOVA', 'QT_MAT','QT_CONC','QT_SIT_TRANCADA','QT_SIT_DESVINCULADO', 'QT_SIT_TRANSFERIDO', 'QT_ING']
lista = [df_12, df_13, df_14,df_15, df_16, df_17, df_18, df_19, df_20, df_21]
for n in lista:
  n.drop(columns=n.columns[~n.columns.isin(lista_colunas)], inplace=True)

#Selecting TP_DIMENSAO = 1 (Face-to-face course) TP_MODALIDADE_ENSINO = 1 (Face-to-face) TP_GRAU_ACADEMICO = (1. Bachelor's degree
#2. Teaching degree 4. Bachelor's and Teaching degree) TP_NIVEL_ACADEMICO = (Undergraduate) TP_CATEGORIA_ADMINISTRATIVA = (1. Public
#2. State, 3. Municipal, 4. For-profit private, 5. Non-profit private
#6. Private - Strictly private) TP_REDE = 1 (Public)
lista = [df_12, df_13, df_14,df_15, df_16, df_17, df_18, df_19, df_20, df_21]
for n in lista:
  #n.query("TP_DIMENSAO == 1", inplace=True)
  n.query("TP_GRAU_ACADEMICO == 1 | TP_GRAU_ACADEMICO == 2 | TP_GRAU_ACADEMICO == 4", inplace=True)
  n.query("TP_MODALIDADE_ENSINO == 1", inplace=True)
  n.query("TP_NIVEL_ACADEMICO == 1", inplace=True)
  n.query("TP_CATEGORIA_ADMINISTRATIVA == 1 | TP_CATEGORIA_ADMINISTRATIVA == 2", inplace=True)
  #n.query("TP_CATEGORIA_ADMINISTRATIVA == 4 | TP_CATEGORIA_ADMINISTRATIVA == 5| TP_CATEGORIA_ADMINISTRATIVA == 6", inplace=True)
  #n.query("TP_REDE == 1", inplace=True) a df_13 tem TP_REDE nula
  n.query("IN_GRATUITO == 1", inplace=True)

#Deleting columns with a single value
lista = [df_12, df_13, df_14,df_15, df_16, df_17, df_18, df_19, df_20, df_21]
for n in lista:
  n.drop('TP_DIMENSAO', axis = 1, inplace =True)
  n.drop('TP_MODALIDADE_ENSINO', axis = 1, inplace =True)
  n.drop('TP_NIVEL_ACADEMICO', axis = 1, inplace =True)
  n.drop('TP_REDE', axis = 1, inplace =True)
  #n.drop('IN_GRATUITO', axis = 1, inplace =True)

"""###Completing null values of QT_VG_NOVA in df_12 and df_13"""

lista = [df_12, df_13, df_14,df_15, df_16, df_17, df_18, df_19, df_20, df_21]
for n in lista:
  print('Média QT_VG_NOVA -', n['QT_VG_NOVA'].mean())
  print('Média QT_VG_TOTAL -', n['QT_VG_TOTAL'].mean())
  print()
#df.groupby('NU_ANO_CENSO')['QT_VG_NOVA'].mean()

df_auxiliar = df_14[['CO_CURSO', 'QT_VG_NOVA']]

df_12 = df_12.merge(df_auxiliar, how='inner', on='CO_CURSO')
df_12.head()

df_12 = df_12.drop(columns = 'QT_VG_NOVA_x')

df_12 = df_12.rename(columns={'QT_VG_NOVA_y': 'QT_VG_NOVA'})

df_12.head()

df_13 = df_13.merge(df_auxiliar, how='inner', on='CO_CURSO')
df_13 = df_13.drop(columns = 'QT_VG_NOVA_x')
df_13 = df_13.rename(columns={'QT_VG_NOVA_y': 'QT_VG_NOVA'})

df_13.head()

"""##Concatenating df"""

#Concatenando dfs pos e pre pandemia
df = pd.concat([df_12, df_13, df_14,df_15, df_16, df_17, df_18, df_19,df_20, df_21])
# df_pre = pd.concat([df_12, df_13, df_14,df_15, df_16, df_17, df_18, df_19])
# df_pos = pd.concat([df_20, df_21])

"""##Data Cleaning"""

df.isna().sum().sort_values(ascending = False)

print(df.duplicated().sum())

#Eliminating rows with QT_ING = 0 and QT_MAT=0
df = df[df['QT_ING']>0]
df = df[df['QT_MAT']>0]

df.shape

"""##Lags"""

df.columns

#For each Course (CO_COURSE) creates lags_1 and lag_2 of
  # ['QT_VG_NOVA', 'QT_INSC_VG_NOVA', 'QT_ING', 'QT_MAT', 'QT_CONC','QT_SIT_LOCKED', 'QT_SIT_DESINNCULADO', 'QT_SIT_TRANSFERIDO']
lista_lags = ['QT_VG_NOVA', 'QT_VG_TOTAL','QT_INSC_VG_NOVA', 'QT_ING', 'QT_MAT', 'QT_CONC','QT_SIT_TRANCADA', 'QT_SIT_DESVINCULADO', 'QT_SIT_TRANSFERIDO']
#colunas_lags1 = ['lag1_QT_VG_NOVA', 'lag1_QT_VG_TOTAL','lag1_QT_INSC_VG_NOVA', 'lag1_QT_ING', 'lag1_QT_MAT', 'lag1_QT_CONC','lag1_QT_SIT_TRANCADA', 'lag1_QT_SIT_DESVINCULADO', 'lag1_QT_SIT_TRANSFERIDO']
colunas_lags2 = ['lag2_QT_VG_NOVA', 'lag2_QT_VG_TOTAL','lag2_QT_INSC_VG_NOVA', 'lag2_QT_ING', 'lag2_QT_MAT', 'lag2_QT_CONC','lag2_QT_SIT_TRANCADA', 'lag2_QT_SIT_DESVINCULADO', 'lag2_QT_SIT_TRANSFERIDO']
colunas_lags3 = ['lag3_QT_VG_NOVA', 'lag3_QT_VG_TOTAL','lag3_QT_INSC_VG_NOVA', 'lag3_QT_ING', 'lag3_QT_MAT', 'lag3_QT_CONC','lag3_QT_SIT_TRANCADA', 'lag3_QT_SIT_DESVINCULADO', 'lag3_QT_SIT_TRANSFERIDO']
#df[colunas_lags1] = df.groupby(['CO_CURSO'])[lista_lags].shift(1)
df[colunas_lags2] = df.groupby(['CO_CURSO'])[lista_lags].shift(2)
df[colunas_lags3] = df.groupby(['CO_CURSO'])[lista_lags].shift(3)
# df_pos[colunas_lags] = df_pos.groupby(['CO_CURSO'])[lista_lags].shift(1)
# df_pre[colunas_lags] = df_pre.groupby(['NU_ANO_CENSO', 'CO_CURSO'])[lista_lags].shift(1)
# df_pos[colunas_lags] = df_pos.groupby(['NU_ANO_CENSO', 'CO_CURSO'])[lista_lags].shift(1)

df = df.dropna()

df.shape

"""##Adding population column and geographic regions

####df_pop22 population
"""

df_pop22 = pd.read_excel('/content/drive/MyDrive/Dados UFs/Censo/POP2022_Municipios.xls')
#df_pop22 = pd.read_excel('/content/Pasta1.xlsx')

df_pop22.head(3)

df_pop22.shape

df_pop22.columns

df_pop22.info()

# df_pop22['COD. MUNIC'] = df_pop22['COD. MUNIC'].astype('str')
# df_pop22['COD. UF'] = df_pop22['COD. UF'].astype('str')

# df_pop22['CO_MUNICIPIO'] = df_pop22['COD. UF']+df_pop22['COD. MUNIC']

# df_pop22['CO_MUNICIPIO'] = df_pop22['CO_MUNICIPIO'].astype('int')

df_pop22 = df_pop22.rename(columns={'COD. MUNIC':'Cod_munic'})
df_pop22 = df_pop22.rename(columns={'COD. UF':'CO_UF'})

df_pop22 = df_pop22.rename(columns={'NOME DO MUNICÍPIO':'NO_MUNICIPIO'})

df_pop22.head(3)

#df_pop22 = df_pop22[['Cod_munic','POPULAÇÃO']]

df_pop22.head()

df_pop22.isna().sum()

df_pop22['POPULAÇÃO'] = df_pop22['POPULAÇÃO'].astype('float')

df_pop22[df_pop22['NO_MUNICIPIO'] == 'Itapuã do Oeste']

"""####df_reg geographic regions"""

df_reg = pd.read_excel('/content/drive/MyDrive/Dados UFs/Censo/regioes_geograficas_composicao_por_municipios_2017_20180911.xlsx')
#df_pop22 = pd.read_excel('/content/Pasta1.xlsx')

df_reg.head(3)

df_reg.info()

df_reg = df_reg.rename(columns = {'nome_mun':'NO_MUNICIPIO'})

df_reg = df_reg.rename(columns = {'CD_GEOCODI':'CO_MUNICIPIO'})

df_reg['NO_MUNICIPIO'] = df_reg['NO_MUNICIPIO'].replace('Moji Mirim', 'Mogi Mirim')
df_reg['NO_MUNICIPIO'] = df_reg['NO_MUNICIPIO'].replace('Restinga Seca', 'Restinga Sêca')

df_reg.shape

df_reg.isna().sum()

df_reg.duplicated().sum()

df_reg.head(2)

# Adaptando a coluna de codigo do municipio do df_reg para ficar com os mesmos valores de df_pop22
df_reg['Cod_munic'] = df_reg['CO_MUNICIPIO'].astype('str')
df_reg['Cod_munic'] = df_reg['Cod_munic'].apply(lambda x: x[2:])
df_reg['Cod_munic'] = df_reg['Cod_munic'].astype('int')

df_reg.head(2)

# Criando a coluna CO_UF e e convertendo de inteiro para string
df_reg['CO_UF'] = df_reg['CO_MUNICIPIO'].astype('str')

# Extrair os dois primeiros dígitos da coluna 'CO_UF'
df_reg['CO_UF'] = df_reg['CO_UF'].str[:2]

# # Converter a coluna 'CO_UF' de volta para inteiro
df_reg['CO_UF'] = df_reg['CO_UF'].astype('int')

df_reg.head()

"""####Join df_pop22 and  df_reg"""

df_pop22.shape

df_pop22 = pd.merge(df_pop22, df_reg, on = ['CO_UF','Cod_munic'], how = 'left')

df_pop22.shape

df_pop22.isna().sum()

df_pop22[df_pop22['Cod_munic']==8254]

df_pop22 = df_pop22.rename(columns = {'NO_MUNICIPIO_y':'NO_MUNICIPIO'})

#Ficando apenas com o 'NO_MUNICIPIO' df_reg, pois tem menos erros.
df_pop22 = df_pop22.drop(columns = 'NO_MUNICIPIO_x')

df_pop22.head()

"""####Created columns with pop_rgi and pop_rgint in df_pop22"""

df_pop22['pop_rgi'] = df_pop22.groupby('cod_rgi')['POPULAÇÃO'].transform(sum)

df_pop22['pop_rgint'] = df_pop22.groupby('cod_rgint')['POPULAÇÃO'].transform(sum)

"""####Join df_pop22"""

df['CO_UF'] = df['CO_UF'].astype('int')

df.head(3)

#df = pd.merge(df, df_pop22, on = ['CO_UF','Cod_munic'], how = 'left')
df = pd.merge(df, df_pop22, on = 'CO_MUNICIPIO', how = 'left')

df.shape

df = df.rename(columns = {'NO_MUNICIPIO_x':'NO_MUNICIPIO'})
df = df.rename(columns = {'CO_UF_x':'CO_UF'})

df = df.drop(columns = 'CO_UF_y')

df = df.drop(columns = 'NO_MUNICIPIO_y')

df = df.drop(columns = 'Cod_munic')
df = df.drop(columns = 'UF')

"""####Creating column population/pop_rgi and population/pop_rgint"""

df['pop/pop_rgi'] = df['POPULAÇÃO']/df['pop_rgi']
df['pop/pop_rgint'] = df['POPULAÇÃO']/df['pop_rgint']

df['pop/pop_rgi'].describe(percentiles=[.1,.2,0.33, 0.66, 0.99])

"""###Creating POP_CAT (0, 1, 2)"""

df.POPULAÇÃO.describe(percentiles=[0.33, 0.66, 0.99])

#pop_33_percentile = np.percentile(df['POPULAÇÃO'], 33)
np.percentile(df['POPULAÇÃO'], 33)

# Function that classifies the population
def classificar_populacao(populacao):
    if populacao <= np.percentile(df['POPULAÇÃO'], 33):
        return 0
    elif np.percentile(df['POPULAÇÃO'], 33) <= populacao <= np.percentile(df['POPULAÇÃO'], 66):
        return 1
    else:
        return 2

# The categorical column with the population classification
df['POP_CAT'] = df['POPULAÇÃO'].apply(classificar_populacao)

df[df.POP_CAT==0]['POPULAÇÃO'].max()

df.pop_rgi.describe(percentiles=[0.33, 0.66, 0.99])

# Function for  pop_rgi
def classificar_populacao(pop_rgi):
    if pop_rgi <= np.percentile(df['pop_rgi'], 33):
        return 0
    elif np.percentile(df['pop_rgi'], 33) <= pop_rgi <= np.percentile(df['pop_rgi'], 66):
        return 1
    else:
        return 2

df['pop_rgi_CAT'] = df['pop_rgi'].apply(classificar_populacao)

df[df.pop_rgi_CAT==0]['pop_rgi'].max()

df.pop_rgint.describe(percentiles=[0.33, 0.66, 0.99])

# Function for pop_rgint
def classificar_populacao(pop_rgint):
    if pop_rgint <= np.percentile(df['pop_rgint'], 33):
        return 0
    elif np.percentile(df['pop_rgint'], 33) <= pop_rgint <= np.percentile(df['pop_rgint'], 66):
        return 1
    else:
        return 2

df['pop_rgint_CAT'] = df['pop_rgint'].apply(classificar_populacao)

"""##Outliers"""

df[['QT_VG_NOVA', 'QT_VG_TOTAL','QT_INSC_VG_NOVA', 'QT_ING', 'QT_MAT', 'QT_CONC',
       'QT_SIT_TRANCADA', 'QT_SIT_DESVINCULADO', 'QT_SIT_TRANSFERIDO']].describe(percentiles=[0.01, 0.05,0.1, 0.5, 0.9, 0.95,0.99])

df.shape

#The function returns a Boolean array of the same size as the input
#array, indicating whether or not the value is an outlier:
def is_outlier(array, extreme=False):
  """Custom function to identify outliers in an array"""
  q1,q3 = np.quantile(array, [.25,.75])
  iqr = q3-q1

  factor = 3. if extreme else 1.5
  upper_outlier = q3 + factor*iqr
  lower_outlier = q1 - factor*iqr

  return (array < lower_outlier) | (array > upper_outlier)

df[['QT_VG_NOVA', 'QT_INSC_VG_NOVA', 'QT_ING', 'QT_MAT', 'QT_CONC',
       'QT_SIT_TRANCADA', 'QT_SIT_DESVINCULADO', 'QT_SIT_TRANSFERIDO']].apply(is_outlier).sum()

#Removing outliers from QT_ING
df = df[~df[['QT_ING']].apply(is_outlier).any(axis=1)]
df = df[~df[['QT_VG_NOVA']].apply(is_outlier).any(axis=1)]

df.shape

"""##Correlations"""

df_corr = df[['QT_ING', 'lag3_QT_VG_NOVA', 'lag3_QT_VG_TOTAL','lag3_QT_INSC_VG_NOVA', 'lag3_QT_ING', 'lag3_QT_MAT',
       'lag3_QT_CONC', 'lag3_QT_SIT_TRANCADA', 'lag3_QT_SIT_DESVINCULADO',
       'lag3_QT_SIT_TRANSFERIDO', 'lag2_QT_VG_NOVA', 'lag2_QT_VG_TOTAL','lag2_QT_INSC_VG_NOVA',
       'lag2_QT_ING', 'lag2_QT_MAT', 'lag2_QT_CONC', 'lag2_QT_SIT_TRANCADA',
       'lag2_QT_SIT_DESVINCULADO', 'lag2_QT_SIT_TRANSFERIDO', 'POPULAÇÃO',
       'pop_rgi', 'pop_rgint', 'pop/pop_rgi', 'pop/pop_rgint']]
# df_corr = df[['QT_ING','QT_VG_NOVA', 'QT_VG_TOTAL','QT_INSC_VG_NOVA', 'QT_MAT', 'QT_CONC',
#        'QT_SIT_TRANCADA', 'QT_SIT_DESVINCULADO', 'QT_SIT_TRANSFERIDO',
#        'lag1_QT_VG_NOVA', 'lag1_QT_VG_TOTAL','lag1_QT_INSC_VG_NOVA', 'lag1_QT_ING', 'lag1_QT_MAT',
#        'lag1_QT_CONC', 'lag1_QT_SIT_TRANCADA', 'lag1_QT_SIT_DESVINCULADO',
#        'lag1_QT_SIT_TRANSFERIDO', 'lag2_QT_VG_NOVA', 'lag2_QT_VG_TOTAL','lag2_QT_INSC_VG_NOVA',
#        'lag2_QT_ING', 'lag2_QT_MAT', 'lag2_QT_CONC', 'lag2_QT_SIT_TRANCADA',
#        'lag2_QT_SIT_DESVINCULADO', 'lag2_QT_SIT_TRANSFERIDO', 'POPULAÇÃO',
#        'pop_rgi', 'pop_rgint', 'pop/pop_rgi', 'pop/pop_rgint']]

len(df_corr.columns)

# Calculates correlations
correlations = df_corr.corr()['QT_ING']
correlations.sort_values(ascending = False)

# Filters columns with correlation greater than 0.15
correlated_columns = correlations[correlations > 0.15].index.tolist()
correlated_columns

"""##Modeling"""

from xgboost import XGBRegressor
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error
from sklearn.metrics import mean_absolute_percentage_error
from sklearn.svm import SVR
from sklearn.tree import DecisionTreeRegressor
from sklearn.linear_model import LinearRegression, Ridge, Lasso, LassoCV, ElasticNet, ElasticNetCV, HuberRegressor, LassoLars, BayesianRidge
from sklearn.linear_model import SGDRegressor
from sklearn.ensemble import RandomForestRegressor, AdaBoostRegressor, GradientBoostingRegressor, ExtraTreesRegressor
from sklearn import metrics
from math import sqrt

#defining df_pos and df_pre pandemic
df_pre = df[df.NU_ANO_CENSO < 2020]
df_pos = df[df.NU_ANO_CENSO >= 2020]

"""###Selecting the best model"""

correlated_columns.remove('QT_ING')
colunas_prev = correlated_columns
# colunas_prev = ['NU_ANO_CENSO', 'CO_REGIAO', 'CO_UF','CO_MUNICIPIO', 'IN_CAPITAL',
#        'TP_CATEGORIA_ADMINISTRATIVA', 'CO_IES', 'CO_CURSO',
#        'CO_CINE_AREA_GERAL', 'TP_GRAU_ACADEMICO',
#        'QT_VG_NOVA', 'QT_INSC_VG_NOVA', 'QT_MAT', 'QT_CONC',
#        'QT_SIT_TRANCADA', 'QT_SIT_DESVINCULADO', 'QT_SIT_TRANSFERIDO',
#        'lag1_QT_VG_NOVA', 'lag1_QT_INSC_VG_NOVA', 'lag1_QT_ING', 'lag1_QT_MAT',
#        'lag1_QT_CONC', 'lag1_QT_SIT_TRANCADA', 'lag1_QT_SIT_DESVINCULADO',
#        'lag1_QT_SIT_TRANSFERIDO', 'lag2_QT_VG_NOVA', 'lag2_QT_INSC_VG_NOVA',
#        'lag2_QT_ING', 'lag2_QT_MAT', 'lag2_QT_CONC', 'lag2_QT_SIT_TRANCADA',
#        'lag2_QT_SIT_DESVINCULADO', 'lag2_QT_SIT_TRANSFERIDO','POP_CAT', 'pop_rgi_CAT',
#        'pop_rgint_CAT']

colunas_prev = correlated_columns

colunas_prev

df_pre.NU_ANO_CENSO.value_counts()

# Separating data into training and testing sets for each cross-validation step
df_train_15 = df_pre[df_pre['NU_ANO_CENSO'].isin([2014, 2015])]
df_train_16 = df_pre[df_pre['NU_ANO_CENSO'].isin([2014, 2015,2016])]
df_train_17 = df_pre[df_pre['NU_ANO_CENSO'].isin([2014, 2015,2016,2017])]

# df_train_2016 = df_pre[df_pre['NU_ANO_CENSO'] == 2016]
# df_train_2017 = df_pre[df_pre['NU_ANO_CENSO'].isin([2016, 2017])]
# df_train_2018 = df_pre[df_pre['NU_ANO_CENSO'].isin([2016, 2017, 2018])]

df_test_16_17 = df_pre[df_pre['NU_ANO_CENSO'].isin([2016, 2017])]
df_test_17_18 = df_pre[df_pre['NU_ANO_CENSO'].isin([2017, 2018])]
df_test_18_19 = df_pre[df_pre['NU_ANO_CENSO'].isin([2018, 2019])]


# df_test_2017 = df_pre[df_pre['NU_ANO_CENSO'] == 2017]
# df_test_2018 = df_pre[df_pre['NU_ANO_CENSO'] == 2018]
# df_test_2019 = df_pre[df_pre['NU_ANO_CENSO'] == 2019]

# Defining the features and target
features = colunas_prev
target = 'QT_ING'

#Initializing the list to store the MAEs
mae_dict = {}


# Models
models = {
    'Linear Regression': LinearRegression(),
    'XGBoost': XGBRegressor(),
    'Random Forest': RandomForestRegressor(),
    'Extra Trees': ExtraTreesRegressor()
}

# Loop through models
for model_name, model in models.items():
    mae_list = []

    # First stage: data from 2016 and 2017
    X_train = df_train_15[features]
    y_train = df_train_15[target]
    X_test = df_test_16_17[features]
    y_test = df_test_16_17[target]

    # Training the model
    model.fit(X_train, y_train)

    # Predicting test set values
    y_pred = model.predict(X_test)

    # Calculating the MAE
    mae = mean_absolute_error(y_test, y_pred)
    mae_list.append(mae)

    # Second stage: data from 2016, 2017 and 2018
    X_train = df_train_16[features]
    y_train = df_train_16[target]
    X_test = df_test_17_18[features]
    y_test = df_test_17_18[target]


    model.fit(X_train, y_train)

    y_pred = model.predict(X_test)

    mae = mean_absolute_error(y_test, y_pred)
    mae_list.append(mae)

    # Third stage: data from 2016, 2017, 2018 and 2019
    X_train = df_train_17[features]
    y_train = df_train_17[target]
    X_test = df_test_18_19[features]
    y_test = df_test_18_19[target]

    model.fit(X_train, y_train)

    y_pred = model.predict(X_test)

    mae = mean_absolute_error(y_test, y_pred)
    mae_list.append(mae)

    # Storing the list of MAEs for the current model
    mae_dict[model_name] = mae_list


# Printing the results
for model_name, maes in mae_dict.items():
    print("Model:", model_name)
    print("MAEs:", maes)
    print("Average MAE:", sum(maes) / len(maes))
    print()

df_test_18_19[target].describe()

"""####Hyperparam tunning Extra Trees"""

from sklearn.model_selection import ParameterGrid

#ref https://www.kaggle.com/code/manovirat/timeseries-using-prophet-hyperparameter-tuning
params_grid = {
    'n_estimators': [250 ,500, 1000],
    'max_features':[None,1]
    # ,'max_depth': [None, 15]
    # ,'min_samples_split': [2, 10],
    #  'min_samples_leaf': [1, 10]
    }
grid = ParameterGrid(params_grid)
cnt = 0
for p in grid:
    cnt = cnt+1

print('Total Possible Models',cnt)

# Separating data into training and testing sets for each cross-validation step
df_train_15 = df_pre[df_pre['NU_ANO_CENSO'].isin([2014, 2015])]
df_train_16 = df_pre[df_pre['NU_ANO_CENSO'].isin([2014, 2015,2016])]
df_train_17 = df_pre[df_pre['NU_ANO_CENSO'].isin([2014, 2015,2016,2017])]


df_test_16_17 = df_pre[df_pre['NU_ANO_CENSO'].isin([2016, 2017])]
df_test_17_18 = df_pre[df_pre['NU_ANO_CENSO'].isin([2017, 2018])]
df_test_18_19 = df_pre[df_pre['NU_ANO_CENSO'].isin([2018, 2019])]



# Features e o target
features = colunas_prev
target = 'QT_ING'

# Initializing the list to store the MAEs
mae_dict = {}

model_parameters = pd.DataFrame(columns = ['Average_MAE','Parameters'])

for p in grid:
    mae_list = []

    # First stage: data from 2016 and 2017
    X_train = df_train_15[features]
    y_train = df_train_15[target]
    X_test = df_test_16_17[features]
    y_test = df_test_16_17[target]

    ext_model = ExtraTreesRegressor(n_estimators = p['n_estimators'],
     max_features = p['max_features']
    # ,max_depth = p['max_depth']
    #  ,min_samples_split = p['min_samples_split'],
    #  min_samples_leaf = p['min_samples_leaf']
                                     )

    # Training the model
    ext_model.fit(X_train, y_train)

    # Predicting test set values
    y_pred = ext_model.predict(X_test)

    # Calculating the MAE
    mae = mean_absolute_error(y_test, y_pred)
    mae_list.append(mae)

    # Second stage: data from 2016, 2017 and 2018
    X_train = df_train_16[features]
    y_train = df_train_16[target]
    X_test = df_test_17_18[features]
    y_test = df_test_17_18[target]

    # Training the model again
    ext_model.fit(X_train, y_train)

    y_pred = ext_model.predict(X_test)

    mae = mean_absolute_error(y_test, y_pred)
    mae_list.append(mae)

    # Third stage: data from 2016, 2017, 2018 and 2019
    X_train = df_train_17[features]
    y_train = df_train_17[target]
    X_test = df_test_18_19[features]
    y_test = df_test_18_19[target]

    ext_model.fit(X_train, y_train)

    y_pred = ext_model.predict(X_test)

    mae = mean_absolute_error(y_test, y_pred)
    mae_list.append(mae)

    Average_MAE = sum(mae_list) / len(mae_list)
    model_parameters = model_parameters.append({'Average_MAE':Average_MAE,'Parameters':p},ignore_index=True)

parameters = model_parameters.sort_values(by=['Average_MAE'])
parameters = parameters.reset_index(drop=True)
parameters

#Melhores parametros
print(parameters.Average_MAE[0])
parameters.Parameters[0]

best_params_xg = parameters.Parameters[0]
best_params_xg

MAE = parameters.Average_MAE[0]
MAE

"""##Counterfactual scenario - Using the best model"""

# Using pre-pandemic data such as df_treino to create the counterfactual scenario (in which the pandemic did not exist),
# using 2020 and 2021 data as df_teste
X_train = df_pre[features]
y_train = df_pre[target]
X_test = df_pos[features]
y_test = df_pos[target]
#Best parameters


best_model_xg = ExtraTreesRegressor(n_estimators =best_params_xg['n_estimators'],
      max_features = best_params_xg['max_features'])


best_model_xg.fit(X_train, y_train)

y_pred = best_model_xg.predict(X_test)

df_contraf = df_pos.copy()

df_contraf['QT_ING_Contraf'] = y_pred

df_contraf.head()

"""The value of QT_ING_Contraf must be at most equal to the value of QT_VG_TOTAL"""

# Replaces values where df_contraf['QT_ING_Contraf'] > df_contraf['QT_VG_TOTAL']
df_contraf.loc[df_contraf['QT_ING_Contraf'] > df_contraf['QT_VG_TOTAL'], 'QT_ING_Contraf'] = df_contraf['QT_VG_TOTAL']

# Replace negative values of QT_ING_Contraf
df_contraf.loc[df_contraf['QT_ING_Contraf'] < 0, 'QT_ING_Contraf'] = 0

"""###Adding the impact_pandemic column"""

#Defining the impact with the difference between the variation in the real scenario and the counterfactual scenario

#Separating df_contraf_20 and df_contraf_21
df_contraf_20 = df_contraf[df_contraf['NU_ANO_CENSO'] == 2020]
df_contraf_21 = df_contraf[df_contraf['NU_ANO_CENSO'] == 2021]

# Filter the df_pre rows for the year 2019 and select only the necessary columns
df_pre_2019 = df_pre[df_pre['NU_ANO_CENSO'] == 2019][['CO_CURSO', 'QT_ING']]

df_contraf_20 = pd.merge(df_contraf_20, df_pre_2019, how='inner', on='CO_CURSO', suffixes=('', '_19'))

#The pandemic_impact in 2020 is defined as::
#((QT_ING_20 - QT_ING_19)/QT_ING_19) - ((QT_ING_Contraf - QT_ING_19)/QT_ING_19) = (QT_ING_20 - QT_ING_Contraf)/QT_ING_19
df_contraf_20 = df_contraf_20[df_contraf_20['QT_ING_19']>1]
df_contraf_20['Impacto_pand'] = (df_contraf_20['QT_ING'] - df_contraf_20['QT_ING_Contraf'])/df_contraf_20['QT_ING_19']

df_contraf_20[['QT_ING','QT_ING_Contraf','QT_ING_19','Impacto_pand']].head(2)

#The Pandemic_Impact in 2021 is defined as (QT_ING_21 - QT_ING_Contraf of 21)/ QT_ING_Contraf of 20
# Filter the df_pre rows for the year 2020 and select only the necessary columns
df_contraf_20_temp = df_contraf_20[['CO_CURSO', 'QT_ING_Contraf']]

df_contraf_21 = pd.merge(df_contraf_21, df_contraf_20_temp, how='inner', on='CO_CURSO', suffixes=('', '_20'))

df_contraf_21 = df_contraf_21[df_contraf_21['QT_ING_Contraf_20']>1]
df_contraf_21['Impacto_pand'] = (df_contraf_21['QT_ING'] - df_contraf_21['QT_ING_Contraf'])/df_contraf_21['QT_ING_Contraf_20']

df_contraf_21[['QT_ING','QT_ING_Contraf','QT_ING_Contraf_20','Impacto_pand']].head(2)

df_contraf = pd.concat([df_contraf_20, df_contraf_21])

df_contraf.isna().sum()

df_contraf = df_contraf.drop(['QT_ING_19', 'QT_ING_Contraf_20'], axis=1)

df_contraf = df_contraf.dropna()

df_contraf.isna().sum()

"""####Selecting courses with a greater impact than the MAE error"""

df_contraf['Impacto_pand_dif'] = df_contraf['QT_ING'] - df_contraf['QT_ING_Contraf']

#MAE of target 'QT_ING' - calculated in the parameter optimization section
MAE

#Selecting only courses with Impacto_pand greater than MAE
#since Impacto_pand = QT_ING' - QT_ING_Contraf
#as QT_ING_Contraf has an error = MAE - the Impacto_pand error is also +- MAE
df_contraf_filtrado = df_contraf[(df_contraf['Impacto_pand_dif']>MAE)|(df_contraf['Impacto_pand_dif']<(-1*MAE))]
df_contraf_filtrado

df_contraf_filtrado[['QT_ING','QT_ING_Contraf','Impacto_pand']].describe(percentiles=[0.01, 0.05,0.1, 0.5, 0.9, 0.95,0.99])

"""####Ouliers"""

#Eliminando os outliers
df_contraf_filtrado = df_contraf_filtrado[(df_contraf_filtrado['Impacto_pand']>df_contraf_filtrado.Impacto_pand.quantile(0.01))
& (df_contraf_filtrado['Impacto_pand']<df_contraf_filtrado.Impacto_pand.quantile(0.99))]

df_contraf_filtrado[['QT_ING','QT_ING_Contraf','Impacto_pand']].describe(percentiles=[0.01, 0.05,0.1, 0.5, 0.9, 0.95,0.99])

def is_outlier(array, extreme=False):
  """Custom function to identify outliers in an array"""
  q1,q3 = np.quantile(array, [.25,.75])
  iqr = q3-q1

  factor = 3. if extreme else 1.5
  upper_outlier = q3 + factor*iqr
  lower_outlier = q1 - factor*iqr

  return (array < lower_outlier) | (array > upper_outlier)

# Calculating number of outliers
df_contraf_filtrado[['QT_ING','QT_ING_Contraf','Impacto_pand']].apply(is_outlier).sum()

df_contraf_filtrado.shape

#Removing outliers from Impacto_pand
df_contraf_filtrado = df_contraf_filtrado[~df_contraf_filtrado[['Impacto_pand']].apply(is_outlier).any(axis=1)]

df_contraf_filtrado.shape

"""####Adding average N of enrolled"""

df_media_mat = df_pre.groupby('CO_CURSO')['QT_MAT'].mean().to_frame().reset_index()
df_media_mat.rename(columns={'QT_MAT': 'QT_MAT_Media'}, inplace=True)
df_media_mat

df_media_mat.describe()

df_media_mat = df_media_mat[df_media_mat.QT_MAT_Media < 1000]

df_contraf_filtrado = df_contraf_filtrado.merge(df_media_mat, on='CO_CURSO', how='left')

df_contraf_filtrado[["NU_ANO_CENSO",'CO_CURSO','QT_MAT_Media']]

df_contraf_filtrado = df_contraf_filtrado.dropna()

df_contraf_filtrado.shape

df_contraf_filtrado.to_csv('df_contraf_filtrado_publica.csv', index=False)

"""###df_dados filtered"""

df_contraf_filtrado = pd.read_csv('/content/drive/MyDrive/Dados UFs/Censo/df_contraf_filtrado_publica.csv')

df_contraf_filtrado.shape

df_contraf_filtrado.head()

"""####Courses and areas

#####Courses
"""

df_contraf_filtrado['NO_CURSO'] = df_contraf_filtrado['NO_CURSO'].str.upper()

df_contraf_filtrado.NO_CURSO.value_counts().describe()

df_contraf_filtrado.NO_CURSO.value_counts().describe()

top_cursos = df_contraf_filtrado.NO_CURSO.value_counts()
total = len(top_cursos)
top_20_percent = top_cursos.head(int(0.1 * total))

top_cursos = top_20_percent
top_cursos

df_contraf_filtrado_top = df_contraf_filtrado[df_contraf_filtrado['NO_CURSO'].isin(top_cursos.index)]

cursos_imp = df_contraf_filtrado_top.groupby('NO_CURSO')['Impacto_pand'].median().sort_values()

cursos_imp.describe()

cursos_imp

cursos_imp.describe()

cursos_imp2 = cursos_imp.to_frame()

cursos_imp2 =cursos_imp2.reset_index()

cursos_imp2

plt.figure(figsize=(10, 6))
plt.barh(cursos_imp2['NO_CURSO'], cursos_imp2['Impacto_pand'], color='blue', alpha=0.7, height=0.9)

# Personalizar o gráfico
plt.xlabel('Impact Factor')
# plt.title('Variação da Oferta Média 20-21 por Curso')
# plt.legend()
# Mostrar o gráfico
plt.tight_layout()
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

import pandas as pd

# Dados fornecidos
data = {
    'NO_CURSO': [
        'MUSIC', 'PHYSICS', 'FOOD ENGINEERING', 'CHEMISTRY', 'MATHEMATICS',
        'LETTERS - PORTUGUESE AND ENGLISH', 'FOREST ENGINEERING', 'BIOLOGICAL SCIENCES',
        'PHILOSOPHY', 'ANIMAL SCIENCE', 'GEOGRAPHY', 'LETTERS - PORTUGUESE LANGUAGE',
        'ENVIRONMENTAL ENGINEERING', 'AGRONOMY', 'SOCIAL SCIENCES', 'HISTORY',
        'LETTERS - PORTUGUESE', 'PRODUCTION ENGINEERING', 'PEDAGOGY',
        'ELECTRICAL ENGINEERING', 'PHYSICAL EDUCATION', 'COMPUTER SCIENCE',
        'ADMINISTRATION', 'NURSING', 'CIVIL ENGINEERING', 'CHEMICAL ENGINEERING',
        'ECONOMICS', 'ACCOUNTING', 'MECHANICAL ENGINEERING',
        'NUTRITION', 'INFORMATION SYSTEMS'
    ],
    'Impacto_pand': [
        -0.453878, -0.4004, -0.350464, -0.318246, -0.304606, -0.30124, -0.299484,
        -0.296781, -0.290972, -0.287421, -0.284776, -0.256542, -0.256458, -0.254187,
        -0.25, -0.246566, -0.242693, -0.229102, -0.216088, -0.197374, -0.193282,
        -0.176411, -0.176111, -0.17479, -0.16334, -0.160437, -0.157684, -0.151128,
        -0.145286, -0.138846, 0.109144
    ]
}

# Criando o DataFrame
df_ingles = pd.DataFrame(data)
df_ingles

import matplotlib.pyplot as plt

plt.figure(figsize=(10, 6))
plt.barh(df_ingles['NO_CURSO'], df_ingles['Impacto_pand'], color='blue', alpha=0.7, height=0.9)

# Personalizar o gráfico
plt.xlabel('Impact Factor')
# plt.title('Variação da Oferta Média 20-21 por Curso')
# plt.legend()
# Mostrar o gráfico
plt.tight_layout()
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

#Top 20 most impacted - how many courses in each POP_CAT
cursos_imp20 = cursos_imp.head(20)
cursos_imp_top20 = df_contraf_filtrado[df_contraf_filtrado['NO_CURSO'].isin(cursos_imp20.index)]
cursos_imp_top20.POP_CAT.value_counts()

df_contraf_filtrado.POP_CAT.value_counts()

#Percentage of top20 most affected in small cities
771/1370

#Percentage of top20 most affected in large cities
320/943

"""#####Courses most affected in small towns (more widespread distribution?)

Hypothesis: The most impacted courses are more impacted because they are widely offered in small cities. Thus, location would be one of the causes for the inequality of impact between courses.
To verify this hypothesis, we can compare the impact distribution between courses in the global scenario with the impact distribution only in small cities (POP_CAT = 0)
"""

df_contraf_filtrado = pd.read_csv('/content/drive/MyDrive/Dados UFs/Censo/df_contraf_filtrado_publica.csv')

df_contraf_filtrado['NO_CURSO'] = df_contraf_filtrado['NO_CURSO'].str.upper()

df_contraf_filtrado_pp = df_contraf_filtrado[df_contraf_filtrado['POP_CAT']==0]

top_cursos_pp = df_contraf_filtrado_pp.NO_CURSO.value_counts()
total = len(top_cursos_pp)
top_20_percent_pp = top_cursos_pp.head(int(0.1 * total))

top_cursos_pp = top_20_percent_pp
top_cursos_pp

df_contraf_filtrado_top_pp = df_contraf_filtrado_pp[df_contraf_filtrado_pp['NO_CURSO'].isin(top_cursos_pp.index)]

cursos_imp_pp = df_contraf_filtrado_top_pp.groupby('NO_CURSO')['Impacto_pand'].median().sort_values()

cursos_imp_pp.describe()

cursos_imp_pp

#Comparing the distributions of general impact (cursos_imp) and impact in small cities (cursos_imp_pp)
#plt.boxplot(cursos_imp);
sns.boxplot(y=cursos_imp)
#sns.boxplot(x="IN_CAPITAL", y='Impacto_pand', data=cursos_imp);

sns.boxplot(y=cursos_imp_pp)

fig, axes = plt.subplots(1, 2, figsize=(12, 6))

sns.boxplot(y=cursos_imp, ax=axes[0])
axes[0].set_title('No city restrictions')
axes[0].set_ylabel('Impact Factor')

sns.boxplot(y=cursos_imp_pp, ax=axes[1])
axes[1].set_title('Small Cities')
axes[1].set_ylabel('Impact Factor')

#IQR
q1 = np.percentile(cursos_imp, 25)  # Primeiro quartil
q3 = np.percentile(cursos_imp, 75)  # Terceiro quartil

iqr = q3 - q1

print("Primeiro Quartil (Q1):", q1)
print("Terceiro Quartil (Q3):", q3)
print("Amplitude Interquartil (IQR):", iqr)

#IQR

q1 = np.percentile(cursos_imp_pp, 25)  # Primeiro quartil
q3 = np.percentile(cursos_imp_pp, 75)  # Terceiro quartil


iqr = q3 - q1

print("Primeiro Quartil (Q1):", q1)
print("Terceiro Quartil (Q3):", q3)
print("Amplitude Interquartil (IQR):", iqr)

"""The variation in Imp_pand in courses offered in small cities is smaller than the variation observed in the general data. This reinforces the hypothesis that the concentration of certain courses in small cities is what causes the variation in impact between areas. Considering only data from small cities, we removed the location factor and compared only the areas. In this case, the areas have a similar impact, with fewer variations.

The distribution with a lower Interquartile Range (IQR) tends to have less variation in values compared to the distribution with a higher IQR. The IQR is a measure of dispersion that represents the difference between the third quartile (Q3) and the first quartile (Q1) of a distribution. The lower the IQR, the smaller the difference between the values in the 25% and 75% quartiles, indicating that most of the data is concentrated in a smaller range.

So if you have two distributions and one of them has a smaller IQR than the other, this suggests that the distribution with the smaller IQR has a greater concentration of its values in a narrower range, which generally indicates less variation in the data compared to the distribution with a higher IQR.

####Campus Features
"""

df_contraf_filtrado_20 = df_contraf_filtrado[df_contraf_filtrado.NU_ANO_CENSO == 2020]
df_contraf_filtrado_21 = df_contraf_filtrado[df_contraf_filtrado.NU_ANO_CENSO == 2021]

"""#####Location and Population"""

fig, axes = plt.subplots(1, 2, figsize=(12, 6))

# Criar o boxplot para cursos_imp no primeiro subplot
sns.boxplot(x="POP_CAT", y='Impacto_pand', data=df_contraf_filtrado_20, ax=axes[0])
axes[0].set_title('Year 2020')
axes[0].set_xlabel('Population Group')
axes[0].set_ylabel('Impact')

# Criar o boxplot para cursos_imp_pp no segundo subplot
sns.boxplot(x="POP_CAT", y='Impacto_pand', data=df_contraf_filtrado_21, ax=axes[1])
axes[1].set_title('Year 2021')
axes[1].set_xlabel('Population Group')
axes[1].set_ylabel('Impact')