## LUCA Telefonica Team¶

#### Carlos Cuervo, Abdul Yaver¶

16 de junio de 2020

Se realiza este análisis respondiendo a la convocatoria del concurso Data Jam CRC y dar respuesta a los dos retos plnateados. Después de analizar la data y los retos, el equipo se plantea las siguientes preguntas de investigación, las cuales se espera ser respondidads haciendo uso de diferentes técnicas de la ciencia de datos:

- ¿Cómo las variables socioeconómicas y macroeconómicas afectan la adopción digital en las diferentes regiones de Colombia?

- ¿Cómo se relacionan las dimensiones más determinantes de la calidad en el servicio de los operadores de telecomunicaciones con la portabilidad numérica neta?

## Procesamiento de data de accesos web:¶

In [114]:

    #Importación de librerías principales a usar:
    %matplotlib inline
    import matplotlib.pyplot as plt
    import numpy as np
    import pandas as pd
    import os
    import seaborn as sns

In [5]:

    os.listdir(os.path.join(os.getcwd()))

Out[5]:

    ['.ipynb_checkpoints',
     'ACCESOS_INTERNET_FIJO.csv',
     'Data_limpia',
     'Data_tableau',
     'LUCA_Telefonica_notebook.ipynb']

In [8]:

    #Importación de la data de accesos web, debido a su tamaño no es adecuado ser procesada en hojas de cálculo:
    list_archivos = os.listdir(os.path.join(os.getcwd()))
    ruta_accesos= os.path.join(os.getcwd(),'ACCESOS_INTERNET_FIJO.csv')

    #Importar y limpiar la data de accesos web por departamento y municipio:
    data_acces= pd.read_csv(ruta_accesos, sep=(';'),engine='python')
    data_acces1 = data_acces.groupby(['ANO','DESC_DEPARTAMENTO','DESC_MUNICIPIO']).aggregate({'CANTIDAD_ACCESOS': np.sum})
    data_acces2 = data_acces1.reset_index()

In [14]:

    #Se observa la data:
    data_acces2.head(10)

Out[14]:

| ANO | DESC_DEPARTAMENTO | DESC_MUNICIPIO | CANTIDAD_ACCESOS
--- | --- | --- | --- | ---
0 | 2010 | AMAZONAS | EL ENCANTO | 122
1 | 2010 | AMAZONAS | LA CHORRERA | 14
2 | 2010 | AMAZONAS | LA PEDRERA | 13
3 | 2010 | AMAZONAS | LETICIA | 1420
4 | 2010 | AMAZONAS | MIRITI - PARANA | 3
5 | 2010 | AMAZONAS | PUERTO ALEGRIA | 3
6 | 2010 | AMAZONAS | PUERTO ARICA | 10
7 | 2010 | AMAZONAS | PUERTO NARIÑO | 58
8 | 2010 | AMAZONAS | PUERTO SANTANDER | 22
9 | 2010 | AMAZONAS | TARAPACA | 39

In [16]:

    data_acces2.describe()

Out[16]:

| ANO | CANTIDAD_ACCESOS
--- | --- | ---
count | 11347.000000 | 1.134700e+04
mean | 2014.481978 | 1.741844e+04
std | 2.844074 | 1.931775e+05
min | 2010.000000 | 0.000000e+00
25% | 2012.000000 | 3.600000e+01
50% | 2014.000000 | 2.920000e+02
75% | 2017.000000 | 2.060000e+03
max | 2019.000000 | 7.596896e+06

In [17]:

    data_acces2.info()

    <class 'pandas.core.frame.DataFrame'>
    RangeIndex: 11347 entries, 0 to 11346
    Data columns (total 4 columns):
    ANO                  11347 non-null int64
    DESC_DEPARTAMENTO    11347 non-null object
    DESC_MUNICIPIO       11347 non-null object
    CANTIDAD_ACCESOS     11347 non-null int64
    dtypes: int64(2), object(2)
    memory usage: 354.7+ KB

In [18]:

    data_acces2.DESC_DEPARTAMENTO.unique()

Out[18]:

    array(['AMAZONAS', 'ANTIOQUIA', 'ARAUCA', 'ATLANTICO', 'BOLIVAR',
           'BOYACA', 'CALDAS', 'CAQUETA', 'CASANARE', 'CAUCA', 'CESAR',
           'CHOCO', 'CORDOBA', 'CUNDINAMARCA', 'GUAINIA', 'GUAVIARE', 'HUILA',
           'LA GUAJIRA', 'MAGDALENA', 'META', 'NARIÑO', 'NORTE DE SANTANDER',
           'PUTUMAYO', 'QUINDIO', 'RISARALDA', 'SAN ANDRES',
           'SANTAFE DE BOGOTA', 'SANTANDER', 'SUCRE', 'TOLIMA',
           'VALLE DEL CAUCA', 'VAUPES', 'VICHADA',
           'ARCHIPIÉLAGO DE SAN ANDRÉS, PROVIDENCIA Y SANTA CATALINA',
           'ATLÁNTICO', 'BOGOTÁ D.C.', 'BOLÍVAR', 'BOYACÁ', 'CHOCÓ',
           'CÓRDOBA', 'GUAINÍA', 'QUINDÍO', 'VAUPÉS', 'COLOMBIA'],
          dtype=object)

Se observa que se debe arreglar la data, ya que a lo largo del tiempo el nombre de algunos departamentos cambió. Es el caso para BOGOTA y para SAN ANDRES. Por otra parte, aparece el departamento "Colombia" el cual caramente no hace referencia a ningún departamento.

In [ ]:

    #Se extrae la información para ser publicada en tableau:
    data_acces2.to_excel(os.path.join(os.getcwd(),'Data_tableau','accesos_web_municipios.xlsx'), index= False)

# Análisis de correlaciones entre variables digitales¶

In [10]:

    #Cargar la data:
    ruta_corr= os.path.join(os.getcwd(),'Data_limpia','Pentracion_digital_correlacion.xlsx')
    data_corr= pd.read_excel(ruta_corr, sheet='Hoja1')

In [12]:

    #Observamos una muestra de los datos:
    data_corr.head(10)

Out[12]:

| departamento | annio | concatenado | poblacion | computador | internet | y_accesos | computador_pobla | internet_pobla | y_accesos_pobla
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
0 | Amazonas | 2018 | Amazonas2018 | 78830 | 1.370000e+04 | 1.460000e+04 | 2516 | 0.173792 | 0.185209 | 0.031917
1 | Antioquia | 2018 | Antioquia2018 | 6691030 | 2.658127e+06 | 4.089047e+06 | 4554537 | 0.397267 | 0.611124 | 0.680693
2 | Arauca | 2018 | Arauca2018 | 270708 | 6.890000e+04 | 1.133000e+05 | 44900 | 0.254518 | 0.418532 | 0.165861
3 | Atlántico | 2018 | Atlántico2018 | 2545924 | 1.003900e+06 | 1.568300e+06 | 1363668 | 0.394317 | 0.616004 | 0.535628
4 | Bogotá D.C. | 2018 | Bogotá D.C.2018 | 8181047 | 4.860400e+06 | 6.111700e+06 | 7367348 | 0.594105 | 0.747056 | 0.900539
5 | Bolívar | 2018 | Bolívar2018 | 2171280 | 8.271000e+05 | 1.123600e+06 | 733946 | 0.380927 | 0.517483 | 0.338025
6 | Boyacá | 2018 | Boyacá2018 | 1282063 | 5.044000e+05 | 6.903000e+05 | 451142 | 0.393428 | 0.538429 | 0.351888
7 | Caldas | 2018 | Caldas2018 | 993866 | 3.487000e+05 | 5.585000e+05 | 487311 | 0.350852 | 0.561947 | 0.490319
8 | Caquetá | 2018 | Caquetá2018 | 496241 | 1.249000e+05 | 1.884000e+05 | 83606 | 0.251692 | 0.379654 | 0.168479
9 | Casanare | 2018 | Casanare2018 | 375249 | 1.465000e+05 | 2.226000e+05 | 136733 | 0.390407 | 0.593206 | 0.364379

In [20]:

    #Coeficiente de correlación entre las variables de computadores, uso de internet y acceso web:
    df_corr= data_corr[['computador_pobla','internet_pobla','y_accesos_pobla']].copy()
    df_corr_nombres= ['computador_pobla','internet_pobla','y_accesos_pobla']
    df_corr.corr(method='pearson')

Out[20]:

| computador_pobla | internet_pobla | y_accesos_pobla
--- | --- | --- | ---
computador_pobla | 1.000000 | 0.924707 | 0.892434
internet_pobla | 0.924707 | 1.000000 | 0.927093
y_accesos_pobla | 0.892434 | 0.927093 | 1.000000

In [21]:

    #Lo vemos matricialmente:
    fig = plt.figure()
    ax = fig.add_subplot(111)
    cax = ax.matshow(df_corr.corr(), interpolation='nearest')
    fig.colorbar(cax)
    ax.set_xticklabels(['']+df_corr_nombres, rotation = 45)
    ax.set_yticklabels(['']+df_corr_nombres)
    plt.show()

In [28]:

    #Se grafica el comportamiento de la relación entre el uso del computador e internet con el acceso web:
    sns.pairplot(data_corr, x_vars=['computador_pobla','internet_pobla'],y_vars='y_accesos_pobla', height=5, kind='reg')

    C:\Users\cuervoca\Anaconda3\lib\site-packages\scipy\stats\stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
      return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval

Out[28]:

    <seaborn.axisgrid.PairGrid at 0x1e1f3e93048>

Como se puede observar la correlación entre las variables es muy fuerte, razón por la cual es posible hacer el supuesto de que según como sea el comportamiento de accesos a la web, de manera similar será la proporción de la población que usará el computador o usará internet, dos indicadores clave para medir la penetración digital. Es así como la data de postdata se vuelve un insumo valiso, considerando además que los datos de uso de computador e internet son pobres en calidad, y disponibilidad para todas los departamentos del pais, mientras que la data de accesos si está disponible. Los accesos a la web sobre población será la variable dependiente de nuestro modelo para medir adopción de la tecnología en los departamentos.

# Modelo de regresión:¶

In [29]:

    #Cargue de información para realizar la regresión
    ruta_reg= os.path.join(os.getcwd(),'Data_limpia','Variables_regresion.xlsx')
    data_reg= pd.read_excel(ruta_reg, sheet='Hoja1')

In [30]:

    #Se crean variables dummy para manejar los departamentos:
    data_reg_dummy = pd.get_dummies(data_reg, columns=['departamento'])

In [31]:

    #Ahora se elimina la dummy del "departamento" Total nacional, ya que será absorbido por el beta cero de la regresión
    #al considerar que si todas las variables dummies son cero, significa que será el valor del total nacional. Se eliminan
    #otras variables que no serán usadas en esta etapa del proceso.
    data_reg_dummy = data_reg_dummy.drop(['anio', 'concatenado','departamento_Total Nacional','y_accesos'], axis=1)

In [32]:

    #Observamos la data:
    data_reg_dummy.head(10)

Out[32]:

| pib_percapita | desempleo | poblacion | ninios_mujer | esperanza_vida | tasa_migracion | mortalidad | natalidad | mortalidad_infantil | gini | ... | departamento_Nariño | departamento_Norte de Santander | departamento_Putumayo | departamento_Quindío | departamento_Risaralda | departamento_San Andrés, Providencia y Santa Catalina | departamento_Santander | departamento_Sucre | departamento_Tolima | departamento_Valle del Cauca
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
0 | 8.976426e+06 | 0.120847 | 5682310 | 0.334538 | 73.41 | 1.32 | 6.061325 | 17.753422 | 19.5 | 0.555 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
1 | 9.948822e+06 | 0.128856 | 5757973 | 0.334538 | 73.41 | 1.32 | 6.061325 | 17.753422 | 19.5 | 0.569 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
2 | 1.105656e+07 | 0.111717 | 5834865 | 0.334538 | 73.41 | 1.32 | 6.061325 | 17.753422 | 19.5 | 0.567 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
3 | 1.157509e+07 | 0.121151 | 5911399 | 0.334538 | 73.41 | 1.32 | 6.061325 | 17.753422 | 19.5 | 0.580 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
4 | 1.190752e+07 | 0.127609 | 5988552 | 0.334538 | 73.41 | 1.32 | 6.061325 | 17.753422 | 19.5 | 0.560 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
5 | 1.252585e+07 | 0.119197 | 6066003 | 0.311028 | 75.07 | 1.26 | 5.892037 | 17.088941 | 17.0 | 0.561 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
6 | 1.392916e+07 | 0.103967 | 6143809 | 0.311028 | 75.07 | 1.26 | 5.892037 | 17.088941 | 17.0 | 0.537 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
7 | 1.483779e+07 | 0.106295 | 6221817 | 0.311028 | 75.07 | 1.26 | 5.892037 | 17.088941 | 17.0 | 0.529 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
8 | 1.564651e+07 | 0.097272 | 6299990 | 0.311028 | 75.07 | 1.26 | 5.892037 | 17.088941 | 17.0 | 0.532 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
9 | 1.674769e+07 | 0.094151 | 6378132 | 0.311028 | 75.07 | 1.26 | 5.892037 | 17.088941 | 17.0 | 0.555 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0

10 rows × 41 columns

In [33]:

    #separamos la data:
    X= data_reg_dummy.copy()
    X= X.drop(['y_accesos_pobla'], axis=1)
    Y= data_reg_dummy['y_accesos_pobla']

In [34]:

    #Verificamos correlaciones de variables independientes para evadir la colinealidad:
    X_corr = X.copy()
    X_corr = X_corr[['pib_percapita', 'desempleo', 'poblacion', 'ninios_mujer','esperanza_vida',
                     'tasa_migracion', 'mortalidad', 'natalidad','mortalidad_infantil', 'gini',
                     'tasa_min_exp', 'trm']]

In [36]:

    dframe_xcorr=X_corr.corr(method='pearson')
    dframe_xcorr

Out[36]:

| pib_percapita | desempleo | poblacion | ninios_mujer | esperanza_vida | tasa_migracion | mortalidad | natalidad | mortalidad_infantil | gini | tasa_min_exp | trm
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
pib_percapita | 1.000000 | -0.139654 | 0.070907 | -0.281939 | 0.088090 | 0.391918 | -0.112512 | -0.234629 | -0.174576 | -0.330453 | -0.205932 | 0.227080
desempleo | -0.139654 | 1.000000 | -0.026879 | 0.118923 | -0.242824 | -0.213383 | 0.218698 | 0.099687 | 0.080413 | 0.277650 | 0.111290 | -0.128820
poblacion | 0.070907 | -0.026879 | 1.000000 | -0.203967 | 0.254513 | 0.154748 | -0.117703 | -0.179300 | -0.245991 | 0.147649 | -0.008380 | 0.012602
ninios_mujer | -0.281939 | 0.118923 | -0.203967 | 1.000000 | -0.683529 | -0.342578 | -0.107237 | 0.979040 | 0.841676 | 0.446362 | 0.146219 | -0.143626
esperanza_vida | 0.088090 | -0.242824 | 0.254513 | -0.683529 | 1.000000 | 0.312393 | -0.282453 | -0.677425 | -0.744201 | -0.260239 | -0.159946 | 0.196374
tasa_migracion | 0.391918 | -0.213383 | 0.154748 | -0.342578 | 0.312393 | 1.000000 | -0.389159 | -0.217229 | -0.322784 | -0.266791 | -0.067404 | 0.079799
mortalidad | -0.112512 | 0.218698 | -0.117703 | -0.107237 | -0.282453 | -0.389159 | 1.000000 | -0.209179 | -0.053189 | 0.020254 | 0.000158 | 0.077178
natalidad | -0.234629 | 0.099687 | -0.179300 | 0.979040 | -0.677425 | -0.217229 | -0.209179 | 1.000000 | 0.852978 | 0.435731 | 0.114993 | -0.152399
mortalidad_infantil | -0.174576 | 0.080413 | -0.245991 | 0.841676 | -0.744201 | -0.322784 | -0.053189 | 0.852978 | 1.000000 | 0.391782 | 0.069501 | -0.090942
gini | -0.330453 | 0.277650 | 0.147649 | 0.446362 | -0.260239 | -0.266791 | 0.020254 | 0.435731 | 0.391782 | 1.000000 | 0.169158 | -0.380703
tasa_min_exp | -0.205932 | 0.111290 | -0.008380 | 0.146219 | -0.159946 | -0.067404 | 0.000158 | 0.114993 | 0.069501 | 0.169158 | 1.000000 | 0.029001
trm | 0.227080 | -0.128820 | 0.012602 | -0.143626 | 0.196374 | 0.079799 | 0.077178 | -0.152399 | -0.090942 | -0.380703 | 0.029001 | 1.000000

El analisis muestra que hay una alta correlación entre las variables natalidad, mortalidad_infantil y niños por mujer. Por lo que será necesario transformarlas o eliminarlas.

In [37]:

    #Eliminamos ninios_mujer y obtenemos las correlaciones:
    X_corr = X_corr.drop(['ninios_mujer'],axis=1)
    dframe_xcorr1=X_corr.corr(method='pearson')
    dframe_xcorr1

Out[37]:

| pib_percapita | desempleo | poblacion | esperanza_vida | tasa_migracion | mortalidad | natalidad | mortalidad_infantil | gini | tasa_min_exp | trm
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
pib_percapita | 1.000000 | -0.139654 | 0.070907 | 0.088090 | 0.391918 | -0.112512 | -0.234629 | -0.174576 | -0.330453 | -0.205932 | 0.227080
desempleo | -0.139654 | 1.000000 | -0.026879 | -0.242824 | -0.213383 | 0.218698 | 0.099687 | 0.080413 | 0.277650 | 0.111290 | -0.128820
poblacion | 0.070907 | -0.026879 | 1.000000 | 0.254513 | 0.154748 | -0.117703 | -0.179300 | -0.245991 | 0.147649 | -0.008380 | 0.012602
esperanza_vida | 0.088090 | -0.242824 | 0.254513 | 1.000000 | 0.312393 | -0.282453 | -0.677425 | -0.744201 | -0.260239 | -0.159946 | 0.196374
tasa_migracion | 0.391918 | -0.213383 | 0.154748 | 0.312393 | 1.000000 | -0.389159 | -0.217229 | -0.322784 | -0.266791 | -0.067404 | 0.079799
mortalidad | -0.112512 | 0.218698 | -0.117703 | -0.282453 | -0.389159 | 1.000000 | -0.209179 | -0.053189 | 0.020254 | 0.000158 | 0.077178
natalidad | -0.234629 | 0.099687 | -0.179300 | -0.677425 | -0.217229 | -0.209179 | 1.000000 | 0.852978 | 0.435731 | 0.114993 | -0.152399
mortalidad_infantil | -0.174576 | 0.080413 | -0.245991 | -0.744201 | -0.322784 | -0.053189 | 0.852978 | 1.000000 | 0.391782 | 0.069501 | -0.090942
gini | -0.330453 | 0.277650 | 0.147649 | -0.260239 | -0.266791 | 0.020254 | 0.435731 | 0.391782 | 1.000000 | 0.169158 | -0.380703
tasa_min_exp | -0.205932 | 0.111290 | -0.008380 | -0.159946 | -0.067404 | 0.000158 | 0.114993 | 0.069501 | 0.169158 | 1.000000 | 0.029001
trm | 0.227080 | -0.128820 | 0.012602 | 0.196374 | 0.079799 | 0.077178 | -0.152399 | -0.090942 | -0.380703 | 0.029001 | 1.000000

In [38]:

    #Aun sigue habiendo demasiada correlación entre natalidad y mortalidad infantil.
    #Así que se buscará combinar estas dos variables como mortalidad_infantil/natalidad:
    X_corr['mortalidad_ratio']= X_corr['mortalidad_infantil']/X_corr['natalidad']
    X_corr = X_corr.drop(['mortalidad_infantil','natalidad'],axis=1)

In [39]:

    #obtenemos las correlaciones:
    dframe_xcorr2=X_corr.corr(method='pearson')
    dframe_xcorr2

Out[39]:

| pib_percapita | desempleo | poblacion | esperanza_vida | tasa_migracion | mortalidad | gini | tasa_min_exp | trm | mortalidad_ratio
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
pib_percapita | 1.000000 | -0.139654 | 0.070907 | 0.088090 | 0.391918 | -0.112512 | -0.330453 | -0.205932 | 0.227080 | -0.120775
desempleo | -0.139654 | 1.000000 | -0.026879 | -0.242824 | -0.213383 | 0.218698 | 0.277650 | 0.111290 | -0.128820 | 0.013684
poblacion | 0.070907 | -0.026879 | 1.000000 | 0.254513 | 0.154748 | -0.117703 | 0.147649 | -0.008380 | 0.012602 | -0.266816
esperanza_vida | 0.088090 | -0.242824 | 0.254513 | 1.000000 | 0.312393 | -0.282453 | -0.260239 | -0.159946 | 0.196374 | -0.674509
tasa_migracion | 0.391918 | -0.213383 | 0.154748 | 0.312393 | 1.000000 | -0.389159 | -0.266791 | -0.067404 | 0.079799 | -0.269143
mortalidad | -0.112512 | 0.218698 | -0.117703 | -0.282453 | -0.389159 | 1.000000 | 0.020254 | 0.000158 | 0.077178 | -0.008997
gini | -0.330453 | 0.277650 | 0.147649 | -0.260239 | -0.266791 | 0.020254 | 1.000000 | 0.169158 | -0.380703 | 0.274927
tasa_min_exp | -0.205932 | 0.111290 | -0.008380 | -0.159946 | -0.067404 | 0.000158 | 0.169158 | 1.000000 | 0.029001 | 0.036567
trm | 0.227080 | -0.128820 | 0.012602 | 0.196374 | 0.079799 | 0.077178 | -0.380703 | 0.029001 | 1.000000 | -0.048878
mortalidad_ratio | -0.120775 | 0.013684 | -0.266816 | -0.674509 | -0.269143 | -0.008997 | 0.274927 | 0.036567 | -0.048878 | 1.000000

In [40]:

    #Ahora se hace el ajuste en X para proceder con el modelo:
    X['mortalidad_ratio']= X['mortalidad_infantil']/X['natalidad']
    X= X.drop(['ninios_mujer','mortalidad_infantil','natalidad'],axis=1)

In [41]:

    #Implementamos la regresion lineal:
    from sklearn.linear_model import LinearRegression

    model = LinearRegression().fit(X,Y)
    r_sq = model.score(X, Y)
    print('coefficient of determination:', round(r_sq,4))

    coefficient of determination: 0.9316

In [42]:

    #Ahora obtenemos el R2 del modelo, los coeficientes de las variables y su significancia:
    import statsmodels.api as sm
    X2 = sm.add_constant(X)
    est = sm.OLS(Y, X2)
    est2 = est.fit()
    est2.summary()

Out[42]:

OLS Regression Results
Dep. Variable: | y_accesos_pobla | R-squared: | 0.932
Model: | OLS | Adj. R-squared: | 0.925
Method: | Least Squares | F-statistic: | 142.0
Date: | Fri, 12 Jun 2020 | Prob (F-statistic): | 4.76e-206
Time: | 18:13:58 | Log-Likelihood: | 695.79
No. Observations: | 435 | AIC: | -1314.
Df Residuals: | 396 | BIC: | -1155.
Df Model: | 38 |  |
Covariance Type: | nonrobust |  |

| coef | std err | t | P>|t| | [0.025 | 0.975]
const | -5.4600 | 0.728 | -7.496 | 0.000 | -6.892 | -4.028
pib_percapita | 8.723e-09 | 8.78e-10 | 9.934 | 0.000 | 7e-09 | 1.04e-08
desempleo | -0.3156 | 0.135 | -2.336 | 0.020 | -0.581 | -0.050
poblacion | 1.136e-08 | 6.01e-09 | 1.890 | 0.059 | -4.54e-10 | 2.32e-08
esperanza_vida | 0.0658 | 0.008 | 8.654 | 0.000 | 0.051 | 0.081
tasa_migracion | -0.0122 | 0.003 | -3.963 | 0.000 | -0.018 | -0.006
mortalidad | 0.1153 | 0.019 | 6.043 | 0.000 | 0.078 | 0.153
gini | -0.4549 | 0.158 | -2.877 | 0.004 | -0.766 | -0.144
tasa_min_exp | -0.2191 | 0.166 | -1.321 | 0.187 | -0.545 | 0.107
trm | 4.71e-05 | 7.61e-06 | 6.192 | 0.000 | 3.21e-05 | 6.21e-05
departamento_Antioquia | 0.6575 | 0.242 | 2.719 | 0.007 | 0.182 | 1.133
departamento_Arauca | 0.7600 | 0.278 | 2.737 | 0.006 | 0.214 | 1.306
departamento_Atlántico | 0.6585 | 0.264 | 2.491 | 0.013 | 0.139 | 1.178
departamento_Bogotá D.C. | 0.6571 | 0.238 | 2.766 | 0.006 | 0.190 | 1.124
departamento_Bolívar | 0.7604 | 0.268 | 2.842 | 0.005 | 0.234 | 1.287
departamento_Boyacá | 0.1981 | 0.277 | 0.714 | 0.476 | -0.347 | 0.744
departamento_Caldas | 0.3151 | 0.279 | 1.129 | 0.259 | -0.233 | 0.864
departamento_Caquetá | 0.8050 | 0.273 | 2.947 | 0.003 | 0.268 | 1.342
departamento_Casanare | 0.7830 | 0.276 | 2.835 | 0.005 | 0.240 | 1.326
departamento_Cauca | 0.8412 | 0.273 | 3.078 | 0.002 | 0.304 | 1.378
departamento_Cesar | 0.7106 | 0.272 | 2.610 | 0.009 | 0.175 | 1.246
departamento_Chocó | 1.0647 | 0.282 | 3.776 | 0.000 | 0.510 | 1.619
departamento_Cundinamarca | 0.5737 | 0.262 | 2.189 | 0.029 | 0.058 | 1.089
departamento_Córdoba | 0.6600 | 0.267 | 2.475 | 0.014 | 0.136 | 1.184
departamento_Huila | 0.6475 | 0.270 | 2.396 | 0.017 | 0.116 | 1.179
departamento_La Guajira | 0.6881 | 0.271 | 2.540 | 0.011 | 0.156 | 1.221
departamento_Magdalena | 0.5316 | 0.273 | 1.945 | 0.053 | -0.006 | 1.069
departamento_Meta | 0.8219 | 0.271 | 3.032 | 0.003 | 0.289 | 1.355
departamento_Nariño | 0.8465 | 0.270 | 3.138 | 0.002 | 0.316 | 1.377
departamento_Norte de Santander | 0.6027 | 0.271 | 2.226 | 0.027 | 0.070 | 1.135
departamento_Putumayo | 0.6480 | 0.276 | 2.351 | 0.019 | 0.106 | 1.190
departamento_Quindío | 0.4517 | 0.279 | 1.617 | 0.107 | -0.098 | 1.001
departamento_Risaralda | 0.5306 | 0.276 | 1.922 | 0.055 | -0.012 | 1.073
departamento_San Andrés, Providencia y Santa Catalina | 0.3687 | 0.280 | 1.314 | 0.189 | -0.183 | 0.920
departamento_Santander | 0.4681 | 0.270 | 1.734 | 0.084 | -0.063 | 0.999
departamento_Sucre | 0.4943 | 0.274 | 1.803 | 0.072 | -0.045 | 1.033
departamento_Tolima | 0.3750 | 0.275 | 1.365 | 0.173 | -0.165 | 0.915
departamento_Valle del Cauca | 0.4572 | 0.254 | 1.802 | 0.072 | -0.042 | 0.956
mortalidad_ratio | -0.3525 | 0.071 | -4.951 | 0.000 | -0.492 | -0.213

Omnibus: | 6.673 | Durbin-Watson: | 0.710
Prob(Omnibus): | 0.036 | Jarque-Bera (JB): | 6.845
Skew: | -0.234 | Prob(JB): | 0.0326
Kurtosis: | 3.397 | Cond. No. | 9.18e+09

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 9.18e+09. This might indicate that there are
strong multicollinearity or other numerical problems.
In [43]:

    #Se observa un R2 bueno, al parecer la variable de tasa minima de expamsion
    #no es significativa, por lo cual se elimina:
    X= X.drop(['tasa_min_exp'],axis=1)

In [44]:

    #Se corre un segundo modelo:
    model1 = LinearRegression().fit(X,Y)
    r_sq1 = model1.score(X, Y)
    print('coefficient of determination:',  round(r_sq1,4))

    coefficient of determination: 0.9313

In [45]:

    #Ahora obtenemos el R2 del modelo, los coeficientes de las variables y su significancia:
    import statsmodels.api as sm
    X2 = sm.add_constant(X)
    est = sm.OLS(Y, X2)
    est2 = est.fit()
    est2.summary()

Out[45]:

OLS Regression Results
Dep. Variable: | y_accesos_pobla | R-squared: | 0.931
Model: | OLS | Adj. R-squared: | 0.925
Method: | Least Squares | F-statistic: | 145.5
Date: | Fri, 12 Jun 2020 | Prob (F-statistic): | 9.30e-207
Time: | 18:14:58 | Log-Likelihood: | 694.83
No. Observations: | 435 | AIC: | -1314.
Df Residuals: | 397 | BIC: | -1159.
Df Model: | 37 |  |
Covariance Type: | nonrobust |  |

| coef | std err | t | P>|t| | [0.025 | 0.975]
const | -5.6664 | 0.712 | -7.958 | 0.000 | -7.066 | -4.266
pib_percapita | 8.884e-09 | 8.7e-10 | 10.209 | 0.000 | 7.17e-09 | 1.06e-08
desempleo | -0.3147 | 0.135 | -2.327 | 0.020 | -0.581 | -0.049
poblacion | 1.143e-08 | 6.01e-09 | 1.901 | 0.058 | -3.92e-10 | 2.33e-08
esperanza_vida | 0.0686 | 0.007 | 9.369 | 0.000 | 0.054 | 0.083
tasa_migracion | -0.0119 | 0.003 | -3.864 | 0.000 | -0.018 | -0.006
mortalidad | 0.1159 | 0.019 | 6.072 | 0.000 | 0.078 | 0.153
gini | -0.4812 | 0.157 | -3.065 | 0.002 | -0.790 | -0.173
trm | 4.288e-05 | 6.91e-06 | 6.206 | 0.000 | 2.93e-05 | 5.65e-05
departamento_Antioquia | 0.6595 | 0.242 | 2.725 | 0.007 | 0.184 | 1.135
departamento_Arauca | 0.7762 | 0.278 | 2.796 | 0.005 | 0.230 | 1.322
departamento_Atlántico | 0.6586 | 0.265 | 2.489 | 0.013 | 0.138 | 1.179
departamento_Bogotá D.C. | 0.6496 | 0.238 | 2.733 | 0.007 | 0.182 | 1.117
departamento_Bolívar | 0.7646 | 0.268 | 2.855 | 0.005 | 0.238 | 1.291
departamento_Boyacá | 0.2022 | 0.278 | 0.728 | 0.467 | -0.344 | 0.748
departamento_Caldas | 0.3216 | 0.279 | 1.152 | 0.250 | -0.227 | 0.871
departamento_Caquetá | 0.8229 | 0.273 | 3.014 | 0.003 | 0.286 | 1.360
departamento_Casanare | 0.7925 | 0.276 | 2.868 | 0.004 | 0.249 | 1.336
departamento_Cauca | 0.8526 | 0.273 | 3.119 | 0.002 | 0.315 | 1.390
departamento_Cesar | 0.7172 | 0.272 | 2.632 | 0.009 | 0.182 | 1.253
departamento_Chocó | 1.0866 | 0.282 | 3.857 | 0.000 | 0.533 | 1.640
departamento_Cundinamarca | 0.5738 | 0.262 | 2.188 | 0.029 | 0.058 | 1.090
departamento_Córdoba | 0.6664 | 0.267 | 2.497 | 0.013 | 0.142 | 1.191
departamento_Huila | 0.6560 | 0.270 | 2.426 | 0.016 | 0.124 | 1.188
departamento_La Guajira | 0.6917 | 0.271 | 2.551 | 0.011 | 0.159 | 1.225
departamento_Magdalena | 0.5372 | 0.274 | 1.963 | 0.050 | -0.001 | 1.075
departamento_Meta | 0.8272 | 0.271 | 3.049 | 0.002 | 0.294 | 1.361
departamento_Nariño | 0.8526 | 0.270 | 3.158 | 0.002 | 0.322 | 1.383
departamento_Norte de Santander | 0.6128 | 0.271 | 2.262 | 0.024 | 0.080 | 1.145
departamento_Putumayo | 0.6616 | 0.276 | 2.400 | 0.017 | 0.120 | 1.204
departamento_Quindío | 0.4571 | 0.280 | 1.635 | 0.103 | -0.093 | 1.007
departamento_Risaralda | 0.5343 | 0.276 | 1.934 | 0.054 | -0.009 | 1.077
departamento_San Andrés, Providencia y Santa Catalina | 0.3722 | 0.281 | 1.326 | 0.186 | -0.180 | 0.924
departamento_Santander | 0.4694 | 0.270 | 1.737 | 0.083 | -0.062 | 1.001
departamento_Sucre | 0.5000 | 0.274 | 1.822 | 0.069 | -0.039 | 1.039
departamento_Tolima | 0.3839 | 0.275 | 1.397 | 0.163 | -0.156 | 0.924
departamento_Valle del Cauca | 0.4582 | 0.254 | 1.804 | 0.072 | -0.041 | 0.957
mortalidad_ratio | -0.3501 | 0.071 | -4.915 | 0.000 | -0.490 | -0.210

Omnibus: | 5.326 | Durbin-Watson: | 0.720
Prob(Omnibus): | 0.070 | Jarque-Bera (JB): | 5.315
Skew: | -0.207 | Prob(JB): | 0.0701
Kurtosis: | 3.350 | Cond. No. | 9.17e+09

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 9.17e+09. This might indicate that there are
strong multicollinearity or other numerical problems.

El R2 sigue siendo bueno, y se observa que todas las variables no instrumentales son significativas.

In [46]:

    #Se compara la proyección con el real:
    Y_pred= model1.predict(X)
    compare= pd.DataFrame()
    compare['curva']=Y_pred
    compare['tipo']='predicción'
    compare['x']=np.arange(0,435)

    compare1= pd.DataFrame()
    compare1['curva']=Y
    compare1['tipo']='real'
    compare1['x']=np.arange(0,435)

    compare= pd.concat([compare,compare1])
    sns.lineplot(x= 'x',y='curva', hue='tipo',data=compare)

Out[46]:

    <matplotlib.axes._subplots.AxesSubplot at 0x1e1f449ac88>

Al comparar la preducción con el real se obser que

In [47]:

    #Ahora varificamos la distribución de los residuales:
    errors= Y_pred - Y
    errors_x= [number for number in range(0,len(errors))]

    #Ploteamos los residuales:
    plt.scatter(errors_x,errors)
    z = np.polyfit(errors_x,errors, 1)
    p = np.poly1d(z)
    plt.plot(errors_x,p(errors_x),"r--")
    plt.show()

Cómo se observa no hay tendencia alguna de los residuales, se distribuyen de manera aleatoria alrededor del cero, por lo que se considera que es una buena aproximación de la regresión y el modelo se ajusta a los datos.

# Modelo de clustering para determinar si existen departamentos similares¶

Después de identificar las variables significativas en el modelo de regresión y comprobar que hay relación entre ellas, ahora se procede a tratar de encontrar una relación entre las regiones que nos indique si hay clusters de tipos de regiones y si ese comportamiento se mantiene por cada año o no.

In [48]:

    data_clusters= data_reg.copy()

In [49]:

    #Elegir sólo data de los departamentos, quitar el total nacional:
    data_clusters=data_clusters.loc[(data_clusters.departamento != 'Total Nacional')]

    #Eliminamos las variables que no son necesarias y las que ya encontramos que no son significativas, así como la
    #creación de una nueva.
    data_clusters['mortalidad_ratio'] = data_clusters['mortalidad_infantil']/data_clusters['natalidad']
    data_clusters= data_clusters.drop(['ninios_mujer','mortalidad_infantil','natalidad','departamento',
                                       'anio','tasa_min_exp','y_accesos'],axis=1)

In [50]:

    #Usados PCA para obtener los componentes más importantes y poder graficar en 2D:
    from sklearn.decomposition import PCA
    data_pca=data_clusters[['pib_percapita', 'desempleo', 'poblacion', 'esperanza_vida', 'tasa_migracion', 'mortalidad',
                  'gini', 'trm','y_accesos_pobla', 'mortalidad_ratio']]

In [51]:

    #Ahora estandarizamos las variables para quitar efecto de magnitudes:
    from sklearn import preprocessing
    data_pca_std = preprocessing.scale(data_pca)
    componts=PCA(n_components=2)
    componts.fit(data_pca_std)
    trans=componts.transform(data_pca_std) #Aplica la reducción de dimensionalidad
    print (componts.explained_variance_ratio_) #Varianza explicada por los componentes

    [0.3463046  0.15195972]

Se observa que con sólo dos componentes se obtiene el 50% de la información.

In [52]:

    #Ploteamos los datos:
    plt.scatter(trans[:, 0], trans[:, 1])

Out[52]:

    <matplotlib.collections.PathCollection at 0x1e1b37e9278>

Pareciera que PCA no es la mejor manera de obtner los componentes más importantes porque la división no es limpia así que se probará con t-SNE:

In [53]:

    from sklearn.manifold import TSNE
    tcomponts= TSNE(n_components=2,perplexity=25).fit_transform(data_pca_std)
    #Se obtiene el número de componentes y como se ven:
    tcomponts.shape
    #Datos transformados para la primera compoente:
    tcomponts[:,0]
    #Graficar la iformación usando T-SNE:
    plt.scatter(tcomponts[:,0],tcomponts[:,1])

Out[53]:

    <matplotlib.collections.PathCollection at 0x1e1b7e83438>

Con t-SNE se obserba más claro que hay ciertos grupos diferentes, pero se debe realizar clustering para poderlos dierenciar claramente.

In [54]:

    #Implementar cluster de skalearn:
    from sklearn.cluster import KMeans
    #Se implementa el método del codo para obtener el No. de K, así como su visualización:
    #Método del codo:
    sse=[] #Se crea una lista que contenga la suma de los errores en cada iteración
    iterar=10 #número de iteraciones a considerar
    rango = range(1, iterar)
    for i in rango:
        kmeans = KMeans(n_clusters = i, init = 'k-means++',random_state = 0)
        kmeans.fit(tcomponts)
        sse.append(kmeans.inertia_) #Con este atributo se obtiene la suma de distancias al cuadrado de cada punto a su centroide.

In [55]:

    #Graficar el codo:
    plt.scatter(rango,sse)
    plt.plot(rango,sse)
    plt.xlabel('K grupos')
    plt.ylabel('SSE')

Out[55]:

    Text(0,0.5,'SSE')

Se debe probar con 4 o 5 grupos:

In [90]:

    #Con 4 grupos:
    n=4
    km=KMeans(n_clusters=n,max_iter=1000).fit(tcomponts)
    y_km=KMeans(n_clusters=n,max_iter=1000).fit_predict(tcomponts)
    km.cluster_centers_
    plt.scatter(tcomponts[:, 0], tcomponts[:, 1],c=y_km)
    plt.title("Clusters generados")

Out[90]:

    Text(0.5,1,'Clusters generados')

In [57]:

    #Con 5 grupos:
    n=5
    km=KMeans(n_clusters=n,max_iter=1000).fit(tcomponts)
    y_km=KMeans(n_clusters=n,max_iter=1000).fit_predict(tcomponts)
    km.cluster_centers_
    plt.scatter(tcomponts[:, 0], tcomponts[:, 1],c=y_km)
    plt.title("Clusters generados")

Out[57]:

    Text(0.5,1,'Clusters generados')

In [58]:

    #K-means no es calro así que hacemos uso de un modelo por densidad:
    from sklearn.cluster import DBSCAN
    neps=5
    sampls=10
    clusters = DBSCAN(eps=neps, min_samples=sampls).fit(tcomponts)
    y_clusters = DBSCAN(eps=neps, min_samples=sampls).fit_predict(tcomponts)
    plt.scatter(tcomponts[:, 0], tcomponts[:, 1],c=y_clusters)
    plt.xlabel('T-SNE1')
    plt.ylabel('T-SNE2')

Out[58]:

    Text(0,0.5,'T-SNE2')

In [59]:

    #Probamos otra configuración:
    neps=6
    sampls=10
    clusters = DBSCAN(eps=neps, min_samples=sampls).fit(tcomponts)
    y_clusters = DBSCAN(eps=neps, min_samples=sampls).fit_predict(tcomponts)
    plt.scatter(tcomponts[:, 0], tcomponts[:, 1],c=y_clusters)
    plt.xlabel('T-SNE1')
    plt.ylabel('T-SNE2')

Out[59]:

    Text(0,0.5,'T-SNE2')

Con esta información creemos que la mejor configuración será de 4 clusters. Así que se vuelve a correr para guardar esa información.

In [91]:

    #Ahora se relaciona la data original con los clusters:
    data_clusters['Cluster']= y_km

In [ ]:

    #Se descarga para ser visualizada en tableau::
    data_clusters.to_excel(os.path.join(os.getcwd(),'Data_tableau','clusters.xlsx'), index= False)

# Arbol de decisión para encontrar las variables más importantes:¶

In [63]:

    #Se importan las librerías pertinentes:
    from sklearn.model_selection import train_test_split
    from sklearn.tree import DecisionTreeClassifier

In [92]:

    #Se selecciona la data a usar:
    data_tree= data_clusters[['pib_percapita', 'desempleo', 'poblacion','esperanza_vida','tasa_migracion',
                              'mortalidad', 'gini', 'trm','mortalidad_ratio','Cluster']]

    #Se separa x y Y(será la clasificación de cluster obtenida anteriormente):
    data_tree_x= data_tree.copy()
    data_tree_x= data_tree_x.drop(['Cluster'],axis=1)
    data_tree_y= data_tree['Cluster']

    #Separamos los datos en entrenamiento y test:
    xtrain,xtest,ytrain,ytest = train_test_split(data_tree_x, data_tree_y)

    #Toma el 75% de los datos para entrenar
    print (xtrain.shape)

    (315, 9)

In [93]:

    data_tree['Cluster'].unique()

Out[93]:

    array([2, 1, 0, 3], dtype=int64)

In [94]:

    #Se implementa el algoritmo:
    arbol=DecisionTreeClassifier()
    arbol.fit(xtrain,ytrain)

Out[94]:

    DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=None,
                max_features=None, max_leaf_nodes=None,
                min_impurity_decrease=0.0, min_impurity_split=None,
                min_samples_leaf=1, min_samples_split=2,
                min_weight_fraction_leaf=0.0, presort=False, random_state=None,
                splitter='best')

In [95]:

    #Se obtiene el resultado sobre el set de test:
    print ("primera iteración sobre test:", round(arbol.score(xtest,ytest),4)*100) #Capacidad predictiva
    #Se verifica sobre el set de entrenamiento:
    print ("primera iteración sobre train:",round(arbol.score(xtrain,ytrain),4)*100)

    primera iteración sobre test: 99.05000000000001
    primera iteración sobre train: 100.0

Se observa que hay sobre ajuste al tener un score del 100% dobre train, así que se modifica la profundidad del arbol:

In [96]:

    #Para evitar el sobre ajuste se indica una profundidad inferior:
    arbol=DecisionTreeClassifier(max_depth=5)
    arbol.fit(xtrain,ytrain)
    #Se obtiene el resultado sobre el set de test:
    print ("segunda iteración sobre test:",round(arbol.score(xtest,ytest),4)*100) #Capacidad predictiva
    #Se verifica sobre el set de entrenamiento:
    print ("segunda iteración sobre train:",round(arbol.score(xtrain,ytrain),4)*100)

    segunda iteración sobre test: 98.1
    segunda iteración sobre train: 99.68

In [97]:

    #Mostrar las características más importantes:
    targetnames= data_tree_x.columns
    caract=data_tree_x.shape[1]
    plt.barh(range(caract),arbol.feature_importances_)
    plt.yticks(np.arange(caract),targetnames)
    plt.xlabel('Importancia de las características')
    plt.ylabel('Características')
    plt.show()

Según el resultado se observa que la mortalidad_ratio, la tasa de migración y el pib_percapita son variables importantes al momento de diferenciar los custers. Además son variables que se podrían afectar de alguna manera para entender como mejorar la penetración digital, cosa que no se podría por ejemplo con la trm, que al parecer afecta, pero en su comportamiento interfieren demasiadas variables exogenas.

# Proyecciones de series temporales con Holt-Winters y regresiones¶

In [98]:

    #Se importa la librería de holt-winters:
    from statsmodels.tsa.holtwinters import ExponentialSmoothing

In [99]:

    #Obtenemos toda la data procesada hasta el momento:
    data_forecast= X.merge(data_reg)

In [101]:

    data_forecast.head(10)

Out[101]:

| pib_percapita | desempleo | poblacion | esperanza_vida | tasa_migracion | mortalidad | gini | trm | departamento_Antioquia | departamento_Arauca | ... | mortalidad_ratio | departamento | anio | concatenado | ninios_mujer | natalidad | mortalidad_infantil | tasa_min_exp | y_accesos | y_accesos_pobla
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
0 | 8.976426e+06 | 0.120847 | 5682310 | 73.41 | 1.32 | 6.061325 | 0.555 | 2321.494767 | 1 | 0 | ... | 1.098380 | Antioquia | 2005 | Antioquia2005 | 0.334538 | 17.753422 | 19.5 | 0.063571 | 1754013 | 0.308680
1 | 9.948822e+06 | 0.128856 | 5757973 | 73.41 | 1.32 | 6.061325 | 0.569 | 2358.959123 | 1 | 0 | ... | 1.098380 | Antioquia | 2006 | Antioquia2006 | 0.334538 | 17.753422 | 19.5 | 0.064990 | 1734404 | 0.301218
2 | 1.105656e+07 | 0.111717 | 5834865 | 73.41 | 1.32 | 6.061325 | 0.567 | 2076.239753 | 1 | 0 | ... | 1.098380 | Antioquia | 2007 | Antioquia2007 | 0.334538 | 17.753422 | 19.5 | 0.087234 | 1773622 | 0.303970
3 | 1.157509e+07 | 0.121151 | 5911399 | 73.41 | 1.32 | 6.061325 | 0.580 | 1967.112268 | 1 | 0 | ... | 1.098380 | Antioquia | 2008 | Antioquia2008 | 0.334538 | 17.753422 | 19.5 | 0.098061 | 1695186 | 0.286766
4 | 1.190752e+07 | 0.127609 | 5988552 | 73.41 | 1.32 | 6.061325 | 0.560 | 2153.297699 | 1 | 0 | ... | 1.098380 | Antioquia | 2009 | Antioquia2009 | 0.334538 | 17.753422 | 19.5 | 0.057263 | 1852058 | 0.309266
5 | 1.252585e+07 | 0.119197 | 6066003 | 75.07 | 1.26 | 5.892037 | 0.561 | 1898.683671 | 1 | 0 | ... | 0.994795 | Antioquia | 2010 | Antioquia2010 | 0.311028 | 17.088941 | 17.0 | 0.031646 | 1846800 | 0.304451
6 | 1.392916e+07 | 0.103967 | 6143809 | 75.07 | 1.26 | 5.892037 | 0.537 | 1846.967096 | 1 | 0 | ... | 0.994795 | Antioquia | 2011 | Antioquia2011 | 0.311028 | 17.088941 | 17.0 | 0.039817 | 2474288 | 0.402729
7 | 1.483779e+07 | 0.106295 | 6221817 | 75.07 | 1.26 | 5.892037 | 0.529 | 1797.785765 | 1 | 0 | ... | 0.994795 | Antioquia | 2012 | Antioquia2012 | 0.311028 | 17.088941 | 17.0 | 0.049662 | 2906978 | 0.467223
8 | 1.564651e+07 | 0.097272 | 6299990 | 75.07 | 1.26 | 5.892037 | 0.532 | 1869.098904 | 1 | 0 | ... | 0.994795 | Antioquia | 2013 | Antioquia2013 | 0.311028 | 17.088941 | 17.0 | 0.034232 | 3281151 | 0.520818
9 | 1.674769e+07 | 0.094151 | 6378132 | 75.07 | 1.26 | 5.892037 | 0.555 | 2000.326822 | 1 | 0 | ... | 0.994795 | Antioquia | 2014 | Antioquia2014 | 0.311028 | 17.088941 | 17.0 | 0.038842 | 3524609 | 0.552608

10 rows × 46 columns

Debido a que hay valores negativos de la serie temporal de tasa de migración, se debe realizar un tratamiento especial para que el algoritmo sea capaz de aceptar el input, pero además procesar el output para que tenga sentido, es así como el input se eleva al cuadrado y la salida se le saca raiz, siempre que el valor anterior sea positivo, si el valor anterior fuera negativo se debe no sólo sacarle la raiz sino además multiplicar por -1 para que siga la tendencia negativa.

In [102]:

    #Se procede a obtener la predicción de la tasa de migración para cada departamento:
    #Se hace una prueba con Arauca que tiene un comportamiento especial:
    df=data_forecast.loc[(data_forecast.departamento=="Arauca"),'tasa_migracion']
    df= df.reset_index(drop=True)

    #Debido a que hay valores negativos de la serie temporal, se eleva al cuadrado para que el algoritmo lo pueda
    #trabajar:
    train= df.iloc[:14]
    test = df.iloc[11:]
    model = ExponentialSmoothing(train**2, seasonal='mul', seasonal_periods=11).fit()
    pred = model.predict(start=14, end= 19)

    #Posteriormente se trabaja la predicción para volverla a dejar en negativo y en la escala adecuada si es el caso:
    if np.sign(train[13]) > 0:
        pred = np.sqrt(pred)
    else:
        pred = np.sqrt(pred)*-1

    #Se muestra un ejemplo de la capacidad predictiva:
    plt.plot(train.index, train, label='Train')
    plt.plot(test.index, test, label='Test')
    plt.plot(pred.index, pred, label='Holt-Winters')
    plt.legend(loc='best')

Out[102]:

    <matplotlib.legend.Legend at 0x1e19c122e80>

Como se observa, así los valores sean negativos, el modelo es capaz de predecir a futuro la serie temporal.

In [103]:

    #Se procede a hacer la proyección de todos los departamentos y todas las variables. Para el pib percapita y la mortalidad ratio
    #debido a su naturaleza "suave" se decide proyectar con regresiones lineales:

    #Se crean los data frames que guardarán la información, para cada departamento, cada año y variable:
    depart= data_forecast.departamento.unique()
    pred_tasa_migracion=pd.DataFrame(columns=['tasa_migracion','departamento'])
    pred_pib_percapita=pd.DataFrame(columns=['pib_percapita','departamento'])
    pred_mortalidad_ratio=pd.DataFrame(columns=['mortalidad_ratio','departamento'])

    pred_tasa_migracion1=pd.DataFrame(columns=['tasa_migracion','departamento'])
    pred_pib_percapita1=pd.DataFrame(columns=['pib_percapita','departamento'])
    pred_mortalidad_ratio1=pd.DataFrame(columns=['mortalidad_ratio','departamento'])

    for i in range(0,len(depart)):
        df=data_forecast.loc[(data_forecast.departamento==depart[i]),'tasa_migracion']
        df= df.reset_index(drop=True)

        #Debido a que hay valores negativos de la serie temporal, se eleva al cuadrado para que el algoritmo lo pueda
        #trabajar:
        train= df.iloc[:15]
        model_tm = ExponentialSmoothing(train**2, seasonal='mul', seasonal_periods=11).fit()
        pred_tm = model_tm.predict(start=15, end= 19)
        #Posteriormente se trabaja la predicción para volverla a dejar en negativo y en la escala adecuada si es el caso:
        if np.sign(train[14]) > 0:
            pred_tm = np.sqrt(pred_tm)
        else:
            pred_tm = np.sqrt(pred_tm)*-1
        pred_tasa_migracion1['tasa_migracion']= pred_tm
        pred_tasa_migracion1['departamento']= depart[i]
        pred_tasa_migracion=pd.concat([pred_tasa_migracion,pred_tasa_migracion1])

        #Modelo para predecir el PIB_percapita:
        df2=data_forecast.loc[(data_forecast.departamento==depart[0]),'pib_percapita']
        df2= df2.reset_index(drop=True)
        train2= df2.iloc[:15]
        x_pp= pd.DataFrame(np.arange(0, 15))
        model_pp = LinearRegression().fit(x_pp,train2)
        pred_pp= model_pp.predict(pd.DataFrame(np.arange(15, 20)))
        pred_pib_percapita1['pib_percapita'] = pred_pp
        pred_pib_percapita1['departamento']= depart[i]
        pred_pib_percapita=pd.concat([pred_pib_percapita,pred_pib_percapita1])

        #Modelo para predecir la mortalidad_ratio:
        df3=data_forecast.loc[(data_forecast.departamento==depart[i]),'mortalidad_ratio']
        df3= df3.reset_index(drop=True)
        train3= df3.iloc[:15]
        x_mr= pd.DataFrame(np.arange(0, 15))
        model_mr = LinearRegression().fit(x_mr,train3)
        pred_mr= model_mr.predict(pd.DataFrame(np.arange(15, 20)))
        pred_mortalidad_ratio1['mortalidad_ratio']= pred_mr
        pred_mortalidad_ratio1['departamento']=depart[i]
        pred_mortalidad_ratio=pd.concat([pred_mortalidad_ratio,pred_mortalidad_ratio1])

In [104]:

    #Se junta toda la data:
    pred_tasa_migracion=pred_tasa_migracion.reset_index()
    pred_pib_percapita=pred_pib_percapita.reset_index()
    pred_mortalidad_ratio=pred_mortalidad_ratio.reset_index()

In [105]:

    data_proyecx= pred_tasa_migracion.merge(pred_pib_percapita[['pib_percapita']],left_index=True,right_index=True)
    data_proyecx= data_proyecx.merge(pred_mortalidad_ratio[['mortalidad_ratio']],left_index=True,right_index=True)
    data_proyecx.head(10)

Out[105]:

| index | tasa_migracion | departamento | pib_percapita | mortalidad_ratio
--- | --- | --- | --- | --- | ---
0 | 15 | 1.210000 | Antioquia | 2.270115e+07 | 0.862815
1 | 16 | 1.155000 | Antioquia | 2.364158e+07 | 0.845682
2 | 17 | 1.155000 | Antioquia | 2.458200e+07 | 0.828548
3 | 18 | 1.155000 | Antioquia | 2.552243e+07 | 0.811414
4 | 19 | 1.155000 | Antioquia | 2.646285e+07 | 0.794280
5 | 15 | -9.790000 | Arauca | 2.270115e+07 | 1.840653
6 | 16 | -9.202066 | Arauca | 2.364158e+07 | 1.847128
7 | 17 | -9.202066 | Arauca | 2.458200e+07 | 1.853602
8 | 18 | -9.202066 | Arauca | 2.552243e+07 | 1.860076
9 | 19 | -9.202066 | Arauca | 2.646285e+07 | 1.866551

In [ ]:

    #Se baja la información para realizar todo el ajuste y pruebas de aumento y disminución de valores de variables a futuro:
    data_proyecx.to_excel(os.path.join(os.getcwd(),'Data_tableau','proyecciones_variables_x.xlsx'), index= False)
    X.to_excel(os.path.join(os.getcwd(),'Data_tableau','X.xlsx'), index= False)

# Realizar proyecciones de penetración digital tomando como input las proyecciones de las variables independientes¶

Utilizar el primer modelo de regresión lineal hecho para proyectar los valores a futuro junto con las posibles combinaciones de variación de variables x. Habrá 343 escenarios diferentes, ya que son 7 posibles estados: 0%,+5,+10%,+15%,-5%,-10%,-15%; para 3 diferentes variables: mortalidad_ratio, tasa_migracion, pib_percapita => 7x7x7 = 343 escenarios.

Para obtener las diferentes combinaciones se usó el siguiente algoritmo:

lista_posi= [0,0.05,0.1,0.15,-0.05,-0.1,-0.15]

lista_completa=[]

for i in range(0,len(lista_posi)):

    for j in range(0,len(lista_posi)):

        for k in range(0,len(lista_posi)):

        lista_completa.append((lista_posi[i],lista_posi[j],lista_posi[k]))

lista dataframe=pd.DataFrame(lista_completa1)

lista_completa1=lista_completa for i in range(0,144): lista_completa1=lista_completa1+lista_completa

In [106]:

    #Se carga la data con los escenarios construidos:
    ruta_proy= os.path.join(os.getcwd(),'Data_limpia','X_escenarios.xlsx')
    data_proy_escenarios= pd.read_excel(ruta_proy, sheet='Sheet1')

In [107]:

    data_proy_escenarios.head(10)

Out[107]:

| anio | tipo | var_pib_percapita | var_tasa_migra | var_mortalidad_ratio | pib_percapita | desempleo | poblacion | esperanza_vida | tasa_migracion | ... | departamento_Nariño | departamento_Norte de Santander | departamento_Putumayo | departamento_Quindío | departamento_Risaralda | departamento_San Andrés, Providencia y Santa Catalina | departamento_Santander | departamento_Sucre | departamento_Tolima | departamento_Valle del Cauca
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
0 | 2005 | real | 0.0 | 0.0 | 0.0 | 8.976426e+06 | 0.120847 | 5682310 | 73.41 | 1.32 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
1 | 2006 | real | 0.0 | 0.0 | 0.0 | 9.948822e+06 | 0.128856 | 5757973 | 73.41 | 1.32 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
2 | 2007 | real | 0.0 | 0.0 | 0.0 | 1.105656e+07 | 0.111717 | 5834865 | 73.41 | 1.32 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
3 | 2008 | real | 0.0 | 0.0 | 0.0 | 1.157509e+07 | 0.121151 | 5911399 | 73.41 | 1.32 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
4 | 2009 | real | 0.0 | 0.0 | 0.0 | 1.190752e+07 | 0.127609 | 5988552 | 73.41 | 1.32 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
5 | 2010 | real | 0.0 | 0.0 | 0.0 | 1.252585e+07 | 0.119197 | 6066003 | 75.07 | 1.26 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
6 | 2011 | real | 0.0 | 0.0 | 0.0 | 1.392916e+07 | 0.103967 | 6143809 | 75.07 | 1.26 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
7 | 2012 | real | 0.0 | 0.0 | 0.0 | 1.483779e+07 | 0.106295 | 6221817 | 75.07 | 1.26 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
8 | 2013 | real | 0.0 | 0.0 | 0.0 | 1.564651e+07 | 0.097272 | 6299990 | 75.07 | 1.26 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
9 | 2014 | real | 0.0 | 0.0 | 0.0 | 1.674769e+07 | 0.094151 | 6378132 | 75.07 | 1.26 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0

10 rows × 43 columns

In [108]:

    #Obtenemos sólo la data que le corresponde al modelo:
    X_escenarios = data_proy_escenarios[X.columns]

In [109]:

    #Hacemos la predicción con el modelo entrenado anteriormente:
    pred_escenarios= model1.predict(X_escenarios)

In [110]:

    #Unir la predicción con la base original:
    data_proy_escenarios['y_accesos_predict']=pred_escenarios

In [111]:

    #Comparar las curvas de predicción de algunos deparamentos:
    #Ejemplo para Antioquia:
    plt.plot(data_proy_escenarios.loc[(data_proy_escenarios['departamento_Antioquia']==1) & (data_proy_escenarios['tipo']=='real'),'y_accesos_predict'], label='prediccion')
    plt.plot(data_proy_escenarios.loc[(data_proy_escenarios['departamento_Antioquia']==1) & (data_proy_escenarios['tipo']=='real'),'y_accesos_pobla'], label='real')
    plt.legend(loc='best')

Out[111]:

    <matplotlib.legend.Legend at 0x1e1a046cd68>

In [112]:

    #Ejemplo para Valle del Cauca:
    plt.plot(data_proy_escenarios.loc[(data_proy_escenarios['departamento_Valle del Cauca']==1) & (data_proy_escenarios['tipo']=='real'),'y_accesos_predict'], label='prediccion')
    plt.plot(data_proy_escenarios.loc[(data_proy_escenarios['departamento_Valle del Cauca']==1) & (data_proy_escenarios['tipo']=='real'),'y_accesos_pobla'], label='real')
    plt.legend(loc='best')

Out[112]:

    <matplotlib.legend.Legend at 0x1e19a5ceda0>

In [113]:

    #Verificamos de Accuracy:
    r_sq2 = model1.score(X_escenarios[0:435], data_proy_escenarios['y_accesos_pobla'][0:435])
    print('coefficient of determination:', round(r_sq2,4))

    coefficient of determination: 0.9313

Se observa que el modelo tiene un buen ajuste a los datos de penetración digital, sin necesidad de ser una línea recta sino que responde al comportamiento de la interacción de diferentes variables, permitiendo así ajustarle a cambios de pendiente.

In [ ]:

    #Por último se descargan los archivos para realizar las gráficas en tableau. Por favor continúe con el proceso
    #e interactúe con los dashboards que permiten predecir el comportamiento futuro de la penetración digital para los
    #siguientes 5 años!:
    data_proy_escenarios.to_excel(os.path.join(os.getcwd(),'Data_tableau','simulacion_escenarios.xlsx'), index= False)