# Cargamos la librería para el LAB2
library('tidyverse')
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Ejercicio 1: El package de R datasets (https://stat.ethz.ch/R-manual/R-devel/library/datasets/
html/00Index.html) pone a nuestra disposición una serie de conjuntos de
datos con los que poder trabajar como, por ejemplo, iris, cars o
Titanic. Escoged un conjunto de datos. Deberéis: • Buscar un resumen
estadístico de las variables del dataset Iris y Orange. • Generar una
tabla de frecuencias absolutas y una tabla de frecuencias relativas con
el dataset Iris. ¿Todas las tablas generadas tienen sentido para
vosotros? • Generar una tabla de frecuencias absolutas con cada una de
las variables del conjunto de datos Orange. ¿Todas las tablas generadas
tienen sentido para vosotros? • Generar una tabla de doble entrada entre
las variables Tree y Age de Orange.
# Cargar los datasets integrados
data(iris)
data(Orange)
# 1. Resumen estadístico de Iris
summary(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
##
# 2. Resumen estadístico de Orange
summary(Orange)
## Tree age circumference
## 3:7 Min. : 118.0 Min. : 30.0
## 1:7 1st Qu.: 484.0 1st Qu.: 65.5
## 5:7 Median :1004.0 Median :115.0
## 2:7 Mean : 922.1 Mean :115.9
## 4:7 3rd Qu.:1372.0 3rd Qu.:161.5
## Max. :1582.0 Max. :214.0
# Tabla de frecuencia absoluta de la especie
table(iris$Species)
##
## setosa versicolor virginica
## 50 50 50
# Tabla de frecuencia relativa (proporciones)
prop.table(table(iris$Species))
##
## setosa versicolor virginica
## 0.3333333 0.3333333 0.3333333
La tabla de frecuencia absoluta nos dice que hay exactamente 50 flores de cada especie. En cambio, la tabla de frecuencia relativa nos muestra que cada especie representa un \(33.3\%\) (\(0.3333\)) del total. Por lo tanto no tiene sentido hacer todas la tablas porque al ser variables continuas con muchos decimales, la tabla es enorme llena de el mismo decimal. En este caso es mejor hacer una tabla de frecuencia cuando las variables son una categoria o un factor en concreto.
# Frecuencias de la variable Tree (árbol)
table(Orange$Tree)
##
## 3 1 5 2 4
## 7 7 7 7 7
# Frecuencias de la variable age (edad en días)
table(Orange$age)
##
## 118 484 664 1004 1231 1372 1582
## 5 5 5 5 5 5 5
Vemos que hay 5 árboles distintos (del 1 al 5) y que de cada uno se tomaron 7 mediciones a distintas edades
# Tabla cruzada entre el árbol y la edad
table(Orange$Tree, Orange$age)
##
## 118 484 664 1004 1231 1372 1582
## 3 1 1 1 1 1 1 1
## 1 1 1 1 1 1 1 1
## 5 1 1 1 1 1 1 1
## 2 1 1 1 1 1 1 1
## 4 1 1 1 1 1 1 1
Ejercicio 2: Copiad y ejecutad el código siguiente de los dos vectores:
vect1 <- c(1,2,1,2,1,2,1,2,1,2,1,1,1,1,2,2,1,1,2,1)
vect2 <- c(1,1,2,2,2,1,2,1,1,2,1,2,1,1,1,2,1,1,1,1)
Responded a los apartados siguientes: a) Usando el vect1 creamos un nuevo vector llamado Bajo_peso que sea un factor con dos niveles. Las etiquetas corresponden a 1 = Bajo peso y 2 = Peso normal.
# Factor Bajo_Peso
# Convertimos el vector numérico en una variable categórica (factor), indicando que 1 es bajo peso y 2 es peso normal
Bajo_Peso <- factor(vect1, levels = c(1, 2), labels = c("Bajo Peso", "Peso Normal"))
# Hacemos lo mismo con el segundo vector, donde 1 fuma y 2 no fuma
Fumador <- factor(vect2, levels = c(1, 2), labels = c("Fuma", "No Fuma"))
Fumador
## [1] Fuma Fuma No Fuma No Fuma No Fuma Fuma No Fuma Fuma Fuma
## [10] No Fuma Fuma No Fuma Fuma Fuma Fuma No Fuma Fuma Fuma
## [19] Fuma Fuma
## Levels: Fuma No Fuma
# Cruzamos los dos factores para ver cuantos individuos caen en cada combinación
tabla <- table(Bajo_Peso, Fumador)
tabla
## Fumador
## Bajo_Peso Fuma No Fuma
## Bajo Peso 8 4
## Peso Normal 5 3
# Evaluamos si existe asociación/dependencia estadística entre la categoría de peso y el hábito de fumar
chisq.test(tabla)
## Warning in chisq.test(tabla): Chi-squared approximation may be incorrect
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tabla
## X-squared = 0, df = 1, p-value = 1
# Como alternativa más precisa para tablas 2 X 2 con muestras pequeñas o frecuencias reducidas, aplicamos el test de Fisher:
fisher.test(tabla)
##
## Fisher's Exact Test for Count Data
##
## data: tabla
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.1200513 10.9278345
## sample estimates:
## odds ratio
## 1.189031
Ejercicio 3: a) Activad el paquete de datos airquality del paquete datasets y generad los siguientes gráficos:
# Cargar los datasets requeridos
data("airquality")
data("airmiles")
# Ponemos la ventana gráfica en una matriz 2 x 2
par(mfrow = c(2, 2))
• Un gráfico de dispersión de la variable Ozone de color azul y con almohadillas (#) en vez de puntos.
# Gráfico de dispersión de Ozone con símbolos '#' en azul
plot(airquality$Ozone, col = "blue", pch = "#", main = "Dispersión de Ozono")
• Un gráfico de caja de color rojo con la variable Temp con el título Temperatura (en grados Farenheit).
# Diagrama de caja de Temperatura en rojo (Apartado a)
boxplot(airquality$Temp, col = "red", main = "Temperatura (en grados Farenheit)")
# Gráfico de líneas de airmiles
plot(airmiles,
main = "Datos de pasajeros en vuelos comerciales (miles)",
xlab = "Miles de pasajeros",
col = "cadetblue2")
• Un histograma de la serie de datos airmiles de color marrón (chocolate2).
# Histograma de airmiles en color marrón (Apartado b)
hist(airmiles, col = "chocolate2", main = "Histograma de vuelos comerciales")
# Ponemos todos los comandos anteriores juntos:
data("airquality")
data("airmiles")
par(mfrow = c(2, 2))
plot(airquality$Ozone, col = "blue", pch = "#", main = "Dispersión de Ozono")
boxplot(airquality$Temp, col = "red", main = "Temperatura (en grados Farenheit)")
plot(airmiles, main = "Datos de pasajeros en vuelos comerciales (miles)", xlab = "Miles de pasajeros", col = "cadetblue2")
hist(airmiles, col = "chocolate2", main = "Histograma de vuelos comerciales")
par(mfrow = c(1, 1))
Ejercicio 4: Intentad reproducir el ejemplo 7 de este LAB con
unos datos simulados por vosotros. No
hace falta que sean parecidos, pero es necesario que podáis hacer
diferentes gráficos
estadísticos con el comando plot()
# Fijar semilla para reproducibilidad
set.seed(333)
# Primero simulamos dos variables normales correlacionadas
x1 <- rnorm(500)
y1 <- x1 + rnorm(500)
# Después dibujamos el gráfico de dispersión con recta de regresión
plot(x1, y1, main = "Gráfico de dispersión de datos simulados")
abline(lm(y1 ~ x1), col = "orange")
# Por último el gráfico de densidad de x1 con la curva de y1 superpuesta
plot(density(x1), main = "Densidad de x1 e y1", col = "red")
lines(density(y1), col = "blue")
Ejercicio 5: Para poder practicar la creación de gráficos con ggplot2, vamos a crear seis gráficos con diferentes características y con diferentes conjuntos de datos de paquetes trabajados anteriormente.
library(ggplot2)
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
# Cargar datasets
data("airquality")
data("birthwt")
a) Con los datos airquality del paquete datasets cread un gráfico de dispersión simple, pero con la recta de regresión ajustada en color rojo (darkred).
# Dispersión con recta de regresión ajustada en color rojo oscuro
ggplot(airquality, aes(Solar.R, Temp)) +
geom_point(size = 3) +
geom_smooth(method = "lm", col = "darkred")
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 7 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).
b) Con los datos airquality del paquete datasets cread un gráfico de dispersión con las variables Solar.R y Temp con sus respectivas etiquetas de eje. Los colores y la forma de los puntos deben ser en función de la variable Month convertida a factor.
# Dispersión variando color y forma según el mes (Month como factor)
ggplot(airquality, aes(x = Solar.R, y = Temp, col = factor(Month), shape = factor(Month))) +
geom_point()
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).
c) Con los datos de Birthwt del paquete MASS mostrad la distribución de la variable ageque tenga el borde de color negro y azul en el interior de las barras.
# Histograma de edad materna con interior azul y bordes negros
ggplot(birthwt, aes(age)) +
geom_histogram(binwidth = 5, fill = "blue", col = "black")
d) Con los datos de Birthwt del paquete MASS mostrad la distribución de la variable age,pero separada en dos gráficos (fumador/no fumador) por la variable smoke.
# Histograma de edad separado en dos paneles según el hábito de fumar (smoke)
ggplot(birthwt, aes(age)) +
geom_histogram(binwidth = 5, fill = "violet", col = "black") +
facet_grid(~smoke)
Ejercicio 6: Cread un gráfico con unos datos extraídos del
paquete Datasets de R y guardadlo como imagen (.jpg) con el nombre
migrafic1 y, también, como documento en PDF con el nombre migrafic2.
Haced una captura de pantalla del fichero generado.
# Primero guardamos como imagen JPG
jpeg("migrafic1.jpg")
data(iris)
boxplot(iris$Petal.Width, col="green")
dev.off()
## quartz_off_screen
## 2
# Segundo guardar como documento PDF
data(iris)
boxplot(iris$Petal.Width, col="green")
dev.off()
## null device
## 1
Ejercicio 7: Para practicar la regresión lineal simple usaremos el conjunto de datos Orange que se encuentra en la librería tidyverse y que tiene información sobre tres variables (árbol, edad en días desde que se sembró el árbol y circunferencia del tronco en centímetros) de 35 naranjos. • Queremos saber qué valor de circunferencia tendrá un árbol seiscientos días después de plantarlo. Cread para ello un modelo lineal y practicad la regresión lineal paso a paso con los pasos que habéis visto en este laboratorio.
library(tidyverse)
# Primero cargamos y exploramos el conjunto de datos
data("Orange")
head(Orange)
## Grouped Data: circumference ~ age | Tree
## Tree age circumference
## 1 1 118 30
## 2 1 484 58
## 3 1 664 87
## 4 1 1004 115
## 5 1 1231 120
## 6 1 1372 142
summary(Orange)
## Tree age circumference
## 3:7 Min. : 118.0 Min. : 30.0
## 1:7 1st Qu.: 484.0 1st Qu.: 65.5
## 5:7 Median :1004.0 Median :115.0
## 2:7 Mean : 922.1 Mean :115.9
## 4:7 3rd Qu.:1372.0 3rd Qu.:161.5
## Max. :1582.0 Max. :214.0
# Después el gráfico de dispersión múltiple para ver relaciones
pairs(Orange)
# Ejecutamos el test de correlación de Pearson entre edad y circunferencia
cor.test(Orange$age, Orange$circumference)
##
## Pearson's product-moment correlation
##
## data: Orange$age and Orange$circumference
## t = 12.9, df = 33, p-value = 1.931e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8342364 0.9557955
## sample estimates:
## cor
## 0.9135189
# Ajustamos el modelo de regresión lineal simple
mod_reg <- lm(circumference ~ age, data = Orange)
summary(mod_reg)
##
## Call:
## lm(formula = circumference ~ age, data = Orange)
##
## Residuals:
## Min 1Q Median 3Q Max
## -46.310 -14.946 -0.076 19.697 45.111
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.399650 8.622660 2.018 0.0518 .
## age 0.106770 0.008277 12.900 1.93e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 23.74 on 33 degrees of freedom
## Multiple R-squared: 0.8345, Adjusted R-squared: 0.8295
## F-statistic: 166.4 on 1 and 33 DF, p-value: 1.931e-14
# Por último, el gráfico de dispersión con la recta de regresión ajustada
plot(Orange$age, Orange$circumference,
xlab = "Edad (días)",
ylab = "Circunferencia (cm)",
main = "Regresión Lineal: Edad vs Circunferencia")
abline(mod_reg, col = "blue")
# Y la predicción para un árbol de 600 días
dias <- 600
medida <- 0.1068 * dias + 17.3997
print(medida)
## [1] 81.4797
Ejercicio 8: Repetid los gráficos (aquellos que podáis) de la regresión lineal simple del ejercicio anterior pero ahora con ggplot2. Nota: en caso de haber elaborado el ejercicio 7 con gráficos de este tipo, creadlos con gráficos base.
library(tidyverse)
# Cargar datos
data("Orange")
# Representamos el gráfico de dispersión con la recta de regresión y una marca vertical
Orange %>%
ggplot(aes(x = age, y = circumference)) +
geom_point() +
geom_abline(intercept = 17.3997, slope = 0.1068, col = 'blue') +
geom_vline(xintercept = 800, col = 'red')
Ejercicio 9: Activad el conjunto de datos PlantGrowth del
paquete datasets de R. Este archivo tiene los recultados de un
experimento para comparar los rendimientos medios por el peso seco de
las plantas (weight) obtenidos bajo un control y dos condiciones de
tratamiento diferentes (group factor).
# Cargar datos
data("PlantGrowth")
# Realizar el ANOVA para comparar el peso seco entre grupos de tratamiento
anova(lm(weight ~ group, data = PlantGrowth))
## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## group 2 3.7663 1.8832 4.8461 0.01591 *
## Residuals 27 10.4921 0.3886
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
• ¿Creéis que hay diferencias entre tratamientos? Sí, existen diferencias estadísticamente significativas en el peso seco medio de las plantas entre los distintos grupos, por lo que se rechaza la hipótesis nula de igualdad de medias
• Se cumplen las condiciones para poder aplicar una ANOVA. ¿Qué pruebas os planteáis? Necesitamos que se cumplan 2 normas para aplicar ANOVA. Los datos tienen que tener forma de campana normal y que los grupos tengan una dispersión parecida Para saber qué grupo gana a otros grupos podemos hacer la prueba de Tukey
Ejercicio 10: Buscad información del paquete Plotly para la creación de gráficos interactivos y generad un histograma o un gráfico de barras interactivo. Explicad qué se puede hacer con este gráfico.
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:MASS':
##
## select
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
# Gráfico de barras interactivo con los datos de iris
plot_ly(data = iris, x = ~Species, y = ~Petal.Length, type = "bar")
En este gráfico podemos: - Ver información al pasar el ratón por las barras y su valor exacto - Permite hacer zoom, desplazar el gráfico e incluso aislar algún dato
Caso práctico A partir de unos datos bioclínicos o biosanitarios que escojáis y que importéis a R, explicad sus variables (mínimo de ocho variables) y también: • Realizad un resumen estadístico completo del dataset y explicad los resultados. • Realizad cinco gráficos básicos con las variables, explicad su significado y guardadlos como imágenes (jpeg o bmp). • Realizad dos gráficos con el comando ggplot(), explicad su significado y guardadlos como imágenes (jpeg o bmp). • Generad una regresión lineal entre dos de sus variables paso a paso y comentad los resultados obtenidos. Nota: os dejamos diferentes enlaces con repositorios de datos biosanitarios en abierto: • Repositorios para aprendizaje y docencia en bioestadística • Datos abiertos biosanitarios. Biblioteca de Berkeley
library(ggplot2)
# Primero cargamos el dataset elegido de EE.UU. (50 estados y 8 variables)
datos_usa <- data.frame(state.x77)
# Hacemos un resumen estadístico
summary(datos_usa)
## Population Income Illiteracy Life.Exp
## Min. : 365 Min. :3098 Min. :0.500 Min. :67.96
## 1st Qu.: 1080 1st Qu.:3993 1st Qu.:0.625 1st Qu.:70.12
## Median : 2838 Median :4519 Median :0.950 Median :70.67
## Mean : 4246 Mean :4436 Mean :1.170 Mean :70.88
## 3rd Qu.: 4968 3rd Qu.:4814 3rd Qu.:1.575 3rd Qu.:71.89
## Max. :21198 Max. :6315 Max. :2.800 Max. :73.60
## Murder HS.Grad Frost Area
## Min. : 1.400 Min. :37.80 Min. : 0.00 Min. : 1049
## 1st Qu.: 4.350 1st Qu.:48.05 1st Qu.: 66.25 1st Qu.: 36985
## Median : 6.850 Median :53.25 Median :114.50 Median : 54277
## Mean : 7.378 Mean :53.11 Mean :104.46 Mean : 70736
## 3rd Qu.:10.675 3rd Qu.:59.15 3rd Qu.:139.75 3rd Qu.: 81162
## Max. :15.100 Max. :67.30 Max. :188.00 Max. :566432
# Representamos cinco gráficos básicos (R base)
hist(datos_usa$Income, col = "blue", main = "1. Histograma de Ingresos", xlab = "Ingresos ($)")
boxplot(datos_usa$Life.Exp, col = "green", main = "2. Esperanza de Vida")
plot(datos_usa$Illiteracy, datos_usa$Murder, main = "3. Analfabetismo vs Homicidios", xlab = "% Analfabetismo", ylab = "Tasa Homicidios")
barplot(head(datos_usa$Population, 10), names.arg = rownames(head(datos_usa, 10)), las = 2, col = "orange", main = "4. Población (Top 10 estados)")
plot(density(datos_usa$HS.Grad), col = "pink", main = "5. Densidad de Graduados Secundaria (%)")
# Hacemos dos gráficos con ggplot2
ggplot(datos_usa, aes(x = Illiteracy, y = Murder)) +
geom_point(color = "red", size = 2.5) +
labs(title = "Relación entre Analfabetismo y Homicidios", x = "% Analfabetismo", y = "Homicidios")
ggplot(datos_usa, aes(x = Income, y = Life.Exp)) +
geom_point(color = "brown") +
geom_smooth(method = "lm", color = "yellow") +
labs(title = "Esperanza de Vida según Ingresos", x = "Ingresos ($)", y = "Esperanza de Vida (años)")
## `geom_smooth()` using formula = 'y ~ x'
# Hacemos una regresión lineal simple: tasa de Homicidios según analfabetismo
mod_usa <- lm(Murder ~ Illiteracy, data = datos_usa)
summary(mod_usa)
##
## Call:
## lm(formula = Murder ~ Illiteracy, data = datos_usa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5315 -2.0602 -0.2503 1.6916 6.9745
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.3968 0.8184 2.928 0.0052 **
## Illiteracy 4.2575 0.6217 6.848 1.26e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.653 on 48 degrees of freedom
## Multiple R-squared: 0.4942, Adjusted R-squared: 0.4836
## F-statistic: 46.89 on 1 and 48 DF, p-value: 1.258e-08
plot(datos_usa$Illiteracy, datos_usa$Murder, main = "Regresión Lineal: Analfabetismo vs Homicidios", xlab = "% Analfabetismo", ylab = "Homicidios")
abline(mod_usa, col = "red", lwd = 2)