Modelamiento y Simulación
1
de Sistemas Ambientales
Semana 6: Programa R
Mg. Yasser Vásquez Baca
Universidad de Huánuco - 2020
2 PROGRAMA R
Es un lenguaje de programación de código abierto
(open source).
RStudio es el Entorno Ingregrado de Desarrollo (IDE). Un
IDE es un software que ayuda a programar de una
manera simple y eficiente, sugiere qué deberíamos
continuar escribiendo, marca dónde nos hemos
confundido y administra nuestros archivos en
proyectos, entre otras funciones.
Requiere una computadora de 4gb de memoria RAM
o superior.
1. Descargar el lenguaje de programación R, la última
versión de R para Windows, Mac o Linux ->
https://cloud.r-project.org/
2. Descargar la interfaz gráfica RStudio Desktop ->
https://rstudio.com/products/rstudio/download/
3 PROGRAMA R
Al abrir R Studio, aparecerán cuatro paneles:
a. Panel de edición de archivos.- permite modificar archivos donde
guardamos nuestro código.
b. Panel de estructuras de datos.- muestra las “estructuras de datos” que
están creadas.
c. Consola o terminal.- Donde se va ejecutado el código.
d. Panel multiuso.- de utilidad para explorar archivos, ver la salida gráfica de
nuestro código, encontrar ayuda sobre nuestros comandos, entre otras
funciones.
4 PROGRAMA R
5 PROGRAMA R
Ejecutar un código en R:
En la consola, ejecutar el siguiente código:
5 * 10 (Enter)
## [1] 50
Otra manera, es escribir exactamente el mismo código, en el Panel de
edición de archivos.
Una vez escrito, seleccionar el código completo con el mouse y presionar
Control + Enter.
En la consola deberían obtener exactamente el mismo resultado.
Al seleccionar el código y apretar Control + Enter, RStudio pasa este código a
la consola y se ejecuta, tal como antes.
6 PROGRAMA R
Guardar un archivo
Presionar Control + S, poner nombre al nuevo archivo, que tendrá la
terminación .R, que es como vamos a identificar a los archivos que
contienen código de R.
Guardar un Nuevo Proyecto (New Project)
RStudio ofrece la posibilidad de guardar todos los archivos vinculados a un
proyecto en una misma dirección de nuestra computadora
Ir a File -> New Project... -> New Directory -> New Project. Allí deben elegir
un nombre para el proyecto y hacer click en Create Project.
7 PROGRAMA R
Obtener la dirección de la carpeta donde está el proyecto:
Desde la consola, ejecutar el código:
getwd()
Devuelve la dirección de la carpeta donde está este proyecto. Lo que
guarden o a lo que intenten acceder desde el proyecto de RStudio se hará
dentro de este directorio.
Importar datos a R
Existen distintas formas de cargar los datos.
8 PROGRAMA R
Archivos separados por comas (CSV - Comma Separated Values)
Uno de los formatos más conocidos para guardar datos. Consiste
simplemente en elementos separados por algún delimitador, en general una
coma, que terminan formando una matriz.
Veamos un ejemplo con el precio de los inmuebles por año y barrio:
precioAvisos <-
read.csv(file="https://raw.githubusercontent.com/martintinch0/CienciaDe
DatosParaCuriosos/master/data/precioBarrios.csv", sep=";", stringsAsFactors
= FALSE)
Asignar resultados de funciones a objetos.
El operador de asignación <- pasa o asigna lo que sea que está a su
derecha a un objeto que está a su izquierda.
En el ejemplo anterior, a la izquierda esta “precioAvisos”, que es como se
llamará el objeto, y a la derecha está nuestra primera función “read.csv()”
9 PROGRAMA R
Organización de datos
El vector es la estructura más básica, es una colección numerada de valores,
y puede accederse a través de índices que denotan el orden de cada valor
dentro del vector. Creemos un vector:
primerVector <- c(20,40,60,80,100)
Colocar la posición entre corchetes [] para observar el valor que representa.
(la numeración arranca desde 1 en R).
primerVector[3]
## [1] 60
10 PROGRAMA R
Los datos pueden ser: números, texto, variables categóricas, variables lógicas.
Todos estos tipos de datos pueden ser representados en vectores:
vectorTexto <- c("Croacia","Argentina","Nigeria","Islandia")
vectorLogico <- c(FALSE, TRUE, TRUE, FALSE)
vectorNumerico <- c(45,90,105)
Al vector numérico también se le dice integer. Podemos hacer múltiples
operaciones aritméticas sobre un vector numérico:
vectorNumerico*2 ; vectorNumerico/5 ; (vectorNumerico*10)-10;
vectorNumerico- vectorNumerico, etc.
11 PROGRAMA R
Podemos conocer medidas de distribución de datos con funciones como
min(), max() y mean(), entre otras. Ej. En el caso de los precios de los
inmuebles ¿Cuál es el valor mínimo y máximo del m2? ¿, ¿Cuál es el precio
promedio?.
Para acceder a una variable específica de un Data Frame, escribir el
operador $
precios2017 <- precioAvisos$USDm2_2017
resumen_2017 <- c(min(precios2017), max(precios2017),
mean(precios2017))
resumen_2017
## [1] 872.7541 6070.3288 2203.1119
12 PROGRAMA R
Otras estructuras de datos importantes en R son las listas y Data Frames:
Las listas contienen un conjunto ordenados de objetos, que pueden ser de
distintas clases.
Los Data Frames son un conjunto de listas, se puede describir como una matriz
en la cual variables (columnas) pueden ser de distintas clases.
Crear la siguiente lista:
lista <- list(Nombres = c("Fernando","Martín"), Apellido="Montané", tienehijos
= FALSE, edad = 26)
lista
13 PROGRAMA R
## $Nombres
## [1] "Fernando" "Martín"
## $Apellido
## [1] "Montané"
## $tienehijos
## [1] FALSE
## $edad
## [1] 26
Se creó una lista que almacena 4 vectores. El primero, es un
vector Nombres que contiene dos elementos de clase character. El segundo,
un vector de un solo elemento, Apellido de clase character. El tercero, un
vector lógico de un elemento (tienehijos). Finalmente, el vector edad, que
posee un elemento númerico.
14 PROGRAMA R
Crear el siguiente data Frame:
df <- data.frame(Nombres = c('Juan','Pedro','Ana','Delfina'),
Edad = c(21,46,58,27),
EstadoCivil = c('Soltero','Casado','Casado','Soltero'),
SecundarioCompleto = c(TRUE, FALSE, TRUE, TRUE))
En las listas y data frames hay tres formas, por lo menos, de acceder a los
objetos:
Usando el signo $ lista1$Nombres
Usando doble corchetes e indicando la posición del objeto buscado
lista1[[1]]
Usando doble corchetes e indicando el nombre del objeto buscado
lista1[['Nombres']]
15 PROGRAMA R
¿Cómo sabemos qué datos tenemos cargados?
Existen al menos dos funciones para inspeccionar rápidamente los datos:
Usar la función View() (Notar la V mayúscula al principio de la función),nos
devuelve la tabla entera en el panel de edición.
View(precioAvisos)
Usar la función str(), nos devuelve las primeras observaciones de cada una de
las variables, también la clase del objeto, la cantidad de observaciones (filas)
y de variables (columnas), así como la clase de cada una de las variables.
str(precioAvisos)
16 PROGRAMA R
Ej. Analizaremos el dataset que contiene el precio promedio del metro
cuadrado en USD para el período 2013-2017 desagregado a nivel de barrios
de la Ciudad de Buenos Aires.
¿Aumentaron los precios de los inmuebles (medido como USD en m2) en el
2017?
Usando el operador $, que nos deja elegir el vector de nuestro data.frame,
solo tenemos que hacer una división.
mean(precioAvisos$USDm2_2017)/mean(precioAvisos$USDm2_2016) - 1
## [1] 0.05299051
Según nuestros datos, aumentaron 5,3% en 2017.
17 PROGRAMA R
¿Todos los barrios aumentaron en la misma proporción o existen
heterogeneidades?
Vamos a agregarle una nueva columna a nuestro data frame
precioAvisos$variacion2017 <-
precioAvisos$USDm2_2017/precioAvisos$USDm2_2016 - 1
View(precioAvisos)
Ordenen desde el visor a las variaciones por barrio y observen las variaciones,
¿Cuál fue el barrio que los precios cayeron un 63%?.
18 PROGRAMA R
Instalando paquetes (packages) en R
Los paquetes son un conjunto de funciones que potencian a R, a partir de la
colaboración de miles de personas. Una vez instalado, queda guardado en la
PC. Para usarlo en cada nueva sesión de R, hay que activar el package con
la función library() o require(), ambas hacen lo mismo.
Se instalan con la función install.packages()
install.packages('tidyverse’)
require(tidyverse)
Tidyverse es un sistema de paquetes para la manipulación, exploración y
visualización de datos que comparten una filosofía de diseño común.
19 PROGRAMA R
Actualizar R
Instalamos el paquete installr
install.packages("installr")
Ejecutamos los siguientes códigos:
library(installr)
updateR()
Primero comprobará si tenemos la última versión. Después diferentes
ventanas de diálogo nos guiarán a lo largo del proceso de instalación,
ofreciéndonos la opción de copiar los paquetes a la nueva ubicación.
20 PROGRAMA R
21
22 PROGRAMA R
Visualizaciones de datos en R
ggplot es un paquete que viene incluido dentro de tidyverse, por lo que cargando
dicho paquete ya van a contar con todas sus funciones dedicadas a crear
decenas de distintos tipos de gráficos.
La nube de puntos, grafico de dispersión o scatter plot es un tipo de gráfico que
por lo general, muestra si hay relación entre dos variables. Ej. Vamos a filtrar los
datos para el año 2007, el último en este dataset de gapminder.
gapminder_df <- read.table(file =
"https://raw.githubusercontent.com/martintinch0/CienciaDeDatosParaCuriosos
/master/data/gapminder.csv",
sep=';',
header = TRUE,
stringsAsFactors = FALSE)
gapminderLastCut <- gapminder_df %>% filter(year==2007)
23 PROGRAMA R
Gráfico de puntos
ggplot empieza escribiendo los datos que
queremos graficar, luego el gráfico que
queremos hacer. Además, necesita saber qué
variable poner en los ejes x e y, estas últimas
indicaciones van dentro de la función aes(), que
es la abreviación de aesthetics.
ggplot(data = gapminderLastCut,
mapping = aes(x=gdpPercap, y = lifeExp)) +
geom_point()
24 PROGRAMA R
Ahora podemos empezar a editar el gráfico. Los títulos de los ejes tienen por
definición el nombre de las variables, vamos a cambiarlo. En ggplot cada
nueva función es agregadas por medio de +.
ggplot(data = gapminderLastCut,
mapping = aes(x=gdpPercap, y = lifeExp)) +
geom_point() +
labs(x = "PIB per cápita", y = "Expectativa de vida al nacer (en años)")
25 PROGRAMA R
Usamos la función labs(), con sus respectivos parámetros x e y. La función labs
también permite agregar título, subtítulos e incluso información sobre la fuente de
nuestro gráfico:
ggplot(data = gapminderLastCut,
mapping = aes(x=gdpPercap, y = lifeExp)) +
geom_point() +
labs(x = "PIB per cápita",
y = "Expectativa de vida al nacer (en años)",
title="A más ingresos mayor tiempo de vida?",
subtitle="Expectativa de vida al nacer según nivel de ingreso",
caption="Fuente: Gapminder")
26 PROGRAMA R
Agregando colores según otras variables
Dentro de la función aes(), también podemos indicar cuál es la variable que
queremos colorear: mediante el argumento color. En el data.frame de gapminder
contamos con la variable: continente.
ggplot(data = gapminderLastCut,
mapping = aes(x=gdpPercap, y = lifeExp,color=continent)) +
geom_point() +
labs(x = "PIB per cápita",
y = "Expectativa de vida al nacer (en años)",
title="A más ingresos mayor tiempo de vida?",
subtitle="Expectativa de vida al nacer según nivel de ingreso",
caption="Fuente: Gapminder")
27
28 PROGRAMA R
Creando paneles
Imaginemos que queremos mostrar el gráfico anterior, pero por cada continente
por separado:
ggplot(data = gapminderLastCut,
mapping = aes(x=gdpPercap, y = lifeExp)) +
geom_point() +
labs(x = "PIB per cápita",
y = "Expectativa de vida al nacer (en años)",
title="A más ingresos mayor tiempo de vida?",
subtitle="Expectativa de vida al nacer según nivel de ingreso",
caption="Fuente: Gapminder") +
facet_wrap(~ continent)
29
30 PROGRAMA R
Gráfico de líneas
En ggplot , analizar la evolución de una determinada variable mediante una línea
de tiempo, se realiza usando geom_line. Ej. del data frame de gapminder,
observar la evolución de la expectativa de vida al nacer promedio por
continente:
promedioContinente <- gapminder_df %>%
group_by(continent,year) %>%
summarise(promedio=mean(lifeExp))
Luego, ya estamos en condiciones de hacer el gráfico:
ggplot(data = promedioContinente,
mapping = aes(x=year, y = promedio,color=continent)) +
geom_line() +
labs(x = "", y = "Expectativa de vida al nacer (en años)",
title="Expectativa de vida al nacer según continente",
caption="Fuente: Gapminder")
31
32 PROGRAMA R
Guardar gráficos de ggplot como objetos:
Sirve para ir agregando “capas” a un gráfico y para exportarlos.
graficoLinea <- ggplot(data = promedioContinente,
mapping = aes(x=year, y = promedio,color=continent)) +
geom_line() + labs(x = "", y = "Expectativa de vida al nacer (en años)",
title="Expectativa de vida al nacer según continente",
caption="Fuente: Gapminder", color =“Continente”)
Exportando gráficos de ggplot
Usar la función ggsave, nos pide un nombre de archivo, y un gráfico a
exportar. Luego, podemos incluir información sobre la calidad del gráfico o
parámetros adicionales.
ggsave(filename = ‘EvolExpContinente.png’,
plot = graficoLinea, dpi = 300)
33 PROGRAMA R
Métodos de análisis de datos
Cargar los datos:
load(url('https://github.com/datalab-
UTDT/GIS2/raw/master/Data/HowellData.RData'))
str(Howell1)
Howell1 es un data.frame que contiene 544 observaciones en cuatro
variables: height (altura), weight (peso), age (edad) y male (genero). Todas
son variables numéricas, incluyendo male, que toma valor 1 cuando la
observación pertenece a un hombre y 0 cuando es a una mujer.
34 PROGRAMA R
Grafiquemos la relación entre la altura y el peso, con ayuda de ggplot2,
paquete de tidyverse para hacer gráficos y también usaremos el
paquete ggthemes, que nos permite hacer gráficos más estéticos.
library(tidyverse)
library(ggthemes)
ggplot(Howell1) +
geom_point(aes(x=weight, y=height)) +
theme_fivethirtyeight() +
labs(x='Peso (kg)', y = 'Altura (m)') +
theme(axis.title = element_text(size=14))
35
Pareciera existir una relación positiva entre altura y peso
36 PROGRAMA R
Regresión lineal
La función lm(), realiza la regresión lineal. Ej. del caso anterior, supongamos
que no existe ninguna relación entre la altura y el peso. Es mas, asumamos
que la altura de una persona solo puede explicarse por un único valor
(constante). Pero antes de hacer eso, saquemos aquellas observaciones
menores a 18 años.
Howell1Adults <- Howell1 %>% filter(age>=18)
regresion1 <- lm(data = Howell1Adults, formula = height ~ 1)
La función lm necesita dos parámetros: los datos y una fórmula. La fórmula es
donde establecemos la relación entre la variable dependiente y las
independientes. A la izquierda de ~ escribimos la variable dependiente, y a su
derecha ponemos las variables independientes.
37 PROGRAMA R
summary(regresion1)
## Call:
## lm(formula = height ~ 1, data = Howell1Adults)
## Residuals:
## Min 1Q Median 3Q Max
## -18.0721 -6.0071 -0.2921 6.0579 24.4729
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 154.5971 0.4127 374.6 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 7.742 on 351 degrees of freedom
38 PROGRAMA R
La función summary nos devuelve un resumen del objeto lm.
Los coeficientes en una relación lineal nos dicen cual es “el cambio
esperado” en la variable dependiente ante una variación de una unidad en
nuestra variable independiente. El coeficiente del ej. tiene un valor de
154.5971.
Ahora modelaremos para observar si existe relación positiva entre altura y
peso.
regresion2 <- lm(data = Howell1Adults, formula = height ~ weight)
summary(regresion2)
Ahora el coeficiente tiene un valor de 0.91, significa que por cada kg
adicional, la altura crece 0.91 cm en promedio.
39 PROGRAMA R
## Call:
## lm(formula = height ~ weight, data = Howell1Adults)
## Residuals:
## Min 1Q Median 3Q Max
## -19.7464 -2.8835 0.0222 3.1424 14.7744
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 113.87939 1.91107 59.59 <2e-16 ***
## weight 0.90503 0.04205 21.52 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 5.086 on 350 degrees of freedom
## Multiple R-squared: 0.5696, Adjusted R-squared: 0.5684
## F-statistic: 463.3 on 1 and 350 DF, p-value: < 2.2e-16
40 PROGRAMA R
El standard error o desviación estándar de un parámetro cuantifica la
incertidumbre, ya que estamos trabajando con una muestra y no con “la
población completa”. La desviación standard comunica si los datos son
representativos del resto de las observaciones.
t. value es el test estadístico que usaremos para testear qué tan probable es
que, dado el error estándar y el valor del coeficiente estimado, este sea en
realidad cero en términos poblacionales.
El Pr(>|t|) testea el valor de nuestro t value, si es muy alto o muy bajo,
entonces la probabilidad de que el parámetro sea cero es muy baja.
No hay interpretación causal en una regresión, como regla general, cada vez
que usen una regresión lineal no van a poder inferir causalidad. Las
regresiones con datos observacionales permiten establecer relaciones entre
las variables.
El Adjusted R-squared (R2), presenta un valor de 0.5684, significa que el
modelo puede explicar o predecir el 57% del movimiento de las alturas.
41 PROGRAMA R
Los intervalos de confianza de la regresión, R los calcula así:
confint(regresion2,level = 0.95)
## 2.5 % 97.5 %
## (Intercept) 110.1207774 117.6380098
## weight 0.8223315 0.9877267
Uno es un límite inferior (2.5%) y otro un límite superior (97.5%) del intervalo de
confianza con un nivel de 95%. Significa que si tomáramos muestras sobre
altura y peso de esta población, en 95 de cada 100 casos, el coeficiente
estaría entre el límite inferior y superior.
42 PROGRAMA R
Grafiquémoslo:
ggplot(Howell1Adults) +
geom_abline(slope = coef(regresion2)[2],intercept = coef(regresion2)[1]) +
geom_point(aes(x=weight, y=height)) +
theme_fivethirtyeight() +
labs(x='Peso (kg)', y = 'Altura (kg)') +
theme(axis.title = element_text(size=14))
43
44 PROGRAMA R
Agreguemos al modelo el sexo de las personas mayores de 18, como variable
independiente
summary(lm(data = Howell1Adults,
formula = height ~ weight + male))
El R2 pasó del 57% al 69% tras la inclusión del sexo ¿Y con el coeficiente de
peso?
45 PROGRAMA R
Predecir un valor particular
En este caso, si tuviésemos que decir cuanto mide una persona que pesa 50kg,
sería multiplicar el coeficiente de peso por 50 y sumarle el intercepto:
prediccion <- coef(regresion2)[1] + coef(regresion2)[2]*50
prediccion
## (Intercept)
## 159.1308
1599.13 cm sería la respuesta ¿Y cuál sería el intervalo de confianza?
Se calcula el intervalo de confianza de la altura puntual con un intervalo del 95%:
predict(regresion2,newdata = data.frame(weight=50), interval ='confidence')
## fit lwr upr
## 1 159.1308 158.4556 159.8061
Si hiciéramos una muestra de 100 personas con peso de 50kg, 95 de ellas caerían
en el intervalo contenido entre 149.1045 y 169.1572.
46
GRACIAS