Utilizando el Software R
R?R es un lenguaje y un ambiente para el manejo de datos, cálculos, y gráficos en código libre.
Además de ser gratuito, los desarrollos en R se actualzan más rápido que cualquier otro de los costosos software comerciales que se encuentran en el mercado.
Ventajas de usar R
… Es software libre y por tanto su costo es nulo.
… Es multiplataforma: existen versiones para LinuX, Mac y Windows. Los procedimientos y análisis
… Capacidad para acceder a datos en múltiples formatos. Dispone de librerías para leer datos desdes SPSS, SAS, Access, MySQL, Excel, etc.
… Generación de gráficos de alta calidad.
… Todos los algoritmos implementados en R pueden ser vistos e interpretados por cualquier usuario, por lo que este puede saber exactamente que es lo que hace el ordenador cuando ejecuta un comando.
Desventajas de Usar R
Hay empresas que por políticas no pueden instalar software libre en sus maquinas cada una tiene su politica, sus software de preferencia, sus necesidades, etc.
Una de las principales desventajas es que hasta hace poco el uso de R estaba limitado a entornos universitarios y de usuarios con gran conocimiento de la estadística y la programación.
No hay nadie a quien reclamar si algo falla, ni hay un departamento de atención al cliente que nos diga qué podemos hacer si algo va mal, si alguién procedimiento nos da un error, o simplemente si no sabemos qué sintaxis utilizar.
RPara realizar la instalación de R y RStudio en windows, Mac, Ubuntu o Linux se debe ingresar a los siguientes sitios web:
| Operador Matemático | Función en R |
|---|---|
| \(\sqrt{x}\) | sqrt(x) |
| \(e^x\) | exp(x) |
| \(x!\) | factorial(x) |
| \(\log(x)\) | log(x) |
| \(\sin(x)\) | sin(x) |
| \(\cos(x)\) | cos(x) |
| \(\tan(x)\) | tan(x) |
| \(\pi\) | pi |
<- asignar valores a una variable.Un vector es una secuencia ordenada de datos, los cuales han de ser del mismo tipo.
Tipos de datos que se pueden almacenar en un vector se destacan los siguientes:
logical (TRUE o FALSE)
integer (números enteros)
numeric (números reales)
character (palabras)
[1] "Fundamentos" "Estadística" "en" "R"
[1] 0.000000 1.000000 1.414214 1.732051 2.000000 2.236068 2.449490 2.645751
[9] 2.828427 3.000000 3.162278
max() y min() calculan sus valores máximos y mínimos respectivamente.
sum() calcula la suma
prod() calcula el producto
mean() calcula la media
sort() ordena los elementos del vector en el orden natural creciente del tipo de datos que lo forman, se puede incluir en su argumento el parámetro decreasing = TRUE.
Realizando operaciones con estas nuevas funciones se tiene:
Una matriz no es más que un arreglo de números en m renglones y n columnas.
La función matrix(vector, nrow = n,) define una matriz de n filas.
La función
rbind(vector1, vector2, vector3, . . . )
permite la construcción de la matriz de filas
cbind(): obtiene la matriz cuyas columnas son los vectores a los que se aplica.
[,1] [,2] [,3]
[1,] 1 2 1
[2,] 0 3 2
[3,] 2 6 0
Un dataframe tiene la apariencia de una matriz, pero con la diferencia que cada columna de un dataframe tiene la apariencia de una matriz. Particularmente cada columna de un dataframe puede contener datos de un tipo diferente.
Se usa la función data.frame(v1, v2, . . .) aplicada a los vectores en el orden en el que queramos disponer las columnas de la tabla.
# Número de polizas
Poliza <- c(1:9)
# Tipo de poliza
Tipo <- c("Hogar", "Auto", "Auto", "Auto", "Hogar", "Hogar", "Auto",
"Auto", "Hogar")
# Fecha de creación de la póliza
Fecha <- c("12/12/2016", "08/02/2014", "10/08/2012", "01/01/2015",
"21/11/2011", "18/01/2016", "12/04/2005", "29/03/2007",
"18/02/2009")
# Coberturas incluidas en la poliza
Coberturas <- c("Incendio", "Robo", "Terceros", "Robo", "Robo",
"Incendio", "Terceros", "R. Civil", "Incendio")
# Datos de los titulares de las polizas
Nombre <- c("Carlos", "Nancy", "Pedro", "Cecilia", "Ricardo", "Sofia",
"Armando", "Vicente", "Fernando")
# Genero, sexo, edad, estado de providencia
Sexo <- c("M", "F", "M", "F", "M", "F", "M", "M", "M")
Edad <- c(25, 35, 45, 47, 24, 43, 33, 31, 40)
Estado <- c("Campeche", "Chiapas", "Ciudad de M ́exico", "Coahuila",
"Durango", "Guanajuato", "Guerrero", "Hidalgo", "Jalisco") Poliza Tipo Fecha Coberturas Nombre Sexo Edad
1 1 Hogar 12/12/2016 Incendio Carlos M 25
2 2 Auto 08/02/2014 Robo Nancy F 35
3 3 Auto 10/08/2012 Terceros Pedro M 45
4 4 Auto 01/01/2015 Robo Cecilia F 47
5 5 Hogar 21/11/2011 Robo Ricardo M 24
6 6 Hogar 18/01/2016 Incendio Sofia F 43
7 7 Auto 12/04/2005 Terceros Armando M 33
8 8 Auto 29/03/2007 R. Civil Vicente M 31
9 9 Hogar 18/02/2009 Incendio Fernando M 40
installed.packages("name") \(\Rightarrow\) Instalador paquetes
library(name) \(\Rightarrow\) Llamar o cargar el paquete
Primero cargamos la librería tidyverse
Importar datos tipo .csv
Exportar datos como .csv
tidyverse no es la única librería que nos permite exportar e importar conjuntos de datos.
library(xlsx) \(\Rightarrow\) R a Excellibrary(haven) \(\Rightarrow\) R a SPSSlibrary(foreign) \(\Rightarrow\) R a STATARDatos sobre los hábitos de sueño de diferentes mamíferos.
Estan contenidos por defecto en R en el paquete de ggplot2.
# A tibble: 6 × 11
name genus vore order conservation sleep_total sleep_rem sleep_cycle awake
<chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 Cheetah Acin… carni Carn… lc 12.1 NA NA 11.9
2 Owl mo… Aotus omni Prim… <NA> 17 1.8 NA 7
3 Mounta… Aplo… herbi Rode… nt 14.4 2.4 NA 9.6
4 Greate… Blar… omni Sori… lc 14.9 2.3 0.133 9.1
5 Cow Bos herbi Arti… domesticated 4 0.7 0.667 20
6 Three-… Brad… herbi Pilo… <NA> 14.4 2.2 0.767 9.6
# ℹ 2 more variables: brainwt <dbl>, bodywt <dbl>
| Columna | Descripción |
|---|---|
| name | nombre común |
| genus | rango taxonómico |
| vore | carnívoro, omnívoro o herbívoro |
| order | Estado de conservación del mamífero |
| sleep_total | cantidad total de sueño, en horas |
| sleep_rem | sueño rem, en horas |
| sleep_cycle | duración del ciclo de sueño |
| awake | cantidad de tiempo que pasa despierto en horas |
| brainwt | peso del cerebro en kilogramos |
| bodywt | peso corporal en kilogramos |
Antes de sumergirnos en esto, recordemos cómo funcionan los histogramas.
| Estadística | Función en R |
|---|---|
| Media | mean() |
| Mediana | median() |
| Moda | mod() |
Calculemos estas estadísticas en R para el número total de horas de sueño de los mamíferos de nuestro dataset.
# A tibble: 65 × 2
sleep_total n
<dbl> <int>
1 12.5 4
2 10.1 3
3 5.3 2
4 6.3 2
5 8.4 2
6 8.7 2
7 9.1 2
8 9.4 2
9 9.8 2
10 10.3 2
# ℹ 55 more rows
Note que existen 4 mamíferos en nuestro dataset que duermen durante 12.5 horas.
Para responder la interrogante, consideremos el caso que tenemos todos los insectivoros de nuestra data:
library(tidyverse)
msleep %>%
filter(vore == "insecti") %>%
summarize(mean_sleep = mean(sleep_total),
median_sleep = median(sleep_total))# A tibble: 1 × 2
mean_sleep median_sleep
<dbl> <dbl>
1 14.9 18.1
Considere ahora que hemos descubierto un nuevo insectívoro que misteriosamente nunca duerme. Con este supuesto suponga que nuestro resultados son ahora
\[ media = 11.5 \hspace{1cm}y\hspace{1cm} mediana = 17.6 \]
La media se puede reducir en más de 3 horas, mientras que la mediana cambió en menos de una hora. Esto se debe a que la media es mucho más sensible a los valores extremos que la mediana. Esto en datos simétricos
Si los datos estan sesgados, lo que significa que no son simétricos como se muestran a continuación.
Cuando los datos están sesgados, la media y la mediana son diferentes. La media se tira en dirección del sesgo, por lo que es más baja que la mediana en los datos sesgados a la izquierda y más alta que la media en los datos sesgados a la derecha.
Dataset food_consumption indice de carbono alimentario 2018 en kilogramos de alimentos consumidor por persona por año en cada país en cada categoría de alimento (consumption).
co2_emissions información sobre la huella de carbino de alguna categoría de alimento consumption medida en kilogramos de dióxido de carbono \(CO_2\), por persona en cada país.
Encontremos medidas de tendencia central para comparar el consumo de alimentos en EE.UU y Bélgica.
# Cargamos el conjunto de datos en R
food_consumption <- readRDS("food_consumption.rds")
head(food_consumption)# A tibble: 6 × 4
country food_category consumption co2_emission
<chr> <fct> <dbl> <dbl>
1 Argentina pork 10.5 37.2
2 Argentina poultry 38.7 41.5
3 Argentina beef 55.5 1712
4 Argentina lamb_goat 1.56 54.6
5 Argentina fish 4.36 6.96
6 Argentina eggs 11.4 10.5
# Creemos dos marcos de datos:
# 1. uno que contenga las filas de consumo de Bélgica
consumo_belgica <- food_consumption %>%
filter(country == "Belgium")
head(consumo_belgica, 5)# A tibble: 5 × 4
country food_category consumption co2_emission
<chr> <fct> <dbl> <dbl>
1 Belgium pork 38.6 137.
2 Belgium poultry 12.2 13.1
3 Belgium beef 15.6 482.
4 Belgium lamb_goat 1.32 46.2
5 Belgium fish 19.0 30.3
# 2. Otro que contenga las filas de consumo para USA
consumo_usa <- food_consumption %>%
filter(country == "USA")
consumo_usa# A tibble: 11 × 4
country food_category consumption co2_emission
<chr> <fct> <dbl> <dbl>
1 USA pork 27.6 97.8
2 USA poultry 50.0 53.7
3 USA beef 36.2 1118.
4 USA lamb_goat 0.43 15.1
5 USA fish 12.4 19.7
6 USA eggs 14.6 13.4
7 USA dairy 255. 363.
8 USA wheat 80.4 15.3
9 USA rice 6.88 8.8
10 USA soybeans 0.04 0.02
11 USA nuts 7.86 13.9
Cálculo de la media y mediana de kilogramos de alimentos consumidos por persona cada año año para ambos países.
food_consumption %>%
# Filtramos para Bélgica y EE.UU
filter(country %in% c("Belgium", "USA")) %>%
# Agrupamos por país
group_by(country) %>%
# Obtenemos resumen de consumo meduo y mediano
summarise(mean_consumption = mean(consumption),
median_consumption = median(consumption))# A tibble: 2 × 3
country mean_consumption median_consumption
<chr> <dbl> <dbl>
1 Belgium 42.1 12.6
2 USA 44.6 14.6
co2_emission para el arroz (rice).Respondamos las siguientes preguntas:
¿Observe el histograma para el arroz (rice) que acabamos de calcular, cuál de los siguientes términos describe mejor los datos?
Sin sesgo
Sesgado a la izquierda
Sesgado a la derecha
¿Qué medida de tendencia central resume mejor los kilogramos de \(CO_2\) para el arroz?
Media
Mediana
Ambos la media y la mediana
La varianza
Mide la distancia promedio de cada punto de datos hasta la media de los datos.
Matemáticamente esto es
\[ s^2 = \frac{1}{n-1}\sum_{i=1}^n (y_i - \bar{y})^2 \Longrightarrow \mbox{Varianza muestral} \]
Desviación Estándar
Es una medida que se utiliza para cuantificar la variación o dispersión de un conjunto de datos.
\[ \sigma = \sqrt{s^2} \Longrightarrow \mbox{Desviación estándar de la muestra} \]
Un \(\sigma\) bajo indica que la mayor parte de los datos de una muestra tienden a estar agrupados cerca de su media.
Si \(\sigma\) es alto indica que los datos se extiende sobre un rango de valores más amplio.
Cuartiles
Divide los datos en cuatro partes iguales
| Estadística o Métrica | Función en R |
|---|---|
| \(s^2\) | var() |
| \(\sigma\) | sd() |
| cuartiles | quantile() |
Veamos como hacer esto en R.
\[P(Evento) = \frac{\mbox{# de formas en que puede ocurrir el evento}}{\mbox{# total de resultados posibles}}\]
La probabilidad siempre está entre cero y 100%.
\[ P(Brian) = \frac{1}{4} = \fbox{25%} \]
Podemos recrear el escenario anterior de la selección de las cuatro personas, utilizando la función sample_n() de dplyr.
sample_n() toma de un marco de datos la cantidad de filas que queremos extrer de manera aleatoria.\[ P(Brian) = \frac{1}{4} = \fbox{25%} \]
\[ P(Claire) = \frac{1}{3} = \fbox{33%} \]
Recreemos nuestros ejemplos en R,
Probabilidades en R
# Datos de venta de Brian
ventas_Brian <- data.frame(product = c("Product A", "Product B", "Product C",
"Product D", "Product E", "Product F",
"Product G", "Product H", "Product I",
"Product J", "Product N"),
n = c(23, 62, 15, 40, 5, 11, 2, 8, 7, 2, 3))
# Creemos una nueva columna llamada prob
ventas_Brian %>%
mutate(prob = n/sum(n)) product n prob
1 Product A 23 0.12921348
2 Product B 62 0.34831461
3 Product C 15 0.08426966
4 Product D 40 0.22471910
5 Product E 5 0.02808989
6 Product F 11 0.06179775
7 Product G 2 0.01123596
8 Product H 8 0.04494382
9 Product I 7 0.03932584
10 Product J 2 0.01123596
11 Product N 3 0.01685393
Se utilizan para modelar variables contables, como el lanzamiento de un dado de seis caras.
Describe la probabilidad de cada resultado posible en un escenario
Es la media de una distribución, esto es
\[E(Y) = \sum_{y}yp(y)\]
Tomemos una muestra de 10 lanzamientos
n
1 3
2 3
3 3
4 4
5 3
6 6
7 5
8 2
9 3
10 6
Veamos la media de la muestra
Una muestra más grande
Visualizamos la muestra
Lo anterior se le conoce como ley de los grandes números, ya que a medida que aumenta el tamaño de la muestra la media de la muestra se acerca a la media teórica.
Describe la probabilidad del número de éxitos en una secuencia de ensayos independientes.
Se dice que una variable aleatoria \(Y\) tiene distribución binomial basada en \(n\) pruebas con probabilidad de éxito p si y solo si
\[ p(Y = y) = \binom{n}{y}p^yq^{n-y}, \hspace{1.2cm} y = 0, 1,2, . . . n \hspace{0.5cm} y \hspace{0.5cm} 0\leq p \leq 1 \]
con
\[ E[Y] = n\times p \]
\[ V[Y] = np(1- p) \]
| Probabilidades o simulaciones | Función en R |
|---|---|
| valores aleatorios de probabilidad de una muestra dada | rbinom(n, muestra, p) |
| \(P(Y \leq y)\) | dbinom(n, muestra, p) |
| \(P(Y> y)\) | pbinom(n, muestra, p, lower.tail = F) |
Es una de las distribuciones de probabilidad más importantes que aprenderá, ya que innumerables métodos estadísticos se basan en ella y se aplica a más situaciones del mundo real que otras distribuciones.
Se dice que una variable aleatoria \(Y\) tiene una distribución normal de probabilidad si y sólo si para \(\sigma > 0\) y \(-\infty < \mu < \infty\),
\[f(y) = \frac{1}{\sigma\sqrt{2\pi}}e^{-(y-\mu)^2/(2\sigma^2)}, \hspace{1cm}-\infty < y < \infty\]
con
\[ E(Y) = \mu \hspace{1cm} y \hspace{1cm} V(Y) = \sigma^2 \]
Gráficamente
\(P(Y \leq y)\): pnorm(y, mean, sd)
\(P(Y > y)\): pnorm(y, mean, sd, lower.tail = F)
\(P(y_1 \leq Y \leq y_2)\): pnorm(y_2, mean, sd)-pnorm(y_1, mean, sd)
Considere la data de la Encuesta Nacional de Examen de Salud y Nutrición (NHANES)
Paquetes
Hipótesis
\(H_0:\) La muestra proviene de una distribución normal
\(H_1:\) La muestra no proviene de una distribución normal
Nivel de Significancia
\(\alpha = 0.05\)
Si \(p < \alpha\) se rechaza \(H_0\)
Si \(p \geq \alpha\) no se rechaza \(H_0\)
¿Nuestros datos provienen de una distribución normal?
Anderson-Darling normality test
data: NHANES$Height
A = 452.86, p-value < 2.2e-16
¿Cuál es la altura media y la desviación estándar?
# A tibble: 1 × 2
media sd
<dbl> <dbl>
1 162. 20.2
¿Qué porcentaje de mujeres mide menos de 154 cm?
¿Qué porcentaje de mujeres mide más de 154 cm?
¿Qué porcentaje de mujeres mide entre 154 - 157 cm?
La correlación muestra la relación entre las variables. Se puede visualizar utilizando diagramas de dispersión.
Eje x: Variable explicativa o independiente
Eje y: Variable de respuesta o dependiente
# A tibble: 6 × 7
Date Open High Low Close `Adj Close` Volume
<date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 2009-01-02 3.07 3.25 3.04 3.24 2.75 746015200
2 2009-01-05 3.33 3.44 3.31 3.38 2.87 1181608400
3 2009-01-06 3.43 3.47 3.30 3.32 2.82 1289310400
4 2009-01-07 3.28 3.30 3.22 3.25 2.76 753048800
5 2009-01-08 3.23 3.33 3.22 3.31 2.81 673500800
6 2009-01-09 3.33 3.34 3.22 3.24 2.75 546845600
Ambas variables están correlacionadas fuertemente.
La correlación entre el peso corporal y el tiempo de vigilia es solo del 0.3 , que es una relación lineal débil.
Call:
lm(formula = grasas ~ edad, data = datos)
Residuals:
Min 1Q Median 3Q Max
-63.478 -26.816 -3.854 28.315 90.881
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 102.5751 29.6376 3.461 0.00212 **
edad 5.3207 0.7243 7.346 1.79e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 43.46 on 23 degrees of freedom
Multiple R-squared: 0.7012, Adjusted R-squared: 0.6882
F-statistic: 53.96 on 1 and 23 DF, p-value: 1.794e-07
\[ grasa = 102.58 + 5.32edad \]
El \(R^2\) o Multiple R-Squared: 0.701 mide la bondad del ajuste de la recta a los datos.
Suponga que quiere predecir la cantidad de grasa para individuos de edades 30, 31, 32, . . . , 50.
1 2 3 4 5 6 7 8
262.1954 267.5161 272.8368 278.1575 283.4781 288.7988 294.1195 299.4402
9 10 11 12 13 14 15 16
304.7608 310.0815 315.4022 320.7229 326.0435 331.3642 336.6849 342.0056
17 18 19 20 21
347.3263 352.6469 357.9676 363.2883 368.6090
Usaremos el conjunto de datos marketing [paquete darium] Contiene el impacto de tres medios publicitarios (youtube, facebook y paper) en las ventas. Los datos son el presupuesto publicitario en miles de dólares junto con las ventas.
youtube facebook newspaper sales
1 276.12 45.36 83.04 26.52
2 53.40 47.16 54.12 12.48
3 20.64 55.08 83.16 11.16
4 181.80 49.56 70.20 22.20
Queremos predicir las ventas futuras sobre la base del presupuesto publicitario gastado en YouTube.
Representación de la recta de mínimos cuadrados
Evaluación del Modelo
Call:
lm(formula = sales ~ youtube, data = marketing)
Residuals:
Min 1Q Median 3Q Max
-10.0632 -2.3454 -0.2295 2.4805 8.6548
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.439112 0.549412 15.36 <2e-16 ***
youtube 0.047537 0.002691 17.67 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.91 on 198 degrees of freedom
Multiple R-squared: 0.6119, Adjusted R-squared: 0.6099
F-statistic: 312.1 on 1 and 198 DF, p-value: < 2.2e-16
\[ ventas = 8.439112 + 0.047537\cdot youtube \]
Pronósticos
Una serie de tiempo es una colección de observaciones de una variable tomadas de forma secuencial y ordenada en el tiempo.
\[ y_t = \{y_1, y_2, . . . . , y_T\} \]
El análisis de series de tiempo se utiliza para muchas aplicaciones:
Previsión Económica
Pronóstico de ventas
Análisis presupuestario
Análisis del mercado de valores
Proyecciones de rendimiento
Proceso y control de calidad
Estudios de inventario
Ejemplo
El objetivo es el análisis, modelado y predicción de la serie temporal \(CO_2\)
# Cargamos los datos co2
co2<-read.csv("co2.csv", header = TRUE, sep = ";")
# Transformamos los datos en una serie temporal
co2ts<-ts(co2$CO2, start = c(1959,1), frequency = 12)
print(co2ts) Jan Feb Mar Apr May Jun Jul Aug Sep Oct
1959 333.13 332.09 331.10 329.14 327.36 327.29 328.23 329.55 330.62 331.40
1960 333.92 333.43 331.85 330.01 328.51 328.41 329.25 330.97 331.60 332.60
1961 334.68 334.17 332.96 330.80 328.98 328.57 330.20 331.58 332.67 333.17
1962 336.82 336.12 334.81 332.56 331.30 331.22 332.37 333.49 334.71 335.23
1963 337.95 338.00 336.37 334.47 332.46 332.29 333.76 334.80 336.00 336.63
1964 339.05 339.27 337.64 335.68 333.77 334.09 335.29 336.76 337.77 338.26
1965 341.47 341.31 339.41 337.74 336.07 336.07 337.22 338.38 339.32 340.41
1966 343.02 342.54 340.88 338.75 337.05 337.13 338.45 339.85 340.90 341.70
1967 344.28 343.42 342.02 339.97 337.84 338.00 339.20 340.63 341.41 342.68
1968 345.92 345.40 344.16 342.11 340.11 340.15 341.38 343.02 343.87 344.59
1969 347.38 346.78 344.96 342.71 340.86 341.13 342.84 344.32 344.88 345.62
1970 348.53 347.87 346.00 343.86 342.55 342.57 344.11 345.49 346.04 346.70
1971 349.93 349.26 347.44 345.55 344.21 343.67 345.09 346.27 347.33 347.82
1972 351.71 350.94 349.10 346.77 345.73
Nov Dec
1959 331.87 333.18
1960 333.57 334.72
1961 334.86 336.07
1962 336.54 337.79
1963 337.93 338.95
1964 340.10 340.88
1965 341.69 342.51
1966 342.70 343.65
1967 343.04 345.27
1968 345.11 347.07
1969 347.23 347.62
1970 347.38 349.38
1971 349.29 350.91
1972
La serie presenta tendencia ascendente (no es estacionaria en media)
Es estacionaria en cuanto a la varianza, ya que no se aprecia gran variabilidad