My words, sometimes technical, sometimes not¶

2019-05-25
in estadistica, R
13 min read

Regresion Lineal - ISLR Capítulo 3

La regresión lineal simple es un método muy directo para estimar una variable cuantitativa Y en base a un solo predictor X. Asume que hay una relación lineal entre X e Y. $$ Y \approx \beta_0 + \beta_1X$$ $\beta_0$ y $\beta_1$ son dos constantes desconocidas que representan al intercepto y a la pendiente del modelo lineal. Son los coeficientes o parámetros. Con nuestros datos podemos estimar coeficientes para predecir futuros valores de Y basados en X y nuestro modelo.

Estimación de Coeficientes

Los coeficientes que buscamos son $\hat \beta_0$ y $\hat \beta_1$ (estimados, por eso el sombrero) son aquellos que generen una recta que pase lo más cerca posible de todos nuestros datos de entrenamiento. Hay varias manera de definir "cerca" pero la más usada es el enfoque de mínimos cuadrados.

Supongamos un ejemplo donde tenemos datos de horas trabajadas por ciertos individuos y la paga que reciben. Supongamos a fines del ejemplo que la relación entre salario y horas es lineal (sabemos que no es real...)

Cuando estimemos $\hat \beta_0$ y $\hat \beta_1$ obtendremos después un valor $\hat y_i$ para cada valor de $x_i$ (cada observación), que será el resultado de la predicción de nuestro modelo para ese valor de horas trabajadas. $$ \hat y_i = \hat \beta_0 + \hat \beta_1 x_i$$ Luego $e_i = y_i - \hat y_i$ representa el residuo, que es la diferencia entre el valor real del salario para esa observación y el valor que predice nuestro modelo.
Una métrica importante a saber es la suma de resiudos al cuadrado (RSS por siglas en inglés) que es: $$ RSS = e_1^2 + e_2^2 + ... + e_n^2$$ o de manera equivalente: $$ RSS = (y_1 - \hat \beta_0 - \hat \beta_1x_1) ^2 + (y_2 - \hat \beta_0 - \hat \beta_1x_2) ^2 + ... + (y_n - \hat \beta_0 - \hat \beta_1x_n) ^2 $$

Que basicamente es la suma de todas las diferencias entre lo predicho por nuestro modelo y el dato real de nuestro set, elevadas al cuadrado. Esto último es principalmente para evitar que se compensen los errores. Sobreestimar por 10 y luego subestimar por 10 tiene como suma de errores 0. Si elevamos esas diferencias al cuadrado, todos los errores serán positivos y se acumularán. En este caso seria $10^2$ + $(-10)^2$, que es 200.
El enfoque de mínimos cuadrados estima $\hat \beta_0$ y $\hat \beta_1$ de tal manera que el RSS sea el mínimo posible dados los datos.

Usando un poco de cálculo se puede demostrar que los parámetros que minimizan RSS son: $$ \hat \beta_1 = \frac{\sum_{i = 1}^n (x_i - \bar{x}) (y_i - \bar{y})}{\sum_{i=1}^n (x_i - \bar{x})^2}$$ $$ \hat \beta_0 = \bar{y} - \hat \beta_1\bar{x}$$ donde $\bar{y}$ y $\bar{x}$ son las respectivas medias muestrales.

En nuesto caso usando este set de datos generado ficticiamente obtenemos $\hat \beta_0$ = 2.000713 × 10⁴
y $\hat \beta_1$ = 299.82

Recordemos que esta es una estimación en base a los datos y no sabemos los verdaderos parámetros de la DGP real(proceso generador de datos). En este caso yo si lo sé porque generé los datos pero en la vida real es inaccesible.
Lo que hicimos fue estimar, a partir de un set de datos, ciertos coeficientes o característica de una población mucho más amplia. (Todos los trabajadores del país..)

Precisión de nuestros estimadores

En este caso teníamos una sola muestra pero podríamos haber tenido muchas muestras (K) de la misma población (muchos sets de datos con horas trabajadas y salarios).
Si estimáramos los coeficientes para cada uno de esos sets obtendríamos K pares de coeficientes, cada uno calculado con las particularidades de esos sets.
Se puede demostrar que el promedio de una cantidad grande de estimadores provenientes de muchas muestras se centra en el verdadero valor poblacional (si el modelo es correcto). Es decir que el promedio de los K $\hat \beta_1$ va a centrarse en el verdadero valor poblacional de $\beta_1$ ( y lo mismo para $\beta_0$).
Pero estos K parámetros centrados en el verdadero valor van a tener cierta dispoersión, es decir, pueden estar todos muy cerca del verdadero o estar muy dispersos pero que en promedio si quede centrado. Esto determina que tan preciso es el coeficiente que estimemos de una muestra. Este desvío estándar de los parámetros (SE) puede estimarse y depende de la varianza del error del modelo.

Puede ser útil para calcular los intervalos de confianza de los parámetros. Estos son intervalos que con X% de probabilidad contienen al verdadero valor del parámetro poblacional. Lo más habitual es calcular el intervalo de confianza al 95%. Para $\hat \beta_1$ esto es aproximadamente: $$ \hat \beta_1 \pm 2 \cdot SE(\hat \beta_1)$$ La interpretación sería que de 100 intervalos que construya de esta manera (de 100 muestras distintas), 95 van a tener al verdadero valor de $\beta_1$.

Por otra parte podemos realizar un test de hipotésis de los coeficientes. El más común es testear la siguiente hipótesis nula:
H_0 : No hay relación entre X e Y
contra la hipótesis alternativa
H_1 : Hay alguna relación entre X e Y
Lo cual se traduce en: $$ H_0 : \beta_1 = 0 $$ $$ H_1 : \beta_1 \neq 0$$

Lo que se hace es determinar si $\hat \beta_1$ está lo suficientemente lejos de 0 como para rechazar la hipótesis nula. Qué tan lejos es suficiente depende en gran parte del desvío estándar (SE) del coeficiente. Si el SE es grande , necesitaremos valores elevado de $\hat \beta_1$ para estar tranquilos con que el valor real no puede ser 0.
Para esto lo que se hace es calcular el estadístico t: $$ t = \frac{\hat \beta_1 - 0}{SE(\hat \beta_1)}$$ Que mide cuantos desvíos estándar $\hat \beta_1$ está alejado de 0. Si no hay relación entre X e Y se espera que el estadístico tenga una distribución t con n - 2 grados de libertad. Dado ese supuesto, lo que se hace es calcular la probabilidad de obtener un valor de t como el de nuestro estadístico, si este proviene de una distribución t con n-2 grados de libertad. Esta probabilidad se la conoce como p valor. Sería qué tan probable es encontrar un valor al menos tan grande como el de t si este proviniera de la distribución t con n-2 grados de libertad. Si esta probabilidad es muy chica (el umbral habitual es 0.05 pero depende del trabajo) uno rechaza la hipótesis nula en favor de la alternativa, suponiendo que sí hay una relación entre X e Y.

Precisión del modelo.

Naturalmente a uno le interesa saber qué tan bien ajusta nuestro modelo a los datos.

El método más habitual para regresión lineal es el R^2. Toma valores entre 0 y 1 porque es la proporción de la varianza de Y explicada por nuestro modelo.

\[ R^2 = \frac{TSS - RSS}{TSS} = 1 - \frac{RSS}{TSS}\]

donde $TSS = \sum (y_i -\bar{y})^2 $ es la suma de cuadrados totales y RSS es la suma de errores cuadrados que ya definimos antes. TSS mide la varianza total de Y y representa la variabilidad total inherent de la variable dependiente antes de correr la regresión. Por el contrario RSS mide la variabilidad que queda sin explicar por nuestro modelo (recuerden que proviene de los residuos). Por lo tanto el numerador TSS - RSS mide la parte de la variabilidad de Y que sí pudo ser explicada por el modelo, y lo divide por la variabilidad total. $R^2$ mide entonces la proporción de la variablidad que pudo ser explicada usando X. Cuanto más cerca de 1, mejor.

Regresión con Múltiples Predictores

Suena mucho más lógico tratar de explicar una variable dependiente no solo por una independiente si no por varias. La regresión lineal simple puede ampliarse a regresión lineal múltiple donde nuestro modelo pasa a ser: $$ \hat y_i = \hat \beta_0 + \hat \beta_1 x_i + ... + \hat \beta_p x_p$$ Y mantenemos un término de error con distribución normal y media 0.

En esencia la idea es la misma, explicar la variabilidad de Y basado en la variabilidad de nuestros predictores. La metodología para estimar los coeficientes suele ser Mínimos Cuadrados como vimos antes, sin embargo la solución no suele ser tan fácil de expresar y es más sencillo verlo en términos matriciales o simplemente ver los resultados desde el programa estadístico que estemos usando. No olvidar que varias regresiones simples no pueden sumarse en una resgresión múltiple, es decir, los coeficientes de las regresiones simples no tienen por que ser los mismos ni por qué mantener el signo cuando se juntan todas las variables en un solo modelo. Esto sucede porque la regresión múltiple estima coeficientes controlando por todas las otras variables, es decir, quitando el efecto de las otras. Por eso es que por separado quizás dos variables son significativas pero en una regresión múltiple solo una de ellas lo es. En general esto viene dado porque están correlacionadas y se comportan de manera similar. Al final del día no es fácilmente distinguible cuál es realmente la que lidera el efecto.

Hay relación entre la dependiente y los predictores? (test F)

En regresión simple vimos el test de hipótesis para ver si el coeficiente de $\hat \beta_1$ era significativamente distinto de 0. En regresión múltiple lo que debemos hacer es chequear si todos nuestros coeficientes son distintos de 0. (Y no uno por uno)

Lo cual se traduce en: $$ H_0 : \beta_1 = \beta_2 = ... = \beta_p = 0 $$ $$ H_1 : \text{Al menos algún } \beta_j \neq 0$$ El test de hipótesis se hace calculando el estadístico F. $$ F = \frac{(TSS - RSS) / p}{RSS / (n - p -1)}$$

Si los supuestos del modelo lineal se cumplen puede probarse que $E[RSS/(n - p-1)] = \sigma^2$ y si $H_0$ es verdadera $E[(TSS-RSS)/p] = \sigma^2$. Por lo tanto si no hay relación el estadístico F será cercano a 1 y si en realidad la hipótesis alternativa es verdadera el numerador será mayor que $\sigma^2$ y por lo tanto F será mayor que 1. Dependiendo de n, p y del nivel de significatividad que busquemos F deberá superar un umbral distinto para poder rechazar la hipótesis nula.
Es inevitable mirar los p-valores individuales sin embargo debemos tener cuidado particularmente cuando tenemos muchas variables. Por definición, algunos coeficientes saldrán significativos por azar aunque no tengan relación con la variable dependiente. En el caso típico de significativdad del 95%, esto sucede en promedio el 5% de las veces. Con muchas variables nuestras posibilidades de encontrarnos con al menos algún falso significativo aumentan notoriamente por lo que hay que mirar con cuidado. Por su parte, el estadístico F corrige en su cálculo por la cantidad de coeficientes y por lo tanto no se ve afectado por este problema.

Selección de variables

Cuando tenemos un set de datos grande es habitual tener que seleccionar cuáles son las variables importantes para el modelo. Más allá del conocimiento del problema (fundamental), idealmente lo mejor es probar una gran cantidad de modelos y con alguna métrica de comparación seleccionar los mejores. El problema es que la cantidad de modelos posible crece exponencialmente con la cantidad de variables y esto no es posible.
En el libro los autores mencionan como alternativas Forward Selection, Backward Selection y Selección mixta. Básicamente son enfoques que prueban una muestra de todos los modelos posibles según la significatividad de las variables. Son métodos iterativos. Habiendo avanzado la disciplina, llegado el caso buscaría otros métodos vigentes para atacar este problema.

Ajsute del modelo

Para verificar el ajuste del modelo se sigue usando el $R^2$ como métrica principal. En este caso es equivalente a la correlación al cuadrado de Y e $\hat Y$. Un punto a tener en cuenta es que el R^2 nunca puede disminuir al agregar variables ya que el peor escenario posible es que la nueva variable tenga coeficiente de 0 y el ajuste quede igual que antes. Lo que se hace para controlar por esto y poder comparar modelos es ajustar el $R^2$ por la cantidad de variables utilizadas o usar el RSE. De cualquier manera lo importante es recordar que el R2 sigue siendo útil en la regresión lineal múltiple.

Otras consideraciones

Las variables independientes admiten variables categóricas! (Binarias o multiclase). Ej: Educación máxima alcanzada. Hay que mirar con atención la interpretación. Alteran el intercepto según la categoría de la observación y puede alterar pendientes si se las incluye en interacción con alguna variable continua.
El modelo que venimos viendo el aditivo y lineal, pero podemos remover esos supuestos. Por ejemplo podemos agregar interacción entre variable y por lo tanto relajar la aditividad. Esto significa que las variables se modelan multiplicadas entre sí por ejemplo.

\[ Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_1X_2 + \epsilon \]

Podemos aproximar relaciones no lineales extendiendo el modelo a regresión polinómica. Suponiendo que los datos provienen de un modelo polinómico, podemos ver en el siguiente gráfico cómo cambia al agregar el término de no linealidad. La línea naranja es la regresión lineal simple y la celeste (que ajusta casi perfecto) es la regresión polinómica que respeta el proceso generador de los datos (la ecuación que se ve en el gráfico). Vemos que la variable Y depende de X linealmente pero también de X al cuadrado, lo que le da la curvatura.

Potenciales problemas

No linealidad de los datos

Si la relación entre nuestras variables independientes y la dependiente no es lineal nuestro modelo va a tener sesgo alto. Para el modelo simple es fácil de ver al graficar X vs Y en un gráfico de puntos pero con muchas variables eso ya no es tan sencillo. Un buen enfoque es realizar una regresión lineal y graficar los residuos contra los valores predichos.

Es un caso un poco extremo pero supongamos que la relación es no lineal, polinómica de orden 2 (como el ejemplo de arriba).
Si nosotros corremos una regresión lineal, nuestros residuos van a seguir un patrón muy obvio.

Claramente en ese gráfico el residuo no está centrado en 0... En cambio, si nosotros corremos una regresión lineal para un modelo realmente lineal, o para este caso, corremos un modelo no lineal, deberíamos ver una nube de puntos dispersa para los resiudos, centrada en 0 y con algún desvío estándar. Idealmente veríamos algo asi.

Correlación de los términos de error

Uno de los supuestos de la regresión lineal es que los errores no están correlacionados, es decir que el error $\epsilon_i$ de una observación no nos aporta información acerca del error $\epsilon_j$ de otra observación. Son independientes.
Si esto no se cumple lo que sucede es que el SE de los coeficientes estimados es menor al real y puede llevarnos a confiar más en nuestro modelo de lo que deberíamos.
Los errores correlacionados suelen suceder más frecuentemente en series de tiempo pero también pueden darse en estudios experimentales mal diseñados.

Varianza de los términos de error no constante. (Heterocedasticidad)

Otro supuesto de la regresión lineal es que la varianza de los errores es constante $Var(\epsilon_i) = \sigma^2$. Esto no siempre es el caso. En este ejemplo vemos como los resiudos siguen centrados en 0 pero con una dispersión mucho mayor a medida que avanzamos en el eje X. Entre las soluciones para este problema se encuentra transformar la variable dependiente - $ln(y)$ por ejemplo, o utilizar Mínimos cuadrados ponderados, que pondera por la inversa de la varianza. El libro no se explaya mucho más al respecto en este capítulo.

Outliers

Los outliers son observaciones cuya variable dependiente tienen valores que se alejan mucho del patrón regular de los datos, por ejemplo debido a un error de medición o problema al registrar la información. Los outliers pueden tener diversas consecuencias en los modelos lineales. Puede afectar la estimación de los parámetros, puede afectar el ajuste del modelo (caída del $R^2$) o puede por ejemplo aumentar los intervalos de confianza ya que el outlier afecta el RSE que es común a todos los intervalos. Todo esto puede ser generado por una sola observación. Generalmente si no se distinguen en el análisis exploratorio pueden saltar a la vista analizando los resiudos del modelo (o los residuos estandarizados).

Puntos con alto "leverage"

Si alguno tiene una traducción satisfactoria bienvenido sea.. Los puntos con alto leverage son aquellos cuyo valor de la variable independiente se aleja del rango estándar. Las observaciones con esta característica tienden a afectar en buena medida a la curva ajustada y por ende a los parámetros de nuestro modelo. Nuestra estimación por mínimos cuadrados puede verse muy influenciada por estos puntos e invalidar el ajuste por eso es muy necesario identificar estas observaciones.
En regresión simple es sencillo de ver porque resaltan si graficamos una nube de puntos pero en regresión múltiple es más difícil de ver ya que debemos encontrar anomalías en el conjunto de todas las variables. Es decir que una observación puede estar en el rango individual de cada variable pero si miramos a nivel conjunto, esa combinacion dentro de los rangos individuales es súper anómala. Con más de dos variables independientes se dificulta identificar visualmente. Para ayudar en estos casos se puede calcular el estádistico de leverage en algún programa estadístico.

Colinealidad

Este problema refiere a la alta correlación entre variables independientes del modelo, es decir que tienden a aumentar o decrecer de manera conjunta. Esto genera que sea difícil (o imposible en el extremo) diferenciar el impacto de cada una de ellas en la variable dependiente.
En una regresión lineal esto se traduce en aumento de la varianza de los estimadores y por ende incertidumbre sobre los parámetros estimados. A modo intuitivo, con variables con alta correlación puede haber una gran cantidad de combinaciones de coeficientes para estas variables que resulten en un mismo ajuste ($R^2$) y por ende mínimos cuadrados es indistinto frente a ellos. Cambiando alguna observación puede que el modelo pase de una combinación a otra muy disinta en ese arco de posibilidades. Otra consecuencia es que el aumento de la varianza de los coeficientes reduce el estadístico t que miramos para la significatividad y puede que lleve a no rechazar una hipótesis nula que debía ser rechazada. La potencia del test de hipótesis se ve disminuida por la colinealidad.
No solo la correlación sirve para detectar colinealidad. Puede existir multicolinealidad en donde varias variables son colineales aún sin tener alta correlación de a pares. Posiblemente se deba a combinación lineal generada por algunas de las variables. Para estos casos lo que se puede mirar es el VIF ( Variance Inflation Factor) en inglés. Este estadístico se calcula para cara variable y compara la varianza del estimador al tener la variable en el modelo versus ajustando un modelo solo con esa variable. Cuanto mayor es el VIF, mayores problemas de colinealidad resalta. Se puede calcular con la siguiente formula donde $R^2_{X_j|X_{-j}}$ es el $R^2$ de la regresión de $X_j$ contra todas las otras variables independientes del modelo. $$ VIF(\hat \beta_j) = \frac{1}{1 - R^2_{X_j|X_{-j}}}$$ La solución a este problema suele ser descartar alguna de las variables o agruparlas de alguna manera para quedarnos con una única variable que represente a ambas.

2019-05-06
in estadistica
8 min read

Aprendizaje Estadístico - ISLR Capitulo 2

Suponemos que las variables que encontramos en un set de datos son generadas a través de un proceso generador de datos (DGP por sus siglas en inglés) cuya expresión es: $$ Y = f(X) + \epsilon $$

Donde Y es la variable, en este caso la dependiente o la que queremos explicar. f(X) es una función respecto a otra/s variable/s X (independientes) y $\epsilon$ es el error irreducible, es decir un valor aleatorio con media 0 pero que no depende de otras variables, es al azar. Puede referir a errores de medición, cambios inmesurables en las situaciones del experimento o simplemente azar en la generación real de los datos. Cabe destacar que f(X) es desconocida para nosotros y justamente lo que queremos explorar con el análisis estadístico. Puede tenerse suposiciones o conocimiento de la forma funcional (lineal, no lineal, etc) pero en principio no tenemos mayores certezas y esperamos aprender a partir de la muestra que analizamos.

Por qué estimar f(X)?

Los dos principales motivos para interesarse en f(X) son predicción (de Y) e inferencia de los parámetros de f(X).

Predicción

Queremos predecir valores de Y para nuevos datos X. Como $\epsilon$ en promedio es 0 podemos aproximar Y de la forma: $$ \hat Y = \hat f(X)$$ La precisión de $\hat Y$ va a depender del error reducible y del error irreducible. El primero depende de qué tan bien nos aproximemos a la verdadera f(X) y puede ser potencialmente reducido si utilizamos las técnicas más adecuadas para el caso. El segundo error es justamente irreducible y es porque nuestra aproximación no puede tener en cuenta a $\epsilon$. El término aleatorio introducido por esa variable no lo podemos estimar para cada observación y por lo tanto debemos convivir con ese margen de error.

Suponiendo que tenemos una estimación $\hat f$ y un set de datos X puede probarse que: $$ E(Y - \hat Y)^2 = E[f(X) + \epsilon - \hat f(X)]^2$$ $$ = \underbrace{[f(X) - \hat f(X)]^2}\text{Reducible} + \underbrace{Var(\epsilon)}\text{Irreducible}$$

Donde $E(Y-\hat Y)^2$ es el promedio o valor esperado de la diferencia al cuadrado del valor real de Y y de la predicción correspondiente. $Var(\epsilon)$ es la varianza del término de error $\epsilon$.

Inferencia

Este enfoque se basa en entender la relación entre las variables de X y la dependiente. Es necesario entender bien la f(X) elegida para poder interpretar sus coeficientes y poder ver qué variables están asociadas con Y, cómo es esa relación, cuál es la forma de la función f(X), etc para poder actuar sobre las variables X o comprender su efecto aunque no siendo tan exigentes con el poder de predicción de nuestro modelo.

Como Estimar f(X)?

Métodos Paramétricos

Los métodos paramétricos se conforman por dos etapas.
La primera es asumir o suponer la forma funcional de f(X). Podes definir por ejemplo que f(X) es una función lineal de la forma $$ f(X) = \beta_0 + \beta_1X_1 + \beta_2X_2 + ... + \beta_pX_p$$

Una vez definida la forma del modelo la segunda etapa consiste en estimar los parámetros con algún método, a partir de los datos de entrenamiento. En este caso sería estimar todos los $\beta$. Por ejemplo para las funciones lineales se suele utilizar el método de mínimos cuadrados ordinario.
La ventaja de definir una forma funcional es que luego es más sencillo estimar sus parámetro y el problema se reduce a eso finalmente. Por el otro lado, posiblemente la forma que elijamos no sea exactamente igual a la real (DGP) y tengamos que aceptar que va a haber errores debido a eso. Si estamos muy lejos de la forma real esos errores serán groseros. Existen modelos flexibles que permiten ajustar modelos con diferentes formas de f(X) pero en general requieren estimar más parámetros y son más propensos a sufrir sobreajuste/overfitting que básicamente es ajustarse mucho al ruido o error ($\epsilon$) en los datos de entrenamiento y luego ajustar mal en testeo.

Métodos No Paramétricos

Los métodos no paramétricos no requieren definir explícitamente una forma funcional de f. Buscan un f que sea lo más cercano posible a los datos sin ser demasiado estricto o flexible. Al no asumir una forma puede cubrir potencialmente un rango mucho mayor. El problema es que al no reducir el problema a estimar parámetros necesitan muchas más observaciones para estimar f de forma medianamente precisa. En general uno tiene que decidir el nivel de "suavidad" del modelo, lo cual afecta que tan variable termina siendo la estimación. Sirve para encontrar el punto de fleixibilidad/rigidez del modelo que queremos para que no sobreajuste (ni falle demasiado).

Predicción vs Interpretabilidad

Uno puede elegir entre modelos flexibles o más rígidos. Con pocas observaciones a veces uno no puede alejarse mucho de los rígidos pero suponiendo que uno tiene muchos datos, a veces puede igualmente elegir rígidez frente a modelos flexibles que permiten ajustar varias formas de f. El motivo es que generalmente los modelos restrictivos son más fáciles de interpretar y se le puede dar un significado claro a sus coeficientes mientras que con formas muy flexibles no es sencillo entender el impacto de las variables de manera individual. La elección va a depender del objetivo del análisis y de qué tan bien o mal nuestros modelos ajustan a los datos.

Modelos Supervisados vs No Supervisados

Nuestros datos pueden tener una variable dependiente que queremos explicar o predecir en base a un set de variables independientes, con algún modelo a definir. Los casos de este estilo son llamados supervisados porque sabemos la "respuesta" (nuestra variable Y) y podemos validar nuestros modelos contra la realidad.
Si los datos no tienen una variable dependiente lo que se puede hacer es un análisis no supervisado donde por ejemplo lo que se puede hacer es agrupar las observaciones en clusters o grupos. Es decir segmentar en distintas clasificaciones y descubrir patrones. El desafío es que no hay en los datos nada contra qué validarlo, aunque sí contra el conocimiento del dominio o de la temática.

Regresión vs Clasificación

En los modelos supervisados, nuestra variable dependiente puede ser cuantitativa o cualitativa.
En el primer caso la variable toma valores númericos, como por ejemplo la altura de una persona, el precio de una propiedad, etc. Son problemas de regresión.
En el segundo caso la variable dependiente puede tomar el valor de una clase o categoría. Por ejemplo, género de una persona, si paga o no paga su deuda, etc. Son problemas de clasificación.

Midiendo el ajuste del modelo

Regresión

Para problemas de regresión una de las medidas más utilizadas es el Error Cuadrático Medio (MSE por sus siglas en inglés). $$ MSE = \frac{1}{n} \sum_{i=1}^n(y_i - \hat f(x_i))^2 $$

Es básicamente la diferencia promedio entre la realidad y lo que predice nuestro modelo elevado al cuadrado. Esto último es para que los errores sean siempre positivos aunque subestimemos o sobreestimemos (y por su comodidad para cálculos matemáticos).

En primero lugar se calcula este valor con los datos de entrenamiento sin embargo lo que realmente importa es como performa el modelo en datos de testeo, es decir en datos que no fueron utilizados para estimar f(X). Podemos decir que cada modelo debería tener un MSE de entrenamiento y un MSE de testeo. Debido a la posibilidad de sobreajuste y a las diferencias en muestras nada garantiza que el modelo que estimemos con menor MSE en entrenamiento también sea el de menor MSE en testeo.

A medida que aumentamos la flexibilidad de un modelo (sus grados de libertad) el MSE en entrenamiento va a disminuir, ya que tiene más herramientas para ajustarse a los datos pero puede que sobreajuste y por lo tanto no se traduzca en un menor MSE en testeo.

El tradeoff entre Sesgo y Varianza

No está demostrado en el libro pero es posible descomponer el MSE esperado de una observación de testeo en sesgo de $\hat f(x_0)$, varianza de $\hat f(x_0)$ y varianza del error irreducible $\epsilon$.

\[ E(y_0 - \hat f(x_0))^2 = Var(\hat f(x_0)) + [Sesgo(\hat f(x_0))]^2 + Var(\epsilon)\]

El lado izquierdo de la ecuación es el MSE esperado y coresponde MSE de testeo promedio que obtendríamos si estimaramos f utilizando una gran cantidad de sets de entrenamiento y testearamos cada uno en $x_0$.

Algunas observaciones:

El MSE nunca puede ser menor que la varianza de $\epsilon$. Es un término fijo y por eso se lo llama error irreducible.
La varianza es cuanto cambiar $\hat f$ si utilizamos otro set de entrenamiento. Siempre va a cambiar con otro set pero idealmente ese cambio no debería ser grande. Modelos muy flexibles tienden a cambiar más frente a distintos sets y son más inestables.
Sesgo es el error provocado por la diferencia entre el modelo elegido y el verdadero proceso generador de los datos. En general modelos más flexibles tienen menor sesgo ya que pueden ajustar mayor variedad de formas funcionales.
Al aumentar la flexibilidad de un modelo en general reducimos el sesgo pero aumentamos la varianza. En general en un primer momento el sesgo suele disminuir a mayor velocidad de lo que aumenta la varianza y por lo tanto el MSE esperado se reduce. Sin embargo llega un punto donde mayor flexibilidad reduce menos el sesgo que lo que aumenta la varianza y el MSE empieza a aumentar. Es el primer indicio de sobreajuste. Por eso se habla de tradeoff o "balance".
En la realidad donde la verdadera f del DGP es inobservable no suele ser posible calcular explícitamente el MSE de testeo, el sesgo o la varianza de un método estadístico pero el proceso de fondo aplica y siempre debemos tener en mente el tradeoff.

Clasificación

Para problemas de clasificación uno de los enfoques más frecuentes para cuantificar la precisión de una función estimada $\hat f$ se suele usar el porcentaje de error en los datos de entrenamiento.

$$ \frac{1}{n} \sum_{i=1}^nI(y_i \neq \hat y_i) $$ Básicamente es el porcentaje de observaciones clasificadas erroneamente. Al igual que con MSE es de gran importance el porcentaje de error en los datos de testeo.

Clasificador de Bayes

No lo demuestra en el libro pero la mejor manera de reducir el porcentaje de error en test es asignar a cada observación la clase con mayor probabilidad (según el DGP ). Es un concepto muy sencillo, dado X, asignar la clase cuya chance de acierto sea mayor.

$$ Pr(Y = j | X = x_0) $$ El porcentaje de error de Bayes (es decir el error luego de clasificar siguiendo esa regla) es análogo al Error Irreducible de regresión. Hay que tener en cuenta que la distribución condicional de Y dado X no lo sabemos en los casos aplicados en la vida real, sería como saber la función f(X) o el DGP y por lo tanto no lo podemos calcular.

K vecinos más cercanos (KNN en inglés)

Idealmente uno querría aplicar el clasificador de Bayes pero es imposible ya que no sabemos la distribución real de los datos (es justamente lo que queremos estimar). KNN intenta aproximarse a la distribución condicional para clasificar las observaciones. Lo que hace este método es, dado un valor de K que elegimos nosotros, clasificar cada nueva observación según la clase mayoritaria entre las K observaciones más cercanas a esta.

$$ Pr(Y = j | X = x_0) = \frac{1}{K} \sum_{i \in N_0} I(y_i = j) $$ El valor que seleccionemos de K afecta en gran medida las predicciones del modelo. Un K menor hace más variable el modelo ya que selecciona menos observaciones y por lo tanto pocos cambios en el set de entrenamiento cambian la clasificación. Suele reducir el sesgo pero ser mas variable. Es análogo a hacer más flexible un modelo en regresión. Valores de K más grandes seleccionan puntos en un entorno más abarcativo y por lo tanto suele ser más constante pero con sesgo superior.
Al igual que en regresión hay que tener cuidado con el sobreajuste. Reducir K garantiza menos errores en los datos de entrenamiento pero pasado un umbral la varianza aumenta en mayor medida y el porcentaje de error en test se incrementa.

Conclusión: Tanto en Clasificación como en Regresión la elección del nivel de flexibilidad es central en el éxito de método de aprendizaje estadístico.

2019-04-20
in matematica, algebra
7 min read

Esencia del Algebra Lineal

El álgebra lineal está por todas partes en estadística y data science. Matrices, vectores y transformaciones son términos que se escuchan seguido y están detrás de muchos de los métodos y algoritmos que se usan hoy por hoy. Aunque no sea necesario saber del tema para correr un modelo empaquetado en una librería de R, es muy útil entender lo que hacemos realmente ya que nos permite ver a los modelos como algo lógico y no una caja negra mágica.

Hay una serie de videos excelente en inglés que mediante visualizaciones y animaciones permite entender la intuición de muchos de los conceptos básicos, que solo con un libro puede ser medio críptico o poco imaginable. Para el que le interese: ESSENCE OF LINEAR ALGEBRA por 3Blue1Brown.

Este post, aunque quizás medio desordenado y sin mucha prolijidad, es una recopilación de algunas notas. Puede que queden términos en inglés intercalados.

Matrices y vectores

Vector vive en n dimensiones.
Suma de vectores es combinación lineal
En $R^{2}$ $\hat \imath = \left< 1, 0 \right> \text{y} \hat \jmath = \left< 0, 1 \right>$ forman una base. Cualquier punto es una combinación lineal de ellos.
Span es el espacio que pueden generar x vectores. $R^{2}$ es el span de $\left< 1, 0 \right> \left< 0, 1 \right>$
Vector puede ser pensado como una flecha desde el origen (0,0) a las coordenadas que lo identifican. O como un punto directo en las coordenadas..
Matriz es una transformación. Lleva un vector a otro punto. Si transformamos cada posible vector de un espacio por la matriz podemos ver como el espacio es transformado. Ej: rotar, invertir, estirar.
Si transformamos una base, cada punto nuevo puede generarse transformando la nueva base.
Por ej: $z = \left< 3, 2 \right> \text{es } 3\begin{bmatrix} 1 \\ 0 \end{bmatrix} + 2\begin{bmatrix} 0 \\ 1 \end{bmatrix} = 3\hat \imath + 2\hat \jmath$
Aplicando la transformación de la matriz A = $\begin{pmatrix} A & B \\ C & D \end{pmatrix} \text{obtenemos los nuevos vectores base } \hat \imath^{*} \text{y } \hat \jmath^{*}$
$z^{*} = 3\hat \imath^{*} + 2\hat \jmath^{*}$
Multiplicar 2 matrices es transformar un espacio con la primera matriz ( desde la derecha) y luego transformar el resultado por la segunda matriz. Ej: Rotar un espacio y luego invertirlo.
AB != BA -> El orden de las transformaciones importa y se lee de derecha a izquierda.
La matriz (transformación) ya dice como van a ser las nuevas bases.
Si la matriz es $\begin{pmatrix} A & B \\ C & D \end{pmatrix}$, el nuevo $\hat \imath^{*}$ es $\begin{bmatrix} A \\ C \end{bmatrix}$ y $\hat \jmath^{*}$ es $\begin{bmatrix} B \\ D \end{bmatrix}$
Ej: $z = \left< 3, 2 \right> z^{*} = \begin{bmatrix} A & B \\ C & D \end{bmatrix}\begin{bmatrix} 3 \\ 2 \end{bmatrix} = \begin{bmatrix} 3A + 2B \\ 3C + 2D \end{bmatrix}$
Se puede ver también como: $\(z = 3\hat \imath + 2\hat \jmath \text{ } z^{*} = 3\hat \imath^{*} + 2\hat \jmath^{*} = 3\begin{bmatrix} A \\ C \end{bmatrix} + 2\begin{bmatrix} B \\ D \end{bmatrix} = \begin{bmatrix} 3A + 2B \\ 3C + 2D \end{bmatrix}$\)
!!!. Las transformaciones afectan el area (en R2, el volumen en R3..) de las figuras en el espacio (todas por igual). El DETERMINANTE de una matriz mide ese cambio.
Si el determinante es 0 significa que se perdió una dimensión o que todo se comprimió. Pasa de $R^{2}$ a una recta (o a un punto!)
Si el determinante ** es < 0** significa que el espacio se invirtió (en sentido.. como dar vuelta una hoja) pero |DET| siguen siendo el cambio en el area.
A^-1A = I -> una transformación que no hace nada.
Si DET(A) = 0 no existe la matriz inversa. Ej. $R^{2}$ -> si det(A) = 0 la transformación lleva el espacio a una recta. No hay función que lleve cada vector de la recta a un punto en $R{2}$. No hay vuelta atrás.
Si una transformación lleva todos los puntos a una recta tiene rango 1, si lleva a un plano rango 2, y así.. RANGO es el número de dimensiones del output. Rango completo es cuando mantiene las dimensiones del input.
El conjunto de posibles outputs de $A\vec v$ es el Column Space = Span de las columnas
Cuando perdés dimensiones por la transformación todo un conjunto de vectores pasa a ser (0,0). Eso se llama Null Space o Kernel
Matrices no cuadradas cambian la dimensión del espacio.
$$ \begin{bmatrix} A & D \ B & E \ C & F \end{bmatrix} \begin{bmatrix} 1 \ 1 \end{bmatrix} = \begin{bmatrix} A + D \ B + E \ C + F \end{bmatrix} $$ Quedan todos los puntos de $R^{2}$ en un plano en el espacio $R^{3}$. De acá viene la restricción para multiplicar matrices. La cantidad de columnas de la transformación tiene que ser igual a la dimensión del input

DOT PRODUCT o PRODUCTO INTERNO

Dot product entre dos vectores equivale a proyectar uno en el otro y multiplicar sus largos. $\vec A \cdot \vec B = |A^{*}| * |B|$
$A^{*}$ es el vector A proyectado en B.
$\vec B$ es un vector 2D pero también se lo puede ver como una matriz 1x2 que lleva del 2D a la recta.
$\vec B \cdot \vec A = B \vec A \text{que sería llevar A al espacio transformado por B.}$
$B = \begin{bmatrix} B_x & B_y \end{bmatrix}$ tiene en sus columnas donde queda $\hat \imath \text{y } \hat \jmath$ (los vectores unitarios) al ser transformados o algun valor escalado de esto.
$\vec A \cdot \vec B$ es el valor de A en la recta a la que te lleva la transformación B.
Es equivalente proyecto B en A.
Si Dot Product > 0, tienen dirección similar.
Si Dot Product = 0, son ortogonales - proyección que cae en el origen.
Si Dot Product < 0, tienen direcciones opuestas.

CROSS PRODUCT

Está definido para vectores en $R^{3}$
El cross product $\vec u \times \vec v$ es el area del paralelograma que se puede imaginar con las paralelas de los vectores (imaginandolo en $R^{2}$. El signo depende de la orientación de los vectores. El vector de la "derecha" tiene que estar primero para que el cross product sea > 0.
En realidad el paralelograma formado por dos vectores en R³ tiene area equivalente al Largo del vector output de su cross product. Ese vector es ortogonal al paralelograma.

CAMBIO DE BASE

Distintos sistemas de coordenadas definen $\hat \imath = \left< 1, 0 \right>, \hat \jmath \left< 0, 1 \right>$ como algo distinto. NO hay una sola "grilla" válida. El espacio no tiene grilla predeterminada.
Un mismo vector tiene distintas coordenadas según el sistema de bases desde donde se lo mire.
Para pasar de una base a otra se aplica una transformación lineal.
Si $\vec v$ es un vector que queremos pasar de una base a otra, lo transformamos por la nueva base.
Y $\hat \imath^{*} = \left< \hat \imath^{*}_1, \hat \imath^{*}_2 \right>, \hat \jmath^{*} = \left< \hat \jmath^{*}_1, \hat \jmath^{*}_2 \right>$ Entonces: $\(\begin{bmatrix} \hat \imath^{*}_1 & \hat \jmath^{*}_1 \\ \hat \imath^{*}_2 & \hat \jmath^{*}_2 \end{bmatrix} \begin{bmatrix} v_1 \\ v_2 \end{bmatrix} = \begin{bmatrix} v^{*}_1 \\ v^{*}_2 \end{bmatrix}$\)
Donde $\begin{bmatrix} v^{*}_1 \\ v^{*}_2 \end{bmatrix}$ es el vector en la nueva base, es decir, serían las coordenadas del vector $\vec v$ en el nuevo sistema de coordenadas y representando ese punto bajo el sistema de coordenadas original. -> Como se vería $\vec v$ en la nueva base? desde un punto de vista de la base original.
La matriz transforma un vector siguiendo en el lenguaje de la base original.
Ej: Si $\vec v$ es (1,2) en el sistema cartesiano típico y aplicamos la matriz de cambio de base, un vector (1,2) bajo otros ejes se ubicaría en otro punto del espacio. Qué punto es ese bajo el sistema cartesiano? Es (1,2) en el nuevo, pero queremos saber su equivalente en el sistema original.
Por otra parte si queremos saber que coordenadas tomaría el vector $\vec v$ bajo otra base debemos multiplicar por la inversa de la transformación. Transforma el vector al lenguaje de la nueva base. Responde a la pregunta. Qué coordenadas toma el punto V_1, V_2 del espacio en el sistema nuevo?
Para aplicar una transformación a otra base conviene llevar el vector a transformar a la base original, transformar y reconvertir a la nueva base. $$ [A]^{-1}[T][A]\vec v = \vec v^{*}$$
A lo expresa en términos de la base original, luego se le aplica la transformación T y luego se lo devuelve al lenguaje de la nueva transformación.

Eigenvalues y Eigenvectors (autovalores y autovectores)

!!! Al aplicar una transformación lineal a un espacio algunos vectores no cambian de dirección, solo se estiran o contraen pero sobre la misma recta. El resto sí se mueve. Los que se mantienen son los eigenvectors, y su factor de expansión o contracción es su eigenvalue.
Si A es la matriz de transformación, $\vec v$ es un eigenvector y $\lambda$ su eigenvalue. $$ A\vec v = \lambda \vec v$$
Si una transformación es una matriz diagonal, lo único que hace es estirar $\hat \imath \text{y } \hat \jmath$ por lo tanto los vectores base son eigenvectors y la diagonal son los eigenvalues.
Si cambiamos la base a una formada por los eigenvectors (que spanean el espacio) de la matriz podemos aplicar la transformación (la matriz original de donde salieron los eigenvectors) a esta nueva base y solo la va a estirar, por lo tanto es una transformación diagonal, que permite calculos mucho más fácil. Después habría que volver a la base original.
A -> Matriz de transformación E -> Matriz de autovectores que forman la nueva base $\begin{bmatrix} e_11 & e_21 \\ e_12 & e_22 \end{bmatrix}$ D -> Matriz Diagonal cuyos valores son los eigenvalues $$ E^{-1}AE=D$$

E cambia la base a eigenvectors (expresado en la base original), A aplica transformación y E^-1 lo lleva al lenguaje de la nueva base (queda expresado en las nuevas coordenadas)

Espacios Vectoriales Abstractos

!!! Ver funciones como un tipo especial de vectores.
Las funciones se pueden sumar y escalar $f(x) + g(x) \text{y } 2f(x)$
Existen transformaciones lineales de funciones, convierten una función en otra. También conocidas como "operadores"
Para que una transformación sea lineal tiene que cumplir aditividad y mulitplicación por escalar
$$ L(\vec v + \vec w) = L(\vec v) + L(\vec w)$$ $$ L(c\vec v) = cL(\vec v)$$
En general cualquier espacio que cumpla los axiomas los espacios vectoriales puede ser considerado uno.

2019-04-03
in estadistica, R
4 min read

Funciones de Probabilidad y Distribucion

Variables Aleatorias

Consideremos un experimento cuyo espacio muestral denominaremos S.
Una funcion valuada en el dominio de los reales definida en S es una variable aleatoria.

En otras palabras es una función que asigna a cada resultado posible de un experimento un valor real.

Por ejemplo:

Si el experimento es lanzar una moneda 10 veces hay 2¹⁰ combinaciones posibles de caras (o) y cruz (x).
Si definimos la variable aleatoria X como cantidad de caras entonces X(s) será la cantidad de caras del experimento.
Si s resulta ser la secuencia ooxxxoxxxo entonces X(s) = 4.

Distribucion de una variable aleatoria

Si tenemos la distribución de probabilidad del espacio muestral del experimento podemos determinar la distribución de probabilidad de cualquier variable aleatoria válida.

Volviendo al ejemplo de la moneda. Dijimos que hay 2¹⁰ combinaciones de cara o cruz. La cantidad de combinaciones de X caras en 10 lanzamientos es $P(X = x) = \binom{n}{x} \frac{1}{2^{10}}$ para $x = 0,1,2,..,10$

Distribuciones Discretas

Una variable aleatoria tiene una distribución discreta si solo puede tomar valores de una secuencia (generalmente finita pero puede no serlo).

La función de probabilidad le otorga una probabilidad puntual a cada valor de esa secuencia.
Los valores por fuera de la secuencia tienen probabilidad = 0
La suma de todas las probabilidades tiene que ser 1

Distribución Uniforme

En el caso de la dsitribución uniforme, supongamos que la variable puede tomar valores de 1 a k. La función de probabilidad será $f(x) = \frac{1}{k}$ para x = 1,2,...,k. Y 0 para todos los otros valores.

si k = 10
Los valores de la variable serán cualquier entero entre 1 y 10
Cada valor tendrá probabilidad $\frac{1}{10}$

Distribución Binomial

En el caso de la dsitribución binomial se asumen dos posibles resultados, uno con probabilidad p y su contraparte con probabilidad 1-p.
Por ejemplo la probabilidad p de que una máquina genere un producto defectuoso y 1-p de que sea no defectuoso.
Si una máquina produce n productos va a generar X productos defectuosos. La variable aleatoria X tendrá una distribución discreta y sus posibles valores irán de 0 a n.
Para cualquier valor de x (entre 0 y n), la probabilidad de que la máquina genere x productos defectuosos entre los n producidos (de una secuencia particular) es $p^{x}q^{(n-x)}$
Como existen $\binom{n}{x}$ distintas secuencias posibles con x defectuosos entre los n productos tenemos que:
$Pr(X = x) = \binom{n}{x}p^{x}q^{(n-x)}$
La función de probabilidad será $f(x) = \binom{n}{x}p^{x}q^{(n-x)}$ para x = 0,1,2,...,n. Y 0 para todos los otros valores.

Para usar esta distribución en R tenemos los siguientes comandos:

Para generar n escenarios al azar donde se producen size productos con probabilidad p de ser defectuosos. El resultado es la variable x por escenario. Es decir la cantidad de defectuosos.
En el primer escenario x = 0, en el segundo x = 1 y así.

set.seed(1)
rbinom(n = 10, size = 5, p = 0.2 )

##  [1] 0 1 1 2 0 2 3 1 1 0

# random binomial

Para saber la probabilidad de obtener x productos defectuosos si una máquina produce size productos y la probabilidad de que produzca un defectuoso es prob.
Hay probabilidad de 0.0264 de obtener 5 defectuosos si producimos 10 con probabilidad 0.2.

dbinom(x = 5, size = 10, prob = 0.2)

## [1] 0.02642412

Para saber la probabilidad acumulada de obtener q o menos productos defectuosos si la máquina fabrica size objetos, con probabilidad de defecto prob. Hay probabiliad de 0.879 de obtener 3 o menos defectuosos si la máquina produce 10 objetos con probabilidad 0.2 de defecto. Es decir, es la suma de obtener exactamanete 0 defectuosos, más exactamente 1 defectuoso, más exactamnente 2 defectuosos, más exactamente 3 defectuosos.

pbinom(q = 3, size = 10, prob = 0.2)

## [1] 0.8791261

Distribuciones Continuas

Una variable aleatoria X tiene una distribución continua si existe una función f definida en los reales tal que para algún intervalo A
$Pr(X \in A) = \int_{A} f(x)$

La función f es la función de densidad de probabilidad. PDF por sus siglas en inglés.
La probabilidad de que X tome algún valor en un intervalo se encuentra integrando f en ese rango.

Por ejemplo para la distribución uniforme en un intervalo (a,b) podemos ver que su pdf (o función de densidad de probabilidad) es
$f(x) = \begin{cases}\frac{1}{b-a} & \text{para } a \leq x \leq b \\ 0 & \text{resto}\\ \end{cases}$

Distribución Normal

Para la distribución Normal tenemos los siguientes comandos:

Para obtener n variables aleatorias provenientes de una normal con media mean y desvío sd

set.seed(1)
rnorm(n = 5, mean = 10, sd = 2)

## [1]  8.747092 10.367287  8.328743 13.190562 10.659016

Para obtener el valor de la pdf de la normal para algún valor de X en particular. Recuerden que no es una probabilidad, solo es el valor de la función. Las probabilidad se encuentra integrando la función en el intervalo deseado.
Si graficáramos los valores de dnorm para el intervalo -3,3 obtendríamos la forma típica de la normal.

dnorm(0.5, mean = 0, sd = 1)

## [1] 0.3520653

Para obtener la probabilidad acumulada hasta determinado punto. También conocido como Función de Distribución o Función de Distribución Acumulada C.D.F. por sus siglas en ingles Por ejemplo, cual es la probabilidad de obtener un valor igual o menos a 1.5 si tomamos una muestra de una normal estándar

$N \sim (0,1)$

pnorm(q = 1.5, mean = 0, sd = 1)

## [1] 0.9331928

Hay 93.31% de chances de obtener un valor inferior a 1.5 si tomamos una muestra al azar de una normal con media 0 y desvío 1.

La inversa también se puede calcular facilmente en R. Que valor debe tomar la variable aleatoria normal si deseo tenes un 93.31% de chances de obtener un valor menor o igual a ese?

qnorm(p = 0.9331, mean = 0, sd = 1)

## [1] 1.499284

La diferencia respecto al código anterior es porque redondeamos la probabilidad.

2019-03-30
in R
6 min read

Introduccion a graficos con mapas

Data

Vamos a ver un ejemplo sencillo para representar información visualmente sobre mapas. En este caso un pequeño dataset de incendios forestales en Argentina de 2012 a 2015. La idea es usar ggplot y mantener el enfoque de gráficos por capas.

Vamos a necesitar.

> tidyverse
> rgdal
> rgeos

Tendremos como input las provincias, departamento, cantidad de focos por incendio, area afectada y año de inicio y fin. Cada observación es un incendio.

Para este ejemplo nos vamos a centrar en las provincias, los focos y su efecto sin importar la fecha.

Empezamos cargando la data.

library(tidyverse)
# Load Raw data
raw <- read.csv("../../static/post/2019-03-30-introduccion-a-graficos-en-mapas/focosincendio.csv", sep = ";")
raw <- as.tibble(raw)

## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## i Please use `as_tibble()` instead.
## i The signature and semantics have changed, see `?as_tibble`.

glimpse(raw)

## Rows: 120
## Columns: 11
## $ pais_id         <int> 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32,~
## $ pais            <chr> "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina~
## $ provincia_id    <int> 6, 14, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 30, 30, 30, 30, 30, 30, 30, 30, 30, ~
## $ provincia       <chr> "Buenos Aires", "Córdoba", "Corrientes", "Corrientes", "Corrientes", "Corrientes", "Corrientes", "Corrientes", "Corrientes", "C~
## $ departamento_id <int> 833, 14, 56, 28, 28, 28, 70, 84, 84, 84, 112, 112, 119, 119, 126, 147, 154, 154, 161, 168, 168, 168, 168, 8, 15, 15, 15, 15, 28~
## $ departamento    <chr> "Tres Arroyos", "Calamuchita", "General Alvear", "Concepción", "Concepción", "Concepción", "Goya", "Itizaingó", "Ituzaingo", "I~
## $ sup_afectada    <dbl> 2400.00, 50.00, 257.00, 130.00, 5.00, 146.00, 30.00, 294.30, 378.00, 158.00, 300.00, 450.00, 450.00, 15.00, 20.00, 141.00, 295.~
## $ uni_med_id      <chr> "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "ha", "~
## $ cant_focos      <int> 1, 1, 1, 2, 1, 1, 1, 3, 1, 3, 1, 1, 5, 2, 2, 1, 3, 1, 1, 18, 2, 3, 2, 1, 0, 0, 1, 7, 0, 2, 1, 0, 0, 0, 5, 10, 10, 3, 1, 1, 1, 1~
## $ año_inicial     <int> 2014, 2015, 2012, 2012, 2013, 2015, 2012, 2012, 2013, 2014, 2012, 2013, 2012, 2013, 2015, 2012, 2012, 2014, 2012, 2012, 2013, 2~
## $ año_final       <int> 2014, 2015, 2012, 2012, 2013, 2015, 2012, 2012, 2013, 2014, 2012, 2013, 2012, 2013, 2015, 2012, 2012, 2014, 2012, 2012, 2013, 2~

Exploramos un poco el dataset. Lo que nos vas a interesar represetnar es la segunda parte del código. Variables agregadas a nivel provincia.

# Generate Summary to Explore
sum_year <- raw %>% group_by(año_inicial) %>%
  summarise(focos = sum(cant_focos), sup = sum(sup_afectada, na.rm = TRUE)) %>%
  mutate(sup_prom = sup/focos)

# Actual data to be plotted
sum_prov <- raw %>% group_by(provincia) %>%
  summarise(focos = sum(cant_focos), sup = sum(sup_afectada, na.rm = TRUE)) %>%
  mutate(sup_prom = sup/focos) %>%
  arrange(desc(sup_prom)) %>%
  mutate(provincia = as.character(provincia))

head(sum_prov)

## # A tibble: 6 x 4
##   provincia    focos   sup sup_prom
##   <chr>        <int> <dbl>    <dbl>
## 1 Buenos Aires     1 2400    2400  
## 2 Corrientes      55 5440.     98.9
## 3 Misiones        10  679      67.9
## 4 Córdoba          1   50      50  
## 5 Entre Ríos      12  594      49.5
## 6 Jujuy           62 2721.     43.9

Tenemos por provincia la cantidad de focos y la superficie afectada. Vamos a usar la superficie promedio por foco para visualizar la magnitud de cada foco de incendio.

Ahora llega lo importante. ¿Cómo representar esta data provincial en un mapa de Argentina?

Mapa

Necesitamos un ShapeFile de Argentina, que basicamente es el tipo de archivo que se usa para representar mapas en gráficos. Contiene divisiones del país (provincias) con sus respectivas coordenadas y nombres.

Utilizaremos data descargada del siguiente link. url <- "http://biogeo.ucdavis.edu/data/diva/adm/ARG_adm.zip"

Lo descargué y deszipee en la computadora. Lo leemos con una librería particular RGDAL.

dsn contiene la ruta a la carpeta con los archivos del shapefile. Layer apunta a al set de archivos que contiene la data que queremos. Generalmente hay otros sets con información no relevante al gráfico.

argentina <- rgdal::readOGR(dsn = "../../static/post/2019-03-30-introduccion-a-graficos-en-mapas/ARG_adm", layer = "ARG_adm1", use_iconv=TRUE, encoding='UTF-8')

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\DataScience\StatsBlog\blogStats\static\post\2019-03-30-introduccion-a-graficos-en-mapas\ARG_adm", layer: "ARG_adm1"
## with 24 features
## It has 9 fields
## Integer64 fields read as strings:  ID_0 ID_1

Es un archivo S4 por lo que se utiliza "@" para acceder a su contenido. Por ejemplo:

head(argentina@data)

##   ID_0 ISO    NAME_0 ID_1                 NAME_1           TYPE_1        ENGTYPE_1 NL_NAME_1
## 0   12 ARG Argentina    1           Buenos Aires        Provincia         Province      <NA>
## 1   12 ARG Argentina    2                Córdoba        Provincia         Province      <NA>
## 2   12 ARG Argentina    3              Catamarca        Provincia         Province      <NA>
## 3   12 ARG Argentina    4                  Chaco        Provincia         Province      <NA>
## 4   12 ARG Argentina    5                 Chubut        Provincia         Province      <NA>
## 5   12 ARG Argentina    6 Ciudad de Buenos Aires Distrito Federal Federal District      <NA>
##                                                                               VARNAME_1
## 0                                                                   Baires|Buenos Ayres
## 1                                                                               Cordova
## 2                                                                                  <NA>
## 3                                                        El Chaco|Presidente Juan Peron
## 4                                                                                  <NA>
## 5 BUENOS AIRES D.F.|Capital Federal|Distretto Federale|Distrito Federal|Federal Capital

Este tipo de archivos tiene una estructura complicada y hay varias librerías útiles. Con el fin de mantenernos dentro del tidyverse usaremos el enfoque de ggplot por capaz para graficar. Primero necesitamos llevar la información del shapefile a un dataframe.

# Transformo a dataframe.
argentina_df <- broom::tidy(argentina)

## Regions defined for each Polygons

# id es la provincia
head(argentina_df)

## # A tibble: 6 x 7
##    long   lat order hole  piece group id   
##   <dbl> <dbl> <int> <lgl> <fct> <fct> <chr>
## 1 -60.2 -33.3     1 FALSE 1     0.1   0    
## 2 -60.2 -33.3     2 FALSE 1     0.1   0    
## 3 -60.2 -33.3     3 FALSE 1     0.1   0    
## 4 -60.2 -33.3     4 FALSE 1     0.1   0    
## 5 -60.2 -33.3     5 FALSE 1     0.1   0    
## 6 -60.2 -33.3     6 FALSE 1     0.1   0

Ahora genero un diccionario de ids con su respectiva provincia para poder linkear mi mapa con la data de incendios.

ids <- cbind.data.frame(provincia = as.character(argentina@data$NAME_1), id = as.character(rownames(argentina@data))) %>%
  mutate(provincia = as.character(provincia), id = as.character(id))

El dataframe generado a partir del shapefile tiene las coordenadas que forman cada provincia y el id. Ahora lo que haremos es pegarle la data de incendios para poder utilizarla sobre el mapa. Está generado el log de los focos porque en una prueba intenté usar la data transformada para suavizar outliers pero no se verá en esta versión.

# Agrego provincia por ID y data de incendios a la data del shapefile.
argentina_df2 <- argentina_df %>% left_join(ids, by = "id") %>%
  left_join(sum_prov, by = "provincia") %>%
  mutate(focos = ifelse(is.na(focos),0.5,focos),
         logfocos = log(focos), 
         logfocos2 = logfocos - min(logfocos))

Por ahora vamos a poder ser capaces de graficar un mapa de argentina con sus provincias delimitadas y pintarlas según la cantidad de focos de incendio por ejemplo. Pero además vamos a querer agregar alguna forma sobre cada provincia. Por ejemplo un punto de distinto tamaño según el area afectada por cada foco en promedio. Para eso necesitamos localizar el centroide de cada poligono, es decir, el centro de cada provincia. Para ellos usamos la librería RGEOS Luego a cada centro le agrego la data que voy a querer usar. En este caso la superficie promedio afectada por foco de cada provincia.

# Calculo el centro de cada poligono (provincias)
# para obtener el "centro" donde iran los puntos o nombres.
centros <- rgeos::gCentroid(argentina, byid = TRUE) %>%
  as.data.frame() %>%
  mutate(id = rownames(.))

# Agrego data relevante para el ploteo (superficie promedio)
centros2 <- centros %>% left_join(ids, by = "id") %>%
  left_join(sum_prov, by = "provincia") %>%
  mutate(focos = ifelse(is.na(focos),0.5,focos),
         sup = ifelse(is.na(sup),0,sup),
         log_sup_prom = log(sup_prom),
         sup_prom_sin_outlier = ifelse(sup_prom > 150, 150,sup_prom) ) # esto es para suavizar el outlier. Buscar otro enfoque.

Grafico!

Ya tenemos todo. Tenemos el mapa, tenemos la cantidad de focos de incendio por provincia y tenemos el centro de cada provincia donde vamos a incluir un punto que muestra la intensidad de los incendios. Simplemente graficamos sigueindo la lógica por capas de ggplot. Las provincias sin puntos son aquellas que no tuvieron ningún incendio.

ggplot() +
  geom_polygon(data = argentina_df2, aes(x=long, y = lat, group = group,  fill = focos), color = "white") + # mapa de argentina
  # coloreado segun cantidad de focos
  coord_fixed(0.8) + # tamañp del mapa
  scale_fill_gradient2("Cantidad de Focos de incendio", low = "white", mid = "lightgreen", high = "darkred") + # escala de colores para focos
  geom_point(data = centros2, aes(x = x, y = y, size = sup_prom_sin_outlier)) + # puntos por provincia con superficie promedio
  scale_size(name = "Superficie Promedio Afectada (ha)",range = c(1,5)) + # escala de los puntos
  guides(fill = guide_legend(order = 1), # Orden de los leyendas a la derecha.
         size = guide_legend(order = 2)) # Por algun motivo esto discretizo la leyenda de focos

2018-11-10
in R
5 min read

Como correr un proyecto con R en Google Cloud

Hay situaciones en que nuestras computadoras no alcanzan para correr ciertos algoritmos por la cantidad de memoria o núcleos que tenemos (y comprar otra no es una opción..). Una solución es correr nuestro proyecto en "la nube", es decir en servidores ajenos mantenidos por empresas. Los servicios de esta índole más conocidos son:

Google Cloud
Amazon AWS
Microsoft Azure

En este post usaremos el primero. Como es de esperar, estos servicios son pagos y si su negocio lo amerita son una gran opción. Igualmente Google Cloud ofrece U$S300 de regalo al crear una cuenta por lo que podrían hacer uso para algún proyecto o pruebas. Les aseguro que no es particularmente bajo el monto. Solo tienen que registrarse y asociar una tarjeta de crédito y no abonar nada.

Google Cloud ofrece un montón de servicios y opciones de las cuales presentaremos lo más básico pero igualmente suficiente para correr un xgboost en gigas y gigas de datos con cientos de variables je.

Crear una máquina virtual instalando R
Crear un Bucket que sirve como Disco duro para guardar data, outputs, etc
Correr un algoritmo

Máquina Virtual

Una vez registrados, lo primero que vamos a hacer es crear la máquina virtual e instalar R y ciertos paquetes. Lamentablemente el proceso es bastante engorroso para quienes no conocen bash ni están familiarizados con Cloud. Es todo por consola y para nada intuitivo sin leer la documentación, que está en inglés. Es un proceso largo y la documentación más clara que encontré se encuentra en este BLOG. No tiene sentido intentar decir lo mismo que él pero peor. Les recomiendo seguirlo y van a obtener una máquina virtual con Rstudio instalado y ciertas dependencias útiles para la mayoría de las librerías que se usan.

Crear Bucket

Ir a Clickear los 3 puntos a la derecha y entrar a "Create key"

Se les descargará un archivo .json que deben renombrar a ** privatekey_inicial.json ** . Luego ir a:

Crear un bucket cuyo nombre no tenga espacios ni caracteres raros. Va a pedirles que sea un nombre que no esté siendo usado por nadie más. Hecho esto el bucket está creado y se le pueden subir archivos. Luego desde los scripts también se le van a poder escribir directamente, es decir, guardar los outputs. Ahora vamos a cambiar los permisos del bucket. Copiarse el account id desde aquí.

Luego ir a storage y en los 3 puntitos de nuesto bucket clickear "Edit Bucket permissions" y pegar el ID.

Crear Imagen

Lo que haremos ahora es crear una imagen de la maquina virtual que generamos (donde instalamos R) para poder levantar futuras máquinas y que directamente tengan instalado R y los paquetes si lo deseamos (invocando a esta imagen). Dejo imágenes de un tutorial que van a ser más claras que yo.

Luego ir a Compute Engine → Images

Status

Resumiendo brevemente:

Tenemos una cuenta y U$S 300 disponibles.
Creamos una VM con R y paquetes.
Generamos una imagen de esa VM para poder retutilizarla.
Creamos un bucket donde almacenar inputs y outputs, linkeado a la imagen.

Estamos casi listos. Falta correr algun script!

Cómo correr un script de R.

Una vez seteado todo lo anterior (y asegurándonos de haber apagado la maquina virtual utilizada) lo que debemos hacer es crear otras VM (instance) con los núcleos y Ram que cremos convenientes - esto depende totalmente de la complejidad del algoritmo que vayan a correr-, asignando una región donde preferentemente sea de noche o fin de semana para que esté menos saturada. Lo más importante es cambiar el boot disk y seleccionar dentro de "custom" la imagen que hayamos generado con R. De esta manera la VM que iniciemos ya tendrá R, Rstudio y librerías instaladas.

Justo debajo de lo que se ve en la imagen hay opciones extras de "Management, security, disks", etc. Depende la importancia de lo que estén haciendo puede ser bueno asegurarse que la opción Preemptibility esté OFF. Si está ON, la VM se correrá en un servidor que máximo puede durar 24hs (posterior a eso se apaga automáticamente la VM) y más importante aún, están sujetos a disponibilidad de Google, es decir que si hay mucha demanda de servidores pueden apagarles el suyo sin consultar. Lo positivo es que son mucho más baratas. Al tenerlo OFF, se aseguran que su VM estará encendida durante todo lo que tarde el script en correr y no va a depender de la demanda. Queda a criterio de cada uno.

Ya estamos listos para correr. Hay dos maneras sencillas:

Directo desde Rstudio en la VM.

Haciendo click en el IP de su instancia. Se abrirá Rstudio y deberán poner Usuario y Contraseña seteados durante la instalación. Ahi pueden trabajar como si fuera directamente R. Recomiendo tener Script listo porque no es muy dinámico trabajar en vivo ahí.

Desde la terminal, llamando a script en el bucket.

Otra opción es subir su script al bucket que generaron y en Rstudio (como en 1.) ir a la terminal y correr:

Rscipt --vanilla ~/cloud/cloud1/pathToScript/script.r

O desde la consola (que basicamente system() simular ser la temrinal) :

system(Rscipt --vanilla ~/cloud/cloud1/pathToScript/script.r)

Desde la terminal de la instancia.

Y ahí deberían poder correr sin problema.

Rscipt --vanilla ~/cloud/cloud1/pathToScript/script.r

En cualquiera de los 3 casos presten atención a las rutas que usan en sus scripts para referenciar al bucket. Si linkearon el bucket a la imagen de la VM como vimos en el post deberían poder usar la siguiente ruta genérica.

"~/cloud/cloud1/RestoDelPath/"

Y eso es todo. Ya pueden levantar una VM en google cloud y correr algoritmos con gigas y gigas de data sin quemar su PC!

2018-11-03
in R
4 min read

Organización de Proyectos

Este va a ser un breve post sobre cómo organizar los proyectos que hagan en R. Es al día de hoy la que utilizo y si googlean van a ver que existe y es usada en el ámbito. Puede ser esta tal cual o alguna alternativa similar.

El enfoque es muy sencillo e intuitivo. La idea es tener por separado cada componente del proyecto y de manera clara y segmentada para poder acceder rápidamente a lo necesario, ya sea código, datos, outputs, etc.

Primero y principal. CREEN UN PROYECTO en R. Esto va a facilitar todo el manejo de rutas, llamados a otros códigos y mismo para compartir si es necesario con otra gente ya que las rutas que utilicemos serán relativas a la ubicación en Disco del proyecto. Por lo tanto si el proyecto está creado en C:\ProyectoR y en algún codigo llamamos a read.csv("/Datos/datos.csv") esto funcionará en cualquier computadora donde el proyecto tenga en su directorio una carpeta "Datos", independientemente de la ruta donde se encuentre. Yo puedo mover toda la carpeta del proyecto a C:\OtraRuta y ejecutar ahí el read.csv sin tener que actualizar la ruta. Es uno de los problemas más básicos y molestos al trabajar con códigos ajenos o mover nuestros proyectos de lugar.

Carpetas

Una vez creado el proyecto, lo que sugerimos es crear una estructura de carpetas como la que se ve en la imagen siguiente.

DATA

Contiene la data que será el input de nuestros proyecto. A su vez podemos guardar archivos intermedios que hayamos ido procesando. Pueden adaptarlo como prefieran pero sugerimos guardar en "raw",la data que será input del análisis, ya sean csvs, txts, htmls, etc. En "working"" ir guardando objetos importantes o que lleven tiempo de procesar asi se pueden leer directamente en vez de tener que correr el código nuevamente en una próxima sesión. Para ellos se usa el comando saveRDS(). En "final" guardar los objetos finales del análisis.

DOCS

Acá guardamos archivos auxiliares útiles como diccionarios de variables, links a webs, consignas, documentación, etc.

OUTPUT

Acá exportamos los resultados del análisis, desde gráficos que vayamos a usar en el reporte, el informe final que hagamos (PDF, HTML,etc), las conclusiones que saquemos, etc. Dependiendo de la complejidad del proyecto puede separar en carpetas al interior si hay outputs muy variados.

SRC

Carpeta para todo nuestro código. Algunos eligen no usarla y dejar los códigos en la ruta del proyecto pero me parece un poco desprolijo. Recomendamos tener muchos scripts con títulos claros y segmentados por lo que hacen. Es decir, uno para levantar la data, otro para análisis exploratorio, otro para feature engineering y así. A su vez, recomendamos tener un script propio para la funciones que definan ustedes y si les resulta cómodo otro para las librerías, de manera tal de tener todo claro, separado y no tener que andar buscando dentro de un grán código lo que necesitan. Además es más sencillo para modificar y arreglar bugs. Recuerden que para invocar código de otro script simplemente lo corren usando source("Script.R"). Dejamos un ejemplo ilustrativo.

GITHUB

En algún otro post lo dijimos pero recomendamos altamente UTILIZAR GITHUB para manejar sus proyectos, tener backups, compartirlos y actualizarlos desde cualquier computadora! Y hacer blogs como este siguiendo este POST

Lo que también sugerimos es no subir la carpeta data a github por dos motivos. Primero por una cuestión de espacio, si tienen data muy pesada Github no les va a permitir incluirla en el repositorio. Por otra parte si la data es confidencial o tiene datos privados mejor que no esté a disposición de cualquiera si tienen cuenta pública. Obviamente queda a criterio y comodidad de cada uno si corresponde subir la data o no. Para evitar que una carpeta sea subida a github solo deben incluirla en el archivo .gitignore de su repositorio.

2018-10-28
in estadistica, R
3 min read

Curvas ROC

La curva ROC y AUC (area bajo la curva) permiten evaluar la eficacia de un modelo clasificador y elegir el mejor umbral de corte donde determinar qué observación es predicha positiva y cual negativa.

Vamos a generar rapidamente un clasificador con regresión logísitca utilizando el dataset mtcars ya provisto por R. Solo a modo ilustrativo utilizaremos AM (caja manual o automática) como la variable a predecir y mpg y drat como independientes. No separamos en train y test dadas las pocas observaciones.

library(tidyverse)
library(modelr)
library(pROC)
df <- mtcars %>% select(am, mpg, drat) %>% mutate(am = as.factor(am))
summary(df)

##  am          mpg             drat      
##  0:19   Min.   :10.40   Min.   :2.760  
##  1:13   1st Qu.:15.43   1st Qu.:3.080  
##         Median :19.20   Median :3.695  
##         Mean   :20.09   Mean   :3.597  
##         3rd Qu.:22.80   3rd Qu.:3.920  
##         Max.   :33.90   Max.   :4.930

# Clase dentro de todo balanceada

mdl.log <- glm(formula = am ~., data = df, family = binomial(link="logit"))
fit <- predict(mdl.log, newdata = df, type = "response")


roc(df[,1],  fit , percent=F,   boot.n=1000, ci.alpha=0.9, stratified=FALSE, plot=TRUE, grid=TRUE, show.thres=TRUE, legacy.axes = TRUE, reuse.auc = TRUE,
    # print.thres = c(0.30,0.35, 0.40, 0.45,0.48, 0.50,0.55, 0.60),#
    print.auc = TRUE, print.thres.col = "blue", ci=TRUE, ci.type="bars", print.thres.cex = 0.7, main = paste("ROC curve using","(N = ",nrow(df),")") )

## 
## Call:
## roc.default(response = df[, 1], predictor = fit, percent = F,     ci = TRUE, plot = TRUE, boot.n = 1000, ci.alpha = 0.9, stratified = FALSE,     grid = TRUE, show.thres = TRUE, legacy.axes = TRUE, reuse.auc = TRUE,     print.auc = TRUE, print.thres.col = "blue", ci.type = "bars",     print.thres.cex = 0.7, main = paste("ROC curve using", "(N = ",         nrow(df), ")"))
## 
## Data: fit in 19 controls (df[, 1] 0) < 13 cases (df[, 1] 1).
## Area under the curve: 0.9433
## 95% CI: 0.8695-1 (DeLong)

Básicamente entrenamos un modelo logístico y graficamos la curva ROC prediciendo sobre el mismo dataset con el que fue entrenado. No es lo adecuado pero dadas las pocas observaciones y el propósito explicativo no lo tomamos como un problema. La curva ROC es la más oscura y como vemos empieza en (0,0) y termina en el (1,1). El eje X es 1 - Especificidad (Falsos Negativos) y el eje Y es Sensitividad (Verdaderos Positivos) por lo tanto lo deseable es estar lo más arriba a la izquierda posible. El punto (0,1) sería óptimo ya que habría 0 falsos negativos y 100% de verdaderos positivos.

Lo que representa la curva es la combinación de Sensitividad y (1 - especificidad) para varios puntos de corte. Recordemos que la regresión logística devuelve un valor entre 0 y 1 por lo tanto hay que determinar en qué valor empezamos a considerar una predicción como positiva o negativa. En este caso positivo sería tener un valor de 1 en am, por lo tanto tener caja automática. Cada punto de la curva corresponde a algún punto de corte. Como decíamos antes, el mejor debería ser el más "arriba a la izquierda" aunque depende el problema eso puede cambiar, dependiendo del costo de equivocarse en uno u otro sentido.

El peor escenario es que la curva siga a la diagonal, lo que equivaldría a ser iguales a un modelo eligiendo siempre la clase mayoritaria, totalmente inútil. Si estuviera por debajo de la diagonal, sería peor aún, pero bastaría con invertir las predicciones para pasar a estar por encima. Un viejo truco no muy científico.

El área bajo la curva (AUC) es una medida resumen de la curva ROC ya que justamente describe el área entre la curva ROC y la diagonal. Valores mayores se corresponden con curvas ROC más alejadas de la diagonal y por lo tanto que separan mejor a la clase dependiente. Es útil para comparar modelos.

2018-10-07
in R, blog
3 min read

Como crear un blog con Blogdown y Netlify

En este post vamos a ver el proceso reusmido para crear un blog donde podemos generar contenido directamente desde RStudio, utilizando el paquete blogdown, Github y Netlify. La idea es que al finalizar la configuración, simplemente creemos el post en Rstudio (un markdown) usando blogdown y que al subirlo a github automaticamente se actualice el blog y se vea reflejado en nuestra página. Es lo que estoy haciendo en este momento.

Lo primero y más importante más allá de este tema en particular es tener una cuenta en github. Si no la tienen se los recomiendo ampliamente para hacer version control - actualizar código de forma segura y con backups constantes en un servidor + compartir proyectos. Es gratis, al menos la versión básica que alcanza y sobra para el uso cotidiano. En nuestro repositorio crearemos un proyecto para el blog y ahi se subirán nuestros posts en formato html. A su vez, estará el theme y otras configuraciones básicas del blog.

Lo segundo es instalar el paquete blogdown en nuestro R. Es una obra maestra de Yihui Xie, ingeniero de RStudio y creador de varios paquetes. En un nuevo proyecto de R, ejecutan el siguiente comando.

blogdown::new_site()

Con eso ya tienen generado la estructura básica de lo que será su blog. Se crearán carpetas y archivos en la ruta del proyecto con contenido de prueba para tener algo funcional. Blogdown es bastante complejo y hay un millón de configuraciones y detalles que uno puede personalizar. No entraremos en eso acá porque se haría super extenso. Para eso está la documentación oficial en ingles.

Eventualmente van a tener que cambiar el archivo config.toml con ciertos pasos de la guía y luego podrán explorar todas las posibilidades que presenta. Entre ellas pueden (y recomiendo) descargar otro theme para cambiar el formato. El que uso actualmente es tranquilpeak.

Para crear un post nuevo simplemente escriben el comando.

blogdown::new_post()

Lo cual genera un script con un YAML, que es la configuración con el título, tags y otros metadatos del post. Simplemente escriben como cualquier markdown debajo. Cuando terminan el post (o mientras para ir visualizando como queda) corren

blogdown::serve_site()

Lo cual generará el archivo html correspondiente que luego será usado en su web y les permite ver el resultado temporal de su post.

Luego lo que deben hacer es pushear la carpeta que se les generó del blog a su repositorio en github. O al menos las carpetas y archivos que se ven en la siguiente imagen. Public no la pusheen.

gitsnap

Llegado a este punto tenemos el contenido inicial del blog, pero no hay sitio web. Ahí es donde entra en juego Netlify. No vamos a entrar en el paso a paso minucioso pero básicamente deben crearse una cuenta, generar una nueva web y linkearla a su repositorio Github. Es bastante lineal. Luego configuran el nombre de la web y otros detalles y en cuestión de minutos ya están publicados! (Y Gratis.) Para ver bien esta etapa les recomiendo la explicación del link de blogdown.

Una vez puesto en marcha solo es cuestión de abrir su proyecto en R (en su pc), crear un nuevo post con blogdown::new_post() y pushear a github! Prueben chusmeando todas las configuraciones para cambiar la estética del blog!

2018-09-01
in R
3 min read

Hola MAP. Chau Apply

Introducción

La idea de este post es introducirlos a la familia de funciones MAP, propias de tidyverse. A grandes rasgos son un remplazo MUY útil a la familia de funciones APPLY, propias de R base. Estas últimas se suelen enseñar en todos los cursos introductorios de R, como la manera correcta de aplicar funciones a listas o columnas de dataframes. No es que no sirvan, pero dado el surgimiento de tantas librerías que facilitan el manejo de la data, no tiene sentido seguir insistiendo con ellas dado que hay nuevas con mayor flexibilidad, muy sencillas de utilizar y mucho más amenas.

Lo mejor que tienen las funciones MAP es:
- Consistencia en los inputs.
- Flexibilidad del output.
- Integración con todo el universo tidyverse y prolijidad.

Empecemos.

library(purrr) # MAP está contenida acá
library(dplyr)

## Warning: package 'dplyr' was built under R version 4.0.5

Como regla general, MAP aplica funciones a elementos de una lista o de un vector. Su output es otra lista. Muy similar a lapply().

l1 <- list( a= c(100,200), b = c(8,10))
map(l1, max)

## $a
## [1] 200
## 
## $b
## [1] 10

A cada lista le calcula el máximo y devuelve una lista con cada elemento siendo el resultado de la función.

Tenemos la flexibilidad para pasarle funciones anónimas..

map(l1, function(x) max(x))

## $a
## [1] 200
## 
## $b
## [1] 10

Aplicando funciones a elementos de un vector. Cada numero de 1 a 5 es usado como primer input de la funcion rnorm, sd y n son otros parámetros de rnorm. El resultado de nuevo es una lista.

set.seed(1)
1:5 %>% map(., rnorm,sd =2, n=5)

## [[1]]
## [1] -0.2529076  1.3672866 -0.6712572  4.1905616  1.6590155
## 
## [[2]]
## [1] 0.3590632 2.9748581 3.4766494 3.1515627 1.3892232
## 
## [[3]]
## [1]  6.023562  3.779686  1.757519 -1.429400  5.249862
## 
## [[4]]
## [1] 3.910133 3.967619 5.887672 5.642442 5.187803
## 
## [[5]]
## [1] 6.837955 6.564273 5.149130 1.021297 6.239651

Consistencia entre variantes

Por ahora solo vimos la versión de lapply en MAP, pero esta familia tiene varios integrantes.

map_if

Ejecuta la función solo si el elemento cumple determinada condición. Devuelve una lista.

l2 <- list(a = 213, b = "string", c = c(1,2))
map_if(l2, is.numeric, function(x) x*2)

## $a
## [1] 426
## 
## $b
## [1] "string"
## 
## $c
## [1] 2 4

El output es la lista original con los elementos correspondientes transformados. Vemos que no hubo ningún problema con "string" ya que fue omitido.

map_at

Ejecuta la función solo en los elementos que seleccionemos. No hace falta que cumplan alguna condición. Misma función de antes pero solo aplicada al tercer elemento. Devuelve una lista.

map_at(l2, c(3), function(x) x*2)

## $a
## [1] 213
## 
## $b
## [1] "string"
## 
## $c
## [1] 2 4

Variantes súper útiles que permiten no utilizar loops y que dan mucho control de manera sencilla sobre las funciones a ejecutar. Por otra parte, en términos de consistencia, la estructura es siempre la misma. El primer argumento es x= y luego viene la función a aplicar. En el caso de map_if y map_at entre medio surge el condicionante. Si recuerdan, la familia apply cambia el orden de los inputs según si es apply, lapply, mapply, sapply...

Flexibilidad del output

Por el momento vimos que todos los outputs eran listas. Lo interesante es que podemos controlar eso y cambiar el formato del resultado, ahorrándonos conversiones molestas con unlist y etc.

l3 <- list(c(1,2,4), c(100,200), c(5000,6000))
map_dbl(l3, max)

## [1]    4  200 6000

Nos devuelve un vector con los resultados de aplicar la función max a cada elemento!

De este mismo tipo esta.

map_chr # vector caracter
map_int # vector de integers
map_lgl # vector de booleanos