Informe.utf8.md

Universidad CAECE

Taller de Lenguajes I

Profesora

Silvia A. Cobialca

Informe por

Andrés Garcia Alves (leg. 1033744)

El objetivo de la materia es el aprendizaje por parte del alumno del Lenguaje R, y su utilización en el análisis estadístido de datos.

Indice

Capítulo 1: Investigación

En la primera parte se investigan conceptos tanto de programación en R, como de Estadística.
Estos se aplicarán en la 2da parte del informe a un proyecto concreto.

Capítulo 2: Análisis de los Datos

En la segunda parte se realiza un análisis estadístico sobre muestras de nivel de Ph, unidas estas a datos climáticos.
Se aplica regresión lineal para intentar determinar posibles asociaciones entre estos.

Anexo 1: Script Datos del Clima (SMN)

Se presenta aquí un script para parsear datos en crudo del clima, disponibles en formato TXT desde la web del Servicio Meteorológico Nacional.

Anexo 2: Script Datos de pH

Se presenta aquí un script para parsear datos en crudo de las mediciones de pH, disponibles en formato CSV.

Anexo 3: Proceso de Medición de pH

Se presenta aquí una pequeña investigación sobre el proceso de medición del pH y cómo funciona un pH-metro.

Historial de Versiones

Detalle de avances del presente informe.

Capítulo 1: Investigación

Markdown

Markdown es un lenguaje de marcado ligero creado por John Gruber en el año 2004 que trata de conseguir la máxima legibilidad y facilidad de publicación, tanto en su forma de entrada como de salida, inspirándose en muchas convenciones existentes para marcar mensajes de correo electrónico usando texto plano. [Wikipedia: Markdown]

Saltos de línea: Los saltos de línea se generan cuando se encuentran dos espacios juntos.

"Este es un texto,  
separado en dos líneas"

Encabezados: Los encabezados se generan cuando se encuentra una almohadilla antes de texto.

# Encabezado h1 
## Encabezado h2
### Encabezado h3
#### Encabezado h4
##### Encabezado h5
###### Encabezado h6

Citas: Para citar solo es necesario escribir una cuña antes del texto.

> "Hay una fuerza motriz más poderosa que el vapor, la electricidad y la energía atómica: la voluntad". -Albert Einstein

Texto con énfasis: envolver entre un asterisco para cursiva y dos para negrita.

*énfasis* (cursiva)
**énfasis fuerte** (negrita)

Subíndices: envolver entre ‘~’.
H~2~O >> H₂O

Superíndices: envolver entre ‘^’ para superíndice.
2^3^=8 >> 2³=8

Código: Se utiliza el acento grave para identificar código, y corchetes para identificar el lenguaje de programación.

ˋˋˋ
Código en 
varias líneas
ˋˋˋ

Listas:

* Un elemento en una lista no ordenada
* Otro elemento en una lista

1. Elemento en una lista enumerada u ordenada.
2. Otro elemento

Enlaces:

[Google](http://www.google.com.ar "Título del enlace")

Imágenes:

![Texto Alt.](path/imagen.png "Titulo")

Letras Griegas:
Rodear su nombre de la siguiente forma $\mu$ –> $\mu$, $\pi$ –> $\pi$.

Tipos de Objetos

La siguiente tabla resume los tipos de objetos y los datos que representan:

Objeto	Tipos	Puede Mezclar Tipos?
vector	numérico, caracter, complejo o lógico	No
arreglo	numérico, caracter, complejo o lógico	No
matriz	numérico, caracter, complejo o lógico	No
factor	numérico o caracter	No
data.frame	numérico, caracter, complejo o lógico	Si
ts	numérico, caracter, complejo o lógico	Si
lista	numérico, caracter, complejo o lógico	Si

Un ‘vector’ es una variable en el significado comunmente asumido. Un ‘factor’ es una variable categórica. Un arreglo es una tabla de dimensión k, y una ‘matriz’ es un caso particular de un arreglo donde k=2. Notar que los elementos en un arreglo o una matriz son del mismo tipo. Un ‘data.frame’ es una tabla compuesta de uno o más vectores y/o factores de la misma longitud pero que pueden ser de diferentes tipos. Un ’ts’ es una serie temporal y como talcontiene atributos adicionales tales como frecuencia y fechas. Finalmente, una ‘lista’ puede contenercualquier tipo de objeto incluyendo otras listas. [CRAN, 'R para Principiantes', E. Paradis]

Vectores

Se le denomina vector, formación o arreglo (en inglés array) a una zona de almacenamiento contiguo que contiene una serie de elementos del mismo tipo. [Wikipedia: Vector]

En específico a R, un vector es una esctructura que representa una secuencia de datos del mismo tipo, identificados por un subíndice con base en uno.

# creación de vectores con rangos numéricos
vector1 <- c(1:3)      # 1:3 es un rango de núm. enteros del 1 al 3
vector2 <- c(4:6)      # 4:6 es un rango de núm. enteros del 4 al 6
vector3 <- c(7:9)      # 1:3 es un rango de núm. enteros del 7 al 9

vector1               # contenido del vector1

## [1] 1 2 3

vector1[2]            # 2do elemento del vector1

## [1] 2

Data Frame

Un data frames es un estructura de datos de dos dimensiones (tablas).

Estos están compuestos por observaciones (representados en las filas), y variables que son medidas sobre estos (representadas en las columnas).

Para más información, véase [R Tutorial], o [R para principiantes].

R provee un amplio conjunto de dataframes por default a fines didácticos. Los mismos se pueden acceder en datasets::nombre-del-dataframe, o en forma abreviada simplemente con nombre-del-dataframe.

# construcción de un data.frame con los vectores de la sección anterior
data.frame(vector1, vector2, vector3)

##   vector1 vector2 vector3
## 1       1       4       7
## 2       2       5       8
## 3       3       6       9

# visualizar el dataframe datasets::mtcars
mtcars

##                      mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive      21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
## Valiant             18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
## Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
## Merc 240D           24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
## Merc 230            22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
## Merc 280            19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
## Merc 280C           17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
## Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
## Merc 450SL          17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
## Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
## Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
## Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
## Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
## Fiat 128            32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
## Honda Civic         30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
## Toyota Corolla      33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
## Toyota Corona       21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
## Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
## AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
## Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
## Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
## Fiat X1-9           27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
## Porsche 914-2       26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
## Lotus Europa        30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
## Ford Pantera L      15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
## Ferrari Dino        19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
## Maserati Bora       15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8
## Volvo 142E          21.4   4 121.0 109 4.11 2.780 18.60  1  1    4    2

Para referenciar a una variable dada del dataframe se utiliza el caracater ‘$’, quedando entonces de la siguiente forma: nombre-del-dataframe$nombre-de-la-variable

mtcars$mpg

##  [1] 21.0 21.0 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 17.8 16.4 17.3 15.2
## [15] 10.4 10.4 14.7 32.4 30.4 33.9 21.5 15.5 15.2 13.3 19.2 27.3 26.0 30.4
## [29] 15.8 19.7 15.0 21.4

Promedio

También conocido como media aritmética, es el valor característico de una serie de datos cuantitativos. [Wikipedia: Media Aritmética]
Se representa con la letra griega Mu ‘$\mu$’.

Para calcular el promedio se utiliza mean(x) sobre una secuencia numérica.

mean(mtcars$mpg)      # promedio

## [1] 20.09062

Desviación Estándar

Es la medida de la variación o dispersión de un conjunto de datos numérico. [Wikipedia: Desviación Típica]
Se representa con la letra griega Sigma ‘$\sigma$’.

Para calcular la desviación estándar se utiliza sd(x) sobre una secuencia numérica.

sd(mtcars$mpg)        # desviacion estandar

## [1] 6.026948

Histograma

Es la representación gráfica de una variable en forma de barras. [Wikipedia: Histograma]

Para visualizar un histograma se utilizará hist(x, ...)
Se puede elegir la cantidad de separaciones deseadas (en cuartiles, deciles, percentiles, etc) a través del 2do parámetro.

hist(mtcars$carb, 10, col=col_caece)  # histograma en deciles

De la gráfica se pueden realizar algunas observaciones:

La frecuencia máxima de aparición es 10 veces, observada en los valores 2 y 4.
Casi todas las mediciones se situan en el rango de valores 1 al 4.

Diagrama de Dispersión

Un diagrama de Dispersiónes un tipo de diagrama matemático, que utiliza las coordenadas cartesianas para mostrar los valores de dos variables para un conjunto de datos. [Wikipedia: Diagrama de Dispersión]

Se generar un diagrama de dispersión se utilza pairs(x).

pairs(mtcars)         # diagrama de dispersión

Estructura

Se puede visualizar la estructrura de cualquier objeto mediante la función str(x).

str(mtcars)           # e(STR)uctura del dataframe mtcars

## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

En la primera línea indica que ese objeto es del tipo ‘data.frame’, con 32 observaciones, con 11 variables en cada una.
En las filas sucesivas muestra los primeros 10 valores y su tipo de dato para cada una de las variables del dataframe.

Sumario

Se puede visualizar un resumen de cualquier objeto mediante la funcion summary(x)

summary(mtcars$carb)  # sumario de mtcars$carb

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   2.000   2.000   2.812   4.000   8.000

Donde para este dataframe$variable arroja que:

El valor mínimo es 1.
El valor máximo es 8.
La mediana es 2.
El promedio es 2,812.
El Cuartil #1 termina en 2.
El Cuartil #3 termina en 4.

Primeros & Ultimos

Se pueden visualizar tanto las n primeras, como n últimas observaciones de un dataframe mediantes el uso de las funciones head(x) y tail(x), respectivamente.
Por default se incluirán hasta 6 elementos, pero se puede modificar este comportamiento mediante el 2do parámetro.

head(mtcars)          # primeras 6 observaciones

##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

tail(mtcars, 10)      # últimas 10 observaciones

##                   mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## AMC Javelin      15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
## Camaro Z28       13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
## Pontiac Firebird 19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
## Fiat X1-9        27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
## Porsche 914-2    26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
## Lotus Europa     30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
## Ford Pantera L   15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
## Ferrari Dino     19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
## Maserati Bora    15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8
## Volvo 142E       21.4   4 121.0 109 4.11 2.780 18.60  1  1    4    2

Valores Atípicos (Outliers)

Es una observación que es numéricamente distante del resto de los datos. Las estadísticas derivadas de los conjuntos de datos que incluyen valores atípicos serán frecuentemente engañosas. [Wikipedia: Valor Atípico]

En otras palabras, los valores atípicos en los datos pueden distorsionar las predicciones y afectar la precisión. Estos deberían ser descartados previo a relizar cualquier análisis.

Diagrama de Caja (Box-Plot)

Es un método estandarizado para representar gráficamente una serie de datos numéricos a través de sus cuartiles. [Wikipedia: Diagrama de Caja]

Se utiliza para la detección de Valores Atípicos, presentados en este gráfico como pequeños círculos en los extremos.

Cuartil #1: desde la línea superior del gráfico hasta el inicio de la caja.
Cuartil #2: desde el inicio de la caja hasta la línea del medio de la caja (la mediana).
Cuartil #3: desde la línea del medio de la caja hasta el fin de la caja.
Cuartil #4: desde el fin de la caja hasta la línea inferior del gráfico.

par(mar=c(2,2,2,2))
boxplot(mtcars$hp, col=col_caece)  # distribucuión en cuartiles

En el margen superior de la gráfica se puede visualizar como el dataframe mtcars, sobre su variable $hp, incluye una observación que presenta un valor atípico.

par(mar=c(2,2,2,2))
boxplot(mtcars, col=col_caece)

De la gráfica anterior se observa que en general no se pueden comparar de forma directa las variables de un dataframe debido a las diferencias de escala. Previo a su comparación se requiere normalizar los datos (ver siguiente sección).

Normalización de Datos

La normalización de índices es el proceso de ajustar los valores medidos en diferentes escalas respecto a una escala común. [Wikipedia: Normalización]

Si bien existen varios procesos, la normalización estándar viene dada por la fórmula: Normalización Estandar

Es decir, cada observación ‘X’ menos la media ‘$\mu$’, dividido sobre la desviación estandar ‘$\sigma$’.

Combinando las funciones vistas hasta el momento, en unas pocas líneas de código se pueden normalizar algunas variables:

# en este ejemplo se normalizarán las variables 'mtcars$hp' y 'mtcars$wt'

mu = mean(mtcars$hp)                    # promedio de mtcars$hp
sigma = sd(mtcars$hp)                   # desviación estandar de mtcars$hp
normal_hp = (mtcars$hp - mu) / sigma    # normalizacion de mtcars$hp

mu = mean(mtcars$wt)                    # promedio de mtcars$wt
sigma = sd(mtcars$wt)                   # desviación estandar de mtcars$wt
normal_wt = (mtcars$wt - mu) / sigma    # normalizacion de mtcars$wt

# construcción de un nuevo dataframe con ambas variables ya normalizadas
datos <- data.frame(normal_hp, normal_wt)
datos                                   # visualizar el dataframe

##      normal_hp    normal_wt
## 1  -0.53509284 -0.610399567
## 2  -0.53509284 -0.349785269
## 3  -0.78304046 -0.917004624
## 4  -0.53509284 -0.002299538
## 5   0.41294217  0.227654255
## 6  -0.60801861  0.248094592
## 7   1.43390296  0.360516446
## 8  -1.23518023 -0.027849959
## 9  -0.75387015 -0.068730634
## 10 -0.34548584  0.227654255
## 11 -0.34548584  0.227654255
## 12  0.48586794  0.871524874
## 13  0.48586794  0.524039143
## 14  0.48586794  0.575139986
## 15  0.85049680  2.077504765
## 16  0.99634834  2.255335698
## 17  1.21512565  2.174596366
## 18 -1.17683962 -1.039646647
## 19 -1.38103178 -1.637526508
## 20 -1.19142477 -1.412682800
## 21 -0.72469984 -0.768812180
## 22  0.04831332  0.309415603
## 23  0.04831332  0.222544170
## 24  1.43390296  0.636460997
## 25  0.41294217  0.641571082
## 26 -1.17683962 -1.310481114
## 27 -0.81221077 -1.100967659
## 28 -0.49133738 -1.741772228
## 29  1.71102089 -0.048290296
## 30  0.41294217 -0.457097039
## 31  2.74656682  0.360516446
## 32 -0.54967799 -0.446876870

Y ahora ya se puede proceder a compararlas:

par(mar=c(2,2,2,2))
boxplot(datos, col=col_caece)

Normalización de Datos, Parte II

Sabiendo ya el proceso de normalización “a mano”, es importante conocer que se le puede pedir a R que lo realice por nosotros, ya sea tanto sobre una variable, como sobre un dataframe completo.

Para realizar la normalización se debe utilizar la función scale(x).

par(mar=c(2,2,2,2))
normal_hp <- scale(mtcars$hp)
normal_wt <- scale(mtcars$wt)
datos <- data.frame(normal_hp, normal_wt)
boxplot(datos, col=col_caece)

Ahora normalizando el dataframe completo:

par(mar=c(2,2,2,2))
datos <- scale(mtcars)
boxplot(datos, col=col_caece)

Importación de Datos

Una forma sencilla de importar datos en R es a través de la clase read.<formato>(path).

read.csv(file="recursos/test.csv", header=TRUE)

##   nombre edad altura
## 1  maria   35   1.65
## 2   josé   32   1.71
## 3  laura   28   1.60
## 4 carlos   25   1.77

Desnormalización de Datos

La desnormalización de datos no es más que aplicar la formula de la normalización a la inversa.

Es decir, cada observación normalizada ‘X_n’ multiplicado por la desviación estandar ‘$\sigma$’, más la media ‘$\mu$’.

mu = mean(mtcars$hp)                    # promedio de mtcars$hp
sigma = sd(mtcars$hp)                   # desviación estandar de mtcars$hp

normal_hp = (mtcars$hp - mu) / sigma    # normalización
de_normal_hp = (normal_hp * sigma) + mu # desnormalización

mtcars$hp                               # los valores originales

##  [1] 110 110  93 110 175 105 245  62  95 123 123 180 180 180 205 215 230
## [18]  66  52  65  97 150 150 245 175  66  91 113 264 175 335 109

normal_hp                               # los valores nomalizados

##  [1] -0.53509284 -0.53509284 -0.78304046 -0.53509284  0.41294217
##  [6] -0.60801861  1.43390296 -1.23518023 -0.75387015 -0.34548584
## [11] -0.34548584  0.48586794  0.48586794  0.48586794  0.85049680
## [16]  0.99634834  1.21512565 -1.17683962 -1.38103178 -1.19142477
## [21] -0.72469984  0.04831332  0.04831332  1.43390296  0.41294217
## [26] -1.17683962 -0.81221077 -0.49133738  1.71102089  0.41294217
## [31]  2.74656682 -0.54967799

de_normal_hp                            # los valores desnormalizados

##  [1] 110 110  93 110 175 105 245  62  95 123 123 180 180 180 205 215 230
## [18]  66  52  65  97 150 150 245 175  66  91 113 264 175 335 109

Correlación Lineal, Parte I: Teoría

Para estudiar la relación lineal existente entre dos variables continuas es necesario disponer de parámetros que permitan cuantificar dicha relación. Uno de estos parámetros es la covarianza, que indica el grado de variación conjunta de dos variables aleatorias.
[R Pubs: Correlación lineal y Regresión lineal simple]

Definición:
Covarianza Muestral

siendo $\bar{x}$ e y $\bar{y}$ la media de cada variable y x_i e y_i el valor de las variables para la observación i.

La covarianza depende de las escalas en que se miden las variables estudiadas, por lo tanto, no es comparable entre distintos pares de variables. Para poder hacer comparaciones se estandariza la covarianza, generando lo que se conoce como coeficientes de correlación. Existen diferentes tipos, de entre los que destacan el coeficiente de Pearson, Rho de Spearman y Tau de Kendall.

Todos ellos varían entre +1 y -1. Siendo +1 una correlación positiva perfecta y -1 una correlación negativa perfecta.

Se emplean como medida de fuerza de asociación (tamaño del efecto):

0.0: asociación nula.
0.1: asociación pequeña.
0.3: asociación mediana.
0.5: asociación moderada.
0.7: asociación alta.
0.9: asociación muy alta.

Además del valor obtenido para el coeficiente de correlación, es necesario calcular su significancia. Solo si el p-value es significativo se puede aceptar que existe correlación, y esta será de la magnitud que indique el coeficiente. Por muy cercano que sea el valor del coeficiente de correlación a +1 o -1, si no es significativo, se ha de interpretar que la correlación de ambas variables es 0, ya que el valor observado puede deberse a simple aleatoriedad.

El test paramétrico de significancia estadística empleado para el coeficiente de correlación es el t-test.

Significancia

El nivel de significación es un concepto asociado a la verificación de una hipótesis. En pocas palabras, se define como la probabilidad de tomar la decisión de rechazar la hipótesis nula cuando ésta es verdadera (decisión conocida como error de tipo I).

La decisión se toma a menudo utilizando el valor-p: si este es inferior al nivel de significación, entonces la hipótesis nula es rechazada. Cuanto menor sea el valor p, más significativo será el resultado. Valores comunes de significación: 0.05, 0.01 y 0.001.

[Wikipedia: Significación estadística] [Wikipedia: Valor P]

Correlación Lineal, Parte II: Práctica

Para calcular en R la Correlación Lineal entre dos variables se emplea el comando cor.test(). En este, además de calcular el valor de la correlación, también indicará su Significancia.

En la siguiente prueba intentaré estudiar el dataframe mtcars y verificar si existe una correlación entre mtcars$wt (peso) y mtcars$mpg (millas por galón) en las observaciones incluidas.

cor.test(x = mtcars$wt, y = mtcars$mpg, alternative = "two.sided", conf.level = 0.95, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  mtcars$wt and mtcars$mpg
## t = -9.559, df = 30, p-value = 1.294e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9338264 -0.7440872
## sample estimates:
##        cor 
## -0.8676594

Dado el cor.test() anterior puedo concluir que gracias al bajo p-value (valor 0.0000000001294), se corrobora la correlación alta (coeficiente de 0.86) entre el peso y las millas por galón, para la muestra dada.

Adicionalmente se podría pedir a R que calcule las correlaciones de todo un dataframe simplemente con el comando cor().
Comando complementario al uso de `pairs(), pero para llevar a cabo un análisis analítico en lugar de visual.

cor(mtcars)

##             mpg        cyl       disp         hp        drat         wt
## mpg   1.0000000 -0.8521620 -0.8475514 -0.7761684  0.68117191 -0.8676594
## cyl  -0.8521620  1.0000000  0.9020329  0.8324475 -0.69993811  0.7824958
## disp -0.8475514  0.9020329  1.0000000  0.7909486 -0.71021393  0.8879799
## hp   -0.7761684  0.8324475  0.7909486  1.0000000 -0.44875912  0.6587479
## drat  0.6811719 -0.6999381 -0.7102139 -0.4487591  1.00000000 -0.7124406
## wt   -0.8676594  0.7824958  0.8879799  0.6587479 -0.71244065  1.0000000
## qsec  0.4186840 -0.5912421 -0.4336979 -0.7082234  0.09120476 -0.1747159
## vs    0.6640389 -0.8108118 -0.7104159 -0.7230967  0.44027846 -0.5549157
## am    0.5998324 -0.5226070 -0.5912270 -0.2432043  0.71271113 -0.6924953
## gear  0.4802848 -0.4926866 -0.5555692 -0.1257043  0.69961013 -0.5832870
## carb -0.5509251  0.5269883  0.3949769  0.7498125 -0.09078980  0.4276059
##             qsec         vs          am       gear        carb
## mpg   0.41868403  0.6640389  0.59983243  0.4802848 -0.55092507
## cyl  -0.59124207 -0.8108118 -0.52260705 -0.4926866  0.52698829
## disp -0.43369788 -0.7104159 -0.59122704 -0.5555692  0.39497686
## hp   -0.70822339 -0.7230967 -0.24320426 -0.1257043  0.74981247
## drat  0.09120476  0.4402785  0.71271113  0.6996101 -0.09078980
## wt   -0.17471588 -0.5549157 -0.69249526 -0.5832870  0.42760594
## qsec  1.00000000  0.7445354 -0.22986086 -0.2126822 -0.65624923
## vs    0.74453544  1.0000000  0.16834512  0.2060233 -0.56960714
## am   -0.22986086  0.1683451  1.00000000  0.7940588  0.05753435
## gear -0.21268223  0.2060233  0.79405876  1.0000000  0.27407284
## carb -0.65624923 -0.5696071  0.05753435  0.2740728  1.00000000

Regresión Lineal, Parte I: Teoría

El objetivo de un modelo de regresión es tratar de explicar la relación que pueda existir entre una variable dependiente Y (variable respuesta) un conjunto de variables independientes X1, X2, … Xn (variables explicativas).
En un modelo de regresión lineal simple se trata de explicar la relación que existe entre la variable respuesta Y y una única variable explicativa X. [Univ. de Santiago de Compostela]

Requisitos:
Para poder crear un modelo de regresión lineal es necesario que se cumpla con los siguientes supuestos:

Linealidad: Que la relación entre las variables sea lineal.
Que los errores en la medición de las variables explicativas sean independientes entre sí.
Homocedasticidad: Que los errores tengan varianza constante.
Distribución normal de los residuos: Que los errores tengan una esperanza matemática igual a cero.
Que el error total sea la suma de todos los errores.
Se debe estudiar con detenimiento los valores atípicos ya que pueden generar una falsa correlación que realmente no existe, u ocultar una existente.

[Wikipedia: Regresión Lineal, sección “Supuestos del modelo de regresión lineal”]

Un modelo de regresión lineal se describe con una ecuación muy similar a una recta:

Siendo β₀ la ordenada en el origen, β₁ la pendiente y ϵ el error aleatorio (residuo). Este último representa la diferencia entre el valor ajustado por la recta y el valor real.

Por norma general, los estudios de correlación lineal preceden a la generación de modelos de regresión lineal. Primero se analiza si ambas variables están correlacionadas y, en caso de estarlo, se procede a generar el modelo de regresión.
[R Pubs: Correlación lineal y Regresión lineal simple]

Regresión Lineal, Parte II: Práctica

Para efectuar en R un modelo de Regresión Lineal se emplea el comando modelo <- lm(var1~var2+var3+varN, dataframe).
Donde var1 será la variable dependiente y desde var2 en adelante las variables independientes.

modelo <- lm(mpg~wt, mtcars)
summary(modelo)

## 
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5432 -2.3647 -0.1252  1.4096  6.8727 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
## wt           -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7446 
## F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

La columna Estimate devuelve el valor estimado para los dos parámetros de la ecuación del modelo lineal (β₀ y β₁) que equivalen a la ordenada en el origen y la pendiente respectivamente. La última columna Pr(>|t|) muestra el p-value (significancia) en ambas variables.

Para el modelo generado, tanto la ordenada en el origen como la pendiente son significativas (p-values < 0.05).

plot(mtcars$wt, mtcars$mpg, pch = 16, col = col_caece, main = "Peso vs. Rendimiento", xlab = "Peso (miles de libras)", ylab = "Millas Por Galón (US)")
abline(modelo)

Se visualiza de la gráfica que a mayor peso (wt), menor millas por galón (mpg).
Referencias: 1 libra = 0.453592 kgs, 1 milla = 1.609344 kms, 1 galón (US) = 3.785412 lts.

Validación del modelo (Residuos & QQ-Normal)

Los residuos (o errores) son la diferencia entre los valores observados y los valores que predice el modelo.
[Maxima Formacion: Como validar tu model de Regresión] [RPus: Regresion Simple]

Residuos = Valores observados - Valores que predice el modelo:
e = y – ŷ

Para validar la homocedasticidad y la linealidad del modelo se suele utilizar un gráfico de residuos.
En este, mediante una inspección visual se debe verificar la aleatoreidad de los mismos.

residuos <- residuals(modelo)

par(mar=c(2,2,2,2))
plot(residuos, pch = 16, col = col_caece)

Luego, para validar la hipótesis de normalidad se suele realizar mediante un gráfico QQ-Normal, también sobre los residuos. Este gráfico mostrará una comparación cuartil a cuartil de los datos recibidos contra una distribución normal teórica.
De verificarse, se debe observar que los residuos siguen aproximadamente la línea recta diagonal.

par(mar=c(2,2,2,2))
qqnorm(residuos, pch = 16, col = col_caece)
qqline(residuos)

Regresión Lineal, Parte III: Práctica

En el apartado Regresión Lineal Pt. II pretendí, como primera aproximación, contrastar una posible relación lineal del rendimiento en función del peso. En aquél caso, la hipótesis nula (H₀) que indicaba la no relación resultó abrumadoramente rechazada, con un p-value de 1.294*10^-10, y aportando de esta manera, evidencia hacia la veracidad de la existencia de la relación (H₁).

Como mejora, en este apartado me propongo validar que la variable wt (peso) NO sea colineal con el resto de las variables del dataframe mtcars. Para eso, procederé a tomar todas las variables del dataframe de a una, correlacionándolas con $mpg (millas por galón), verificando si cada una de estas correlaciones resulta mayor a 0.7 (asociación alta), y en ese caso incluyéndola en el análisis de regresión lineal del lado de las variables independientes (y~x1+x2+xN, mtcars).

cor.test(x = mtcars$wt, y = mtcars$mpg, alternative = "two.sided", conf.level = 0.95, method = "pearson") # SI, el valor base !!

## 
##  Pearson's product-moment correlation
## 
## data:  mtcars$wt and mtcars$mpg
## t = -9.559, df = 30, p-value = 1.294e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9338264 -0.7440872
## sample estimates:
##        cor 
## -0.8676594

cor.test(x = mtcars$wt, y = mtcars$cyl, alternative = "two.sided", conf.level = 0.95, method = "pearson") # SI

## 
##  Pearson's product-moment correlation
## 
## data:  mtcars$wt and mtcars$cyl
## t = 6.8833, df = 30, p-value = 1.218e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5965795 0.8887052
## sample estimates:
##       cor 
## 0.7824958

cor.test(x = mtcars$wt, y = mtcars$disp, alternative = "two.sided", conf.level = 0.95, method = "pearson") # SI

## 
##  Pearson's product-moment correlation
## 
## data:  mtcars$wt and mtcars$disp
## t = 10.576, df = 30, p-value = 1.222e-11
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7811586 0.9442902
## sample estimates:
##       cor 
## 0.8879799

cor.test(x = mtcars$wt, y = mtcars$hp, alternative = "two.sided", conf.level = 0.95, method = "pearson") # NO

## 
##  Pearson's product-moment correlation
## 
## data:  mtcars$wt and mtcars$hp
## t = 4.7957, df = 30, p-value = 4.146e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4025113 0.8192573
## sample estimates:
##       cor 
## 0.6587479

cor.test(x = mtcars$wt, y = mtcars$drat, alternative = "two.sided", conf.level = 0.95, method = "pearson") # SI

## 
##  Pearson's product-moment correlation
## 
## data:  mtcars$wt and mtcars$drat
## t = -5.5608, df = 30, p-value = 4.784e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8499795 -0.4839784
## sample estimates:
##        cor 
## -0.7124406

cor.test(x = mtcars$wt, y = mtcars$qsec, alternative = "two.sided", conf.level = 0.95, method = "pearson") # Doble NO (p-value 34%)

## 
##  Pearson's product-moment correlation
## 
## data:  mtcars$wt and mtcars$qsec
## t = -0.97191, df = 30, p-value = 0.3389
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4933536  0.1852649
## sample estimates:
##        cor 
## -0.1747159

cor.test(x = mtcars$wt, y = mtcars$vs, alternative = "two.sided", conf.level = 0.95, method = "pearson") # NO

## 
##  Pearson's product-moment correlation
## 
## data:  mtcars$wt and mtcars$vs
## t = -3.6535, df = 30, p-value = 0.0009798
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7571117 -0.2556982
## sample estimates:
##        cor 
## -0.5549157

cor.test(x = mtcars$wt, y = mtcars$am, alternative = "two.sided", conf.level = 0.95, method = "pearson") # NO

## 
##  Pearson's product-moment correlation
## 
## data:  mtcars$wt and mtcars$am
## t = -5.2576, df = 30, p-value = 1.125e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8386752 -0.4532461
## sample estimates:
##        cor 
## -0.6924953

cor.test(x = mtcars$wt, y = mtcars$gear, alternative = "two.sided", conf.level = 0.95, method = "pearson") # NO

## 
##  Pearson's product-moment correlation
## 
## data:  mtcars$wt and mtcars$gear
## t = -3.9332, df = 30, p-value = 0.0004587
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7744638 -0.2944887
## sample estimates:
##       cor 
## -0.583287

cor.test(x = mtcars$wt, y = mtcars$carb, alternative = "two.sided", conf.level = 0.95, method = "pearson") # NO

## 
##  Pearson's product-moment correlation
## 
## data:  mtcars$wt and mtcars$carb
## t = 2.5909, df = 30, p-value = 0.01464
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.09273981 0.67556998
## sample estimates:
##       cor 
## 0.4276059

De los análisis de correlación se puede observar que:

Variable ($wt vs…)	Observaciones
$mpg	SI, COR: 0.87, P-Value: 1.294*10^-10 (el valor original)
$cyl	SI, COR: 0.78, P-Value: 1.218*10^-7
$disp	SI, COR: 0.89, P-Value: 1.222*10^-11
$hp	NO, COR: 0.66
$drat	SI, COR: 0.71, P-Value: 4.784*10^-6
$qsec	NO, COR: 0.17, (doble NO, con P-value del 34%)
$vs	NO, COR: 0.55
$am	NO, COR: 0.69
$gear	NO, COR: 0.58
$carb	NO, COR: 0.43

Volviendo a executar la regresión lineal…

modelo <- lm(mpg~wt+cyl+disp+drat, mtcars)
summary(modelo)

## 
## Call:
## lm(formula = mpg ~ wt + cyl + disp + drat, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4067 -1.4096 -0.4954  1.3346  6.0729 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 41.160271   7.304325   5.635 5.56e-06 ***
## wt          -3.638075   1.102460  -3.300  0.00272 ** 
## cyl         -1.786074   0.634821  -2.814  0.00903 ** 
## disp         0.007472   0.012062   0.619  0.54080    
## drat        -0.010492   1.337929  -0.008  0.99380    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.642 on 27 degrees of freedom
## Multiple R-squared:  0.8326, Adjusted R-squared:  0.8078 
## F-statistic: 33.57 on 4 and 27 DF,  p-value: 4.057e-10

De los resultados se puede concluir que si bien el mpg sigue estando en función de wt, con el agregado al análisis de cyl, la fortaleza de esta función, medida por el P-value, se han reducido (el indicador de calidad pasó de ‘***’ a ‘**’).

Capítulo 2: Análisis de los Datos

Tal lo dicho en la introducción del presente informe, en esta segunda parte intentaré realizar un análisis estadístico sobre muestras de nivel de pH, unidas estas a datos climáticos.

Será mi hipótesis que existe una relación del tipo lineal entre la variables Temperatura y las mediciones de pH.

Para iniciar con el análisis, procedo a tomar los archivos ‘mediciones del clima.csv’ y ‘mediciones de ph.csv’ generados por los script de los Anexos I & II, ubicándolos en la sub-carpeta ‘recursos’ del presente informe. Estos archivos contienen los datos ya parseados desde sus distintas fuentes. Y luego de esto suprimiendo las columnas que a efecto de este análisis en particular no resultan necesarias.

library("dplyr")              # cargar la librería "dplyr"

# load & filtro de los datos
medicionesClima <- read.csv(file=paste(getwd(), '/recursos/mediciones del clima.csv', sep=''), header=TRUE)
medicionesClima <- select(medicionesClima, -X, -ManianaTemp, -ManianaHumedad, -ManianaPresion, -TardeTemp, -TardeHumedad, -TardePresion)
medicionesClima <- rename(medicionesClima, 'Temp' = 'NocheTemp', 'Humedad' = 'NocheHumedad', Presion = 'NochePresion')

medicionesPh <- read.csv(file=paste(getwd(), '/recursos/mediciones de pH.csv', sep=''), header=TRUE)
medicionesPh <- select(medicionesPh, -X)
medicionesPh <- rename(medicionesPh, 'pH' = 'A.Medir')

A continuación realizo el merge de los datos, prestando especial atención a que en la agrupación de los mismos, se conserven la fecha y el pH medido, agregándole a estos los datos climáticos.

# merge de los datos
medicionesClima <- medicionesClima %>% filter(Fecha %in% medicionesPh$Fecha)
datos <- medicionesPh %>% left_join(medicionesClima, by = 'Fecha')

##         Fecha   pH       Delta    Temp Humedad  Presion
## 1  2019-09-04 1.00 -0.07666667  6.9625  51.625 1026.888
## 2  2019-09-11 5.50  0.05500000 11.8875  85.250 1020.237
## 3  2019-09-25 1.00 -0.16666667 16.3625  52.125 1020.900
## 4  2019-09-25 1.50 -0.05666667 16.3625  52.125 1020.900
## 5  2019-09-25 2.00  0.16666667 16.3625  52.125 1020.900
## 6  2019-09-25 3.00  0.09000000 16.3625  52.125 1020.900
## 7  2019-09-25 3.50  0.60666667 16.3625  52.125 1020.900
## 8  2019-10-02 6.00  0.20000000 12.5625  53.125 1013.962
## 9  2019-10-02 7.00 -0.20500000 12.5625  53.125 1013.962
## 10 2019-10-02 8.00 -0.42500000 12.5625  53.125 1013.962
## 11 2019-10-02 8.50 -0.70000000 12.5625  53.125 1013.962
## 12 2019-10-02 9.00 -0.16000000 12.5625  53.125 1013.962
## 13 2019-10-09 2.66  0.48000000 18.0875  56.250 1016.975
## 14 2019-10-09 4.00  0.13666667 18.0875  56.250 1016.975
## 15 2019-10-09 7.00 -0.04333333 18.0875  56.250 1016.975
## 16 2019-10-09 8.00 -0.04666667 18.0875  56.250 1016.975
## 17 2019-10-16 6.00  0.04000000 12.7625  90.375 1025.575
## 18 2019-10-16 7.00 -0.16000000 12.7625  90.375 1025.575
## 19 2019-10-16 8.00 -0.07333333 12.7625  90.375 1025.575
## 20 2019-10-16 8.50 -0.15333333 12.7625  90.375 1025.575
## 21 2019-10-16 9.00  0.22000000 12.7625  90.375 1025.575
## 22 2019-10-18 1.50  0.23000000 15.1375  91.125 1016.212
## 23 2019-10-23 2.66  0.16000000 14.1625  67.375 1015.263
## 24 2019-10-23 4.00 -0.09000000 14.1625  67.375 1015.263
## 25 2019-10-23 5.50 -0.03500000 14.1625  67.375 1015.263
## 26 2019-10-23 6.00  0.09333333 14.1625  67.375 1015.263
## 27 2019-10-23 7.00 -0.04000000 14.1625  67.375 1015.263
## 28 2019-10-23 8.00 -0.11000000 14.1625  67.375 1015.263
## 29 2019-10-23 8.50  0.15500000 14.1625  67.375 1015.263
## 30 2019-10-23 9.00  0.69000000 14.1625  67.375 1015.263
## 31 2019-10-30 1.00  0.13000000 16.8000  86.750 1010.562
## 32 2019-10-30 2.00  0.16333333 16.8000  86.750 1010.562
## 33 2019-10-30 3.00 -0.18666667 16.8000  86.750 1010.562
## 34 2019-10-30 3.50 -0.65500000 16.8000  86.750 1010.562
## 35 2019-10-30 4.00 -0.32833333 16.8000  86.750 1010.562

Procedo a la normalización de los mismos, y el posterior chequeo y eliminación de los valores atípicos.

# normalización
datos_n <- data.frame('Fecha' = datos$Fecha, 'pH' = datos$pH, 'Delta' = scale(datos$Delta), 'Temp' = scale(datos$Temp), 'Humedad' = scale(datos$Humedad), 'Presion' = scale(datos$Presion))
hist(datos_n$Delta, 10, col=col_caece)  # histograma en deciles

# check y filtro de valores atípicos
par(mar=c(2,2,2,2), mfrow=c(1,2))
boxplot(datos_n, col=col_caece)
# eliminar valores atípicos
datos_n <- filter(datos_n, Delta > -1.7 & Delta < 1.7)
boxplot(datos_n, col=col_caece)

Llegando a este punto, se puede empezar con el contraste de las mediciones de pH, unidas a los datos climáticos.

# inspección visual de correlaciones
pairs(datos_n)

# análisis de correlaciones
cor.test(x = datos_n$Temp, y = datos_n$Delta, alternative = "two.sided", conf.level = 0.95, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  datos_n$Temp and datos_n$Delta
## t = 1.3509, df = 29, p-value = 0.1872
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1214969  0.5502199
## sample estimates:
##       cor 
## 0.2433184

cor.test(x = datos_n$Humedad, y = datos_n$Delta, alternative = "two.sided", conf.level = 0.95, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  datos_n$Humedad and datos_n$Delta
## t = 0.15635, df = 29, p-value = 0.8768
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3286990  0.3794593
## sample estimates:
##       cor 
## 0.0290213

cor.test(x = datos_n$Presion, y = datos_n$Delta, alternative = "two.sided", conf.level = 0.95, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  datos_n$Presion and datos_n$Delta
## t = 0.24753, df = 29, p-value = 0.8062
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3135253  0.3938481
## sample estimates:
##        cor 
## 0.04591587

De los análisis de correlación se puede observar que en ningún caso se encontraron correlaciones significativas entre las variables en juego, ni tampoco se pudo refutar la hipótesis nula en ningún caso.

Variable	Observaciones
$Temp	NO, COR: 0.24, P-Value: 0.1872
$Humedad	NO, COR: 0.02, P-Value: 0.8768
$Presion	NO, COR: 0.04, P-Value: 0.8062

Conclusiones:

Llegado a este punto, no tiene sentido siquiera intentar una regresión lineal como se hizo en el Capítulo 1: investigación.

Dejo como inquietud, que habiéndose tomado las mediciones de pH dentro de la universidad, está mediciones puedieron verse afectadas negativamente por el ambiente controlado de la misma (aire acondicionado), y debido a esto podrían no verse reflejadas asociaciones entre las mediciones de pH y las variables Temperatura y Humedad. Respecto a la presión barométrica podría que esta -en efecto- no influya en las mediciones de pH.

Anexo 1: Script Datos del Clima (SMN)

Los datos del clima han sido tomados todos desde la web del Servicio Meteorológico Nacional, sección “SMN: Datos Abiertos”.

Específicamente, el SMN suministra diariamente un archivo de texto con las diferentes variables climáticas, agrupado por dentro con el detalle correspondiente de cada una se las estaciones meteorológicas disponibles.

Para descargar, por ejemplo, el archivo del día 31 de Octubre de 2019, se debe visitar:
https://ssl.smn.gob.ar/dpd/descarga_opendata.php?file=observaciones/datohorario20191031.txt

FECHA     HORA  TEMP   HUM   PNM    DD    FF     NOMBRE                                             
         [HOA]  [ºC]   [%]  [hPa]  [gr] [km/hr]                                                     
31102019     0  17.5   73  1020.3  140   13     AEROPARQUE AERO                                     
31102019     1  16.8   78  1020.7  140   15     AEROPARQUE AERO                                     
31102019     2  16.9   78  1020.4  140   22     AEROPARQUE AERO                                     
31102019     3  16.4   71  1020.5  140   15     AEROPARQUE AERO                                     
31102019     4  16.1   71  1020.4  140   17     AEROPARQUE AERO                                     
31102019     5  15.4   74  1020.1  140   13     AEROPARQUE AERO                                     
31102019     6  15.5   72  1020.6  140   15     AEROPARQUE AERO                                     
31102019     7  16.3   75  1021.2  140   15     AEROPARQUE AERO                                     
31102019     8  18.3   63  1021.6  110   28     AEROPARQUE AERO                                     
...

Presento a continuación un script en R capaz de parsear estos archivos, extrayendo de estos los datos relevantes aplicados al Capítulo 2: Análisis, del presente informe.

El único requisito para el funcionamiento del script es que este se encuentre en la misma carpeta junto a la colección de archivos de texto a procesar.


# ==================================================
# Script      'Consolidacion Clima.R'
# Propósito   Procesar los datos del clima obtenidios desde el SMN, sección 'Datos Abiertos'.
#             Ubicar este script dentro del directorio con los archivos de input.
# URL         https://www.smn.gob.ar/descarga-de-datos
#             https://ssl.smn.gob.ar/dpd/descarga_opendata.php?file=observaciones/datohorario20191101.txt
#
# Input       Archivos 'datohorarioYYYYMMDD.txt' del SMN.
# Output      Archivo 'mediciones del clima.csv' con los datos del clima consolidados.
# ==================================================

# install.packages("dplyr")   # instalar la librería "dplyr"
library("dplyr")              # cargar la librería "dplyr"

# aumentar el default para que se muestren todos los datos del dataframe
options(max.print=15000)

# setear el working directory a la ubicación actual del script
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))


# ==================================================
# Construir un rango de fechas en formato texto.
# ==================================================
buildDateRange <- function(daysBack = 2) {
  dates <- format(seq(from = as.Date("2019-08-01"), to = Sys.Date()-daysBack, by = "day"), "%Y%m%d")
}


# ==================================================
# Construir nombres de archivos para el rango de fechas recibido.
# ==================================================
buildFileNames <- function(dateRange) {
  basepath = paste(getwd(), '/', sep='')    # ruta dinámina al documento actual
  basename = 'datohorario'
  baseext = '.txt'

  fileNames <- paste(basepath, basename, dateRange, baseext, sep='')
}


# ==================================================
# Leer y combinar el contenido de los nombre de archivos recibidos.
# Se filtra el resultado a los datos de la estación meteorológica Aeroparque.
# ==================================================
getFilesContent <- function(fileNames) {
  filesContent <- NULL
  colWidths <- c(8, 6, 6, 5, 8, 5, 5, 57)
  colNames <- c('FechaOrig', 'Hora', 'Temp', 'Humedad', 'Presion', 'DD', 'FF', 'Estacion')
  colClass <- c('character', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'numeric', 'character')
  
  for (filename in fileNames)
  {
    tryCatch (
    {
      fileContent <- read.fwf(file=filename, skip=2, widths=colWidths, col.names=colNames, colClasses=colClass)
      filesContent <- rbind(filesContent, fileContent) # el merge() no funciona
      # print(fileContent)
    },
    error = function(e) { message(paste('- error:', e)) },
    warning = function(w) { message(paste('- warning:', w)) }
    )
  }

  filesContent
}


# ==================================================
# Eliminar los datos irrelevantes.
# ==================================================
removeUselessData <- function(filesContent) {
  filterString <- "     AEROPARQUE AERO                                     "
  
  filteredContent <- filter(filesContent, Estacion == filterString)
  filteredContent <- mutate(filteredContent, Fecha = as.Date(FechaOrig, "%d%m%Y"))
  filteredContent <- select(filteredContent, -FechaOrig, -DD, -FF, -Estacion)
  filteredContent <- select(filteredContent, Fecha, everything())
  # filteredContent$DD <- NULL
  # print(filteredContent)
  
  filteredContent
}


# ==================================================
# Construir la estructura para los datos del clima
# Noche 0-7hs, Mañana 8-15hs, Tarde 16-23hs
# ==================================================
buildWeatherData <- function(filteredContent) {
  morning <- c(8:15)
  evening <- c(16:23)
  night <- c(0:7)
  
  dfNight <- filteredContent %>%
    filter(Hora %in% night) %>%
    group_by(Fecha) %>%
    summarize(NocheTemp = mean(Temp), NocheHumedad = mean(Humedad), NochePresion = mean(Presion))

  dfMorning <- filteredContent %>%
    filter(Hora %in% morning) %>%
    group_by(Fecha) %>%
    summarize(MañanaTemp = mean(Temp), MañanaHumedad = mean(Humedad), MañanaPresion = mean(Presion))

  dfEvening <- filteredContent %>%
    filter(Hora %in% evening) %>%
    group_by(Fecha) %>%
    summarize(TardeTemp = mean(Temp), TardeHumedad = mean(Humedad), TardePresion = mean(Presion))

  weatherData <- data.frame(dfNight, dfMorning, dfEvening)
  weatherData <- select(weatherData, -Fecha.1, -Fecha.2)
  weatherData <- arrange(weatherData, Fecha)
  
  weatherData
}


# ==================================================
# Grabar un archivo con los datos del clima
# ==================================================
saveWeatherData <- function(weatherData) {
  
  filename = 'mediciones del clima.csv'
  fullname = paste(getwd(), '/', filename, sep='')    # ruta dinámina al documento actual
  
  write.csv(weatherData, file=fullname)
}


# ==================================================
# MAIN
# ==================================================
print("Step #01: buildDateRange() ...")
dateRange <- buildDateRange(1)
# print(dateRange)

print("Step #02: buildFileNames() ...")
fileNames <- buildFileNames(dateRange)
# print(fileNames)

print("Step #03: getFilesContent() ...")
filesContent <- getFilesContent(fileNames)
# print(head(filesContent, 200))

print("Step #04: removeUselessData() ...")
filteredContent <- removeUselessData(filesContent)
# print(head(filteredContent, 200))

print("Step #05: buildWeatherData() ...")
weatherData <- buildWeatherData(filteredContent)
# print(weatherData)

print("Step #06: saveWeatherData() ...")
saveWeatherData(weatherData)

print('DONE!')

Anexo 2: Script Datos de pH

Script para procesar los datos en crudo con las mediciones de pH, desde el excel suministrado por el grupo de trabajo.

A.Medir  Calibracion.1  Calibracion.2  Real   Tiempo  Intentos  Fecha
   0.00              1              2  0.00       NA        NA  04/09
   0.00              1              2  0.00       NA        NA  04/09
   0.00              3              1  0.10  8' 30''        NA  04/09
   0.00              3              2  0.80  3' 39''        NA  04/09
   0.00              1              4  0.49     30''         1  30/10
   0.00              2              3  0.45  3' 40''         1  30/10
   0.00              2              4  1.14     50''         1  30/10
   0.00              1              3  0.54     30''         1  30/10
   0.00              2              1  0.17     40''         1  30/10
...


# ==================================================
# Script      'Consolidacion pH.R'
# Propósito   Procesar las mediciones de pH según el template recibido.
#             Ubicar este script dentro del directorio con el archivo de input.
#
# Input       Archivo 'Mediciones.csv' exportado desde su archivo excel homónimo.
# Output      Archivo 'mediciones de ph.csv' con los datos del clima consolidados.
# ==================================================


# install.packages("dplyr")   # instalar la librería "dplyr"
library("dplyr")              # cargar la librería "dplyr"

# aumentar el default para que se muestren todos los datos del dataframe
options(max.print=15000)

# setear el working directory a la ubicación actual del script
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))


# ==================================================
# Construir el nombre de archivo a leer.
# ==================================================
buildFileName <- function() {
  
  filename = 'Mediciones.csv'
  fullname = paste(getwd(), '/', filename, sep='')    # ruta dinámina al documento actual
}


# ==================================================
# Leer el contenido del nombre de archivo recibido.
# ==================================================
getFileContent <- function(fileName) {
  
  colNames <- c('A Medir', 'Calibracion 1', 'Calibracion 2', 'Real', 'Tiempo', 'Intentos', 'FechaOrig')
  colClass <- c('numeric', 'numeric', 'numeric', 'numeric', 'character', 'numeric', 'character')
  
  tryCatch (
    {
      fileContent <- read.csv(file=fileName, col.names=colNames, colClasses=colClass)
    },
    error = function(e) { message(paste('- error:', e)) },
    warning = function(w) { message(paste('- warning:', w)) }
  )

  fileContent
}


# ==================================================
# Eliminar los datos irrelevantes.
# ==================================================
removeUselessData <- function(fileContent) {
  
  filteredContent <- filter(fileContent, is.na(A.Medir) == FALSE)
  filteredContent <- filter(fileContent, A.Medir != 0)
  filteredContent <- mutate(filteredContent, Fecha = as.Date(FechaOrig, "%d/%m"))
  filteredContent <- select(filteredContent, -FechaOrig, -Tiempo, -Intentos, -Calibracion.1, -Calibracion.2)
  filteredContent <- mutate(filteredContent, Delta = A.Medir - Real)
  filteredContent <- select(filteredContent, Fecha, everything())

  filteredContent
}


# ==================================================
# Construir la estructura para los datos de pH
# ==================================================
buildPhData <- function(filteredContent) {

  phData <- filteredContent %>%
    group_by(Fecha, A.Medir) %>%
    summarize(Delta = mean(Delta))
  phData <- arrange(phData, Fecha)

  phData
}


# ==================================================
# Grabar un archivo con los datos de pH
# ==================================================
savePhData <- function(phData) {
  
  filename = 'mediciones de ph.csv'
  fullname = paste(getwd(), '/', filename, sep='')    # ruta dinámina al documento actual

  write.csv(phData, file=fullname)
}


# ==================================================
# MAIN
# ==================================================
print("Step #01: buildFileName() ...")
fileName <- buildFileName()

print("Step #02: getFileContent() ...")
fileContent <- getFileContent(fileName)

print("Step #03: removeUselessData() ...")
filteredContent <- removeUselessData(fileContent)

print("Step #04: buildPhData() ...")
phData <- buildPhData(filteredContent)

print("Step #05: savePhData() ...")
savePhData(phData)

print('DONE!')

Anexo 3: Proceso de Medición de pH

El pH (Potencial de Hidrógeno) es una medida que la actividad del hidrógeno en una solución dada.

Lo que un Ph-metro hace es medir con un electrodo (punta de vidrio) la diferencia de potencial entre la solución medida y otro electrodo de referencia (asumo que interno en el aparato) conectando este último a plata o cloruro de plata: materiales de los cuales se conoce su estado natural de ionización y que se ven relativamente poco afectados por las cambios de temperatura ambiente.

Midiendo la conductividad del medio (agua o alguna otra solución del tipo acuosa) se puede deducir su pH.
Valores típicos de conductividad se expresan como mili o micro Siemens/cm (mS/cm, uS/cm respectivamente).

Se conoce como soluciones buffer a mezclas para contrarrestar un nivel de pH elevado o bajo en una solución.
Un uso común de esto es en piscinas, donde el cloro solo actúa si el pH del agua se encuentra entre 6.5 y 8.

¿Cómo funciona un electrodo de pH?

“El electrodo envía un voltaje al medidor que se genera a partir del pH de la muestra en la que se encuentre sumergido.”

"Limpiar el vidrio de detección del electrodo de pH puede producir una carga estática (piense en el frotar un globo en la superficie y la carga estática que se acumula en él).
“La carga estática interfiere con la lectura de voltaje del electrodo.”
“Cuando la lectura de voltaje es incorrecta, el valor de pH que interpreta el voltaje también.”

“Es más, la capa de vidrio hidratada interrumpida al limpiar el bulbo con una toalla de papel. En lugar de limpiar el bulbo del electrodo, simplemente enjuague el electrodo con agua destilada o desionizada.”
[Fuente: #7]

Aspectos técnicos y calibración

"Todos los electrodos de pH se basan en el principio conocido como la ecuación de Nernst.

La ecuación de Nernst toma una lectura de voltaje (mV) y la correlaciona con la concentración del ion (o pH). Esta correlación es una línea recta. Para los electrodos de pH, el valor teórico en mV a pH 7 es 0 mV (neutral) y la pendiente de la línea es 59.16 mV."

"Esto significa que en teoría, el electrodo va a cambiar su salida en 59.16 mV por cada unidad de pH (por ejemplo pH 6 a pH 7 sería 59.16 mV / unidad de pH). Todo esto teóricamente, dado que los electrodos van a cambiar su pendiente y offset a medida que envejecen."*

"En realidad, el electrodo podría funcionar ligeramente diferente al comportamiento teórico (por ejemplo 58,2 mV pendiente y 8mV offset).
La calibración compensa esto mediante la determinación de la pendiente real y el offset de su electrodo usando soluciones buffer conocidas y actualizando el algoritmo en el medidor como corresponde."

"Para obtener los mejores resultados, usted debe asegurarse de que está calibrando con soluciones buffer ‘adecuadas’ para su muestra.
La solución buffer de pH 7 siempre debe ser incluida para obtener el punto de desplazamiento (neutro). Esto significa que si la muestra está alrededor de pH 8.6, se deben utilizar soluciones buffer pH 7 y pH 10."
[Fuente: #7]

FUENTES CONSULTADAS:

[#1 PH]
https://es.wikipedia.org/wiki/PH

[#2 Solucion Patron]
https://es.wikipedia.org/wiki/Solución_patrón

[#3 Potencial normal de electrodo]
https://es.wikipedia.org/wiki/Potencial_normal_de_electrodo

[#4 Electrodo de referencia]
https://es.wikipedia.org/wiki/Electrodo_de_referencia

[#5 Autoionización del agua]
https://es.wikipedia.org/wiki/Autoionización_del_agua

[#6 Conductividad eléctrica (CE)]
https://www.infoagro.com/instrumentos_medida/doc_conductividad_electrica.asp?k=53

[#7 ¿Está cometiendo errores en la medición de pH?]
https://www.hannacolombia.com/blog/post/88/esta-cometiendo-errores-en-la-medicion-ph

[#8 Indicador de pH]
https://es.wikipedia.org/wiki/Indicador_de_pH
https://es.wikipedia.org/wiki/Anexo:Indicadores_de_pH

[#9 Valoración Acido-Base]
https://es.wikipedia.org/wiki/Valoración_ácido-base

[#10 Soluciones Patrón Comerciales]
http://www.merckmillipore.com/AR/es/products/analytics-sample-prep/reference-materials/certipur-reference-material-for-reliable-calibration-in-pH-measurement/ready-to-use-buffer-solutions/PH2b.qB.T0IAAAE_d.x3.Lxi,nav
http://www.crisoninstruments.com/es/laboratorio/disoluciones-/patrones-conductividad/disolucion-patron-de-conductividad-147-s-cm

[#11 Ecuación de Nernst]
https://es.wikipedia.org/wiki/Ecuación_de_Nernst

[#12 Reducción-oxidación (Redox)]
https://es.wikipedia.org/wiki/Reducción-oxidación

Historial de Versiones

Versión #01

NEW: investigar que es markdown.
NEW: documentar en markdown con código fuente + resultados.
mean (mu), desviación estandar (sigma), pairs, hist.
NEW: investigar str(mtcars) + summary(mtcars).
NEW: investigar head(mtcars) TOP 6.
NEW: investigar tail(mtcars) LAST 6.
NEW: investigar detección de outliers.
NEW: investigar grafico de boxplot.

Versión #02

NEW: indicador visual para diferenciar si un sección corresponde a ‘sintaxis de R’, o un ‘concepto estadístico’.
NEW: uso de color verde “CAECE” a las gráficas.
NEW: sección Markdown con referencia a su sintaxis.
NEW: sección Tipos de Objetos (arreglo, matriz, factor, lista, etc).
FIX: comentarios en las cajas de código fuente.
FIX: cambié las definiciones de conceptos en palabras propias por definiciones formales.
NEW: enlace a las fuentes en las definiciones de conceptos.
FIX: mayor detalle en la sección Vectores.
FIX: en la sección Promedio aclaro que este aplica sobre variables.
NEW: gráfico de boxplot al mtcars completo y descripción como intro al concepto de Normalización.
NEW: sección Normalización: “a mano”.
NEW: sección Normalización, parte II: uso de scale(mtcars).
NEW: mini-sección: Importación de Datos.

Versión #03

NEW: sección desnormalización de datos.
NEW: sección correlación lineal, parte I - teoría.
NEW: sección significancia.
NEW: sección correlación lineal, parte II - práctica.
NEW: sección regresión lineal, parte I - teoría.
NEW: sección regresión lineal, parte II - práctica.

Versión #04

NEW: búsqueda de fuentes de datos sobre mediciones climáticas (humedad, temperatura, presion ambiental)
fuente elegida: Serv. Meteorológico Nacional, https://www.smn.gob.ar/descarga-de-datos
FIX: sección regresión lineal, parte II - práctica.
NEW: sección validación del modelo (residuos & QQ-Normal).
NEW: sección regresión lineal, parte III - práctica.

Versión #05

FIX: separación en capítulos, estructura final del informe.
NEW: capítulo anexo 1, con script para consolidación de los datos del clima del Servicio Metereológico Nacional.

Versión #06

NEW: capítulo anexo 2, con script para consolidación de las medición del pH.
NEW: capítulo anexo 3, con investigación sobre el proceso de medición del pH y cómo funciona un pH-metro.

Versión #07

NEW: capítulo 2, análisis de los datos.