Software para análisis de datos - Ejercicios unidad 2

Ejercicio 1

Elijo el conjunto de datos “iris”.

• 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.

data("iris")
str(iris)

## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

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  
##                 
##                 
## 

data("Orange")
str(Orange)

## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame':   35 obs. of  3 variables:
##  $ Tree         : Ord.factor w/ 5 levels "3"<"1"<"5"<"2"<..: 2 2 2 2 2 2 2 4 4 4 ...
##  $ age          : num  118 484 664 1004 1231 ...
##  $ circumference: num  30 58 87 115 120 142 145 33 69 111 ...
##  - attr(*, "formula")=Class 'formula'  language circumference ~ age | Tree
##   .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv> 
##  - attr(*, "labels")=List of 2
##   ..$ x: chr "Time since December 31, 1968"
##   ..$ y: chr "Trunk circumference"
##  - attr(*, "units")=List of 2
##   ..$ x: chr "(days)"
##   ..$ y: chr "(mm)"

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

library(knitr)

## Warning: package 'knitr' was built under R version 4.4.3

# kable(table(iris$Sepal.Length))
# kable(table(iris$Sepal.Width))
# kable(table(iris$Petal.Length))
# kable(table(iris$Petal.Width))
# kable(table(iris$Species))

# kable(prop.table(iris$Sepal.Length))
# kable(prop.table(iris$Sepal.Width))
# kable(prop.table(iris$Petal.Length))
# kable(prop.table(iris$Petal.Width))
# kable(prop.table(table(iris$Species)))

# kable(table(Orange$Tree))
# kable(table(Orange$age))
# kable(table(Orange$circumference))

# kable(table(Orange$Tree, Orange$age))

Ejercicio 2

1. 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.
2. Usando el vect2 creamos un nuevo vector llamado Fumador que sea un factor con dos niveles. Las etiquetas corresponden a 1 = Fuma y 2 = No fuma.
   
3. Creamos una tabla de contingencia con las dos variables anteriores con el nombre Tabla.
   
4. Miramos la relación de las variables anteriores con la prueba del ji cuadrado.
   
5. Miramos también cómo resulta el test de Fisher.

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)

bajo_peso <- factor(vect1, levels=c(1,2), labels = c("peso_bajo", "peso_normal"))
fumador <- factor(vect2, levels=c(1,2), labels = c("fuma", "no_fuma"))
tabla = table(bajo_peso, fumador)
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

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

1. Activad el paquete de datos airquality del paquete datasets y generad los siguientes
   gráficos:
   

• Un gráfico de dispersión de la variable Ozone de color azul y con almohadillas (#) en vez de puntos.
• Un gráfico de caja de color rojo con la variable Temp con el título Temperatura (en grados Farenheit).

2. Activad el paquete de datos airmiles del paquete datasets y generad los siguientes gráficos:
   

• Un gráfico de líneas de la serie de datos airmiles con el título Datos de pasajeros en vuelos comerciales (en miles) y de color azul (cadetblue2) y la etiqueta del eje x con Miles de pasajeros.

3. Un histograma de la serie de datos airmiles de color marrón (chocolate2). Representad los cuatro gráficos en una única imagen donde los veamos juntos.

data("airquality")
str(airquality)

## 'data.frame':    153 obs. of  6 variables:
##  $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
##  $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
##  $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
##  $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
##  $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...

plot(airquality$Ozone, col="blue", pch="#")

boxplot(airquality$Temp, col="red", main="Temperatura (Farenheit)")

data("airmiles")

str(airmiles)

##  Time-Series [1:24] from 1937 to 1960: 412 480 683 1052 1385 ...

plot(airmiles, col="cadetblue2", main="Datos de pasajeros en vuelos comerciales", xlab="Miles de pasajeros")

hist(airmiles , col="chocolate2")

par(mfrow=matrix(2,2))
plot(airquality$Ozone, col="blue", pch="#")
boxplot(airquality$Temp, col="red", main="Temperatura (Farenheit)")
plot(airmiles, col="cadetblue2", main="Datos de pasajeros en vuelos comerciales", xlab="Miles de pasajeros")
hist(airmiles , col="chocolate2")

## 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().

set.seed(1234)                  
x1 <- rnorm(1000) #simulamos una variable x1    
y1 <- x1 + rnorm(1000) #simulamos una variable y1   
plot(x1, y1, main="Gráfico de dispersión", col="salmon") #creamos el gráfico de dispersión de las variables 
abline(lm(y1~x1), col="red") #creamos la recta de regresión de color verde  

plot(density(x1), main="Densidad de x1", col="aquamarine4") 
lines(density(y1), col="darkred")   

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.
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).
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. c) Con los datos de Birthwt del paquete MASS mostrad la distribución de la variable age que tenga el borde de color negro y azul en el interior de las barras. 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.

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.4.3

ggplot(data = airquality, aes(Wind, Ozone)) +
  geom_point()+
  geom_smooth(method = "lm", color = "darkred")

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot (airquality, aes(Solar.R, Temp, col=factor(Month), shape=factor(Month))) +   
  geom_point() +
  xlab("Solar Radiation")+
  ylab("Temperature")

## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).

birthwt <- MASS::birthwt
ggplot (birthwt, aes(age)) +    
  geom_histogram(fill="blue",col="black")

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

# par(mfrow = c(1,2))
# ggplot (birthwt[birthwt$smoke == 0,], aes(age)) +     
#   geom_histogram(fill="salmon",col="black")
# ggplot (birthwt[birthwt$smoke == 1,], aes(age)) +     
#   geom_histogram(fill="green",col="black")

ggplot (birthwt, aes(age, fill=factor(smoke))) +    
  geom_histogram(binwidth = 5, 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.

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.

# Hemos importado previamente los datos;los analizamos
str(Orange)

## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame':   35 obs. of  3 variables:
##  $ Tree         : Ord.factor w/ 5 levels "3"<"1"<"5"<"2"<..: 2 2 2 2 2 2 2 4 4 4 ...
##  $ age          : num  118 484 664 1004 1231 ...
##  $ circumference: num  30 58 87 115 120 142 145 33 69 111 ...
##  - attr(*, "formula")=Class 'formula'  language circumference ~ age | Tree
##   .. ..- attr(*, ".Environment")=<environment: R_EmptyEnv> 
##  - attr(*, "labels")=List of 2
##   ..$ x: chr "Time since December 31, 1968"
##   ..$ y: chr "Trunk circumference"
##  - attr(*, "units")=List of 2
##   ..$ x: chr "(days)"
##   ..$ y: chr "(mm)"

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

#Hacemos las gráficas de dispersión, considerando que una var. es un factor
pairs(Orange)

ggplot(data=Orange, aes(x=age, y=circumference, colour = Tree))+
  geom_point()

ggplot(data=Orange, aes(y=circumference))+
  geom_boxplot()+
  facet_grid(~Tree)

ggplot(data=Orange, aes(y=circumference))+
  geom_boxplot()+
  facet_grid(~age)

#Estudiamos la correlacion (Pearson coef. de var cuantitativas)
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

#Modelo de regresión simple
model_rl <- lm(data = Orange, formula = circumference ~ age )
summary(model_rl)

## 
## 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

#Gráfica de puntos + recta de regresion
plot(x=Orange$age, y=Orange$circumference, col="lightblue4")
abline(model_rl, col="blue4")

#Predecir con (b0=17.399650 , b1= 0.106770) para age = 600 dias
circ_600 <- predict(model_rl, newdata = data.frame(age=600))
print(circ_600)

##        1 
## 81.46185

Ejercicio 8

Repetimos los gráficos, utilizando una librería distinta a la usada en ej 7.

#Hacemos las gráficas de dispersión, considerando que una var. es un factor
# par(mfrow = c(3,1))
plot(x= Orange$age, y=Orange$circumference, col=Orange$Tree)

boxplot(data = Orange, circumference ~ Tree)

boxplot(data = Orange, circumference ~ age)

#Gráfica de puntos + recta de regresion
ggplot(Orange, aes(x=age, y=circumference))+
  geom_point(col ="lightblue4" )+
  geom_smooth(method="lm", col ="blue4")

## `geom_smooth()` using formula = 'y ~ x'

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).
- ¿Creéis que hay diferencias entre tratamientos?
- Se cumplen las condiciones para poder aplicar una ANOVA. ¿Qué pruebas os planteáis?

#Cargamos los datos
data("PlantGrowth")
str(PlantGrowth)

## 'data.frame':    30 obs. of  2 variables:
##  $ weight: num  4.17 5.58 5.18 6.11 4.5 4.61 5.17 4.53 5.33 5.14 ...
##  $ group : Factor w/ 3 levels "ctrl","trt1",..: 1 1 1 1 1 1 1 1 1 1 ...

levels(PlantGrowth$group)

## [1] "ctrl" "trt1" "trt2"

#Resumen estadistico de los datos
summary(PlantGrowth)

##      weight       group   
##  Min.   :3.590   ctrl:10  
##  1st Qu.:4.550   trt1:10  
##  Median :5.155   trt2:10  
##  Mean   :5.073            
##  3rd Qu.:5.530            
##  Max.   :6.310

# summary(PlantGrowth[PlantGrowth$group == "ctrl","weight"])
# summary(PlantGrowth[PlantGrowth$group == "trt1","weight"])
# summary(PlantGrowth[PlantGrowth$group == "trt2","weight"])
library(psych)

## Warning: package 'psych' was built under R version 4.4.3

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

describeBy(PlantGrowth$weight, PlantGrowth$group, mat = TRUE)

##     item group1 vars  n  mean        sd median trimmed      mad  min  max range
## X11    1   ctrl    1 10 5.032 0.5830914  5.155 5.00500 0.719061 4.17 6.11  1.94
## X12    2   trt1    1 10 4.661 0.7936757  4.550 4.62375 0.533736 3.59 6.03  2.44
## X13    3   trt2    1 10 5.526 0.4425733  5.435 5.50375 0.363237 4.92 6.31  1.39
##          skew  kurtosis        se
## X11 0.2311020 -1.116799 0.1843897
## X12 0.4744580 -1.104735 0.2509823
## X13 0.4847153 -1.160361 0.1399540

boxplot(PlantGrowth$weight ~ PlantGrowth$group , col=palette.colors())

# Valores atípicos
rapportools::rp.outlier(PlantGrowth$weight)

## numeric(0)

#Normalidad de los datos
by(PlantGrowth$weight, PlantGrowth$group, shapiro.test)

## PlantGrowth$group: ctrl
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.95668, p-value = 0.7475
## 
## ------------------------------------------------------------ 
## PlantGrowth$group: trt1
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.93041, p-value = 0.4519
## 
## ------------------------------------------------------------ 
## PlantGrowth$group: trt2
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.94101, p-value = 0.5643

#Homocedasticidad
bartlett.test(PlantGrowth$weight, PlantGrowth$group)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  PlantGrowth$weight and PlantGrowth$group
## Bartlett's K-squared = 2.8786, df = 2, p-value = 0.2371

Si se cumplen los supuestos para aplicar ANOVA ya que p>0.05 en los dos tests.

#Cargamos los datos
anova_mod <- aov(data= PlantGrowth, weight ~group)
# anova_mod <- anova(lm(weight ~ group, data = PlantGrowth))    
summary(anova_mod)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## group        2  3.766  1.8832   4.846 0.0159 *
## Residuals   27 10.492  0.3886                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(anova_mod)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = weight ~ group, data = PlantGrowth)
## 
## $group
##             diff        lwr       upr     p adj
## trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
## trt2-ctrl  0.494 -0.1972161 1.1852161 0.1979960
## trt2-trt1  0.865  0.1737839 1.5562161 0.0120064

Las diferencia significativas se ven entre el trat 1 y el trat 2.

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)

## Warning: package 'plotly' was built under R version 4.4.3

## 
## Attaching package: 'plotly'

## 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

help(package=plotly)
#Grafico interactivo para iris . Seleccionamos la especie y vemos un histrograma del Petal.Length
shiny::selectInput("species", "Especie:", choices = unique(iris$Species), selected = "setosa")

Especie:

setosa versicolor virginica

renderPlotly({
  datos_sel <- subset(iris, Species == input$species)
  plot_ly(data = datos_sel,
          x = ~ Petal.Length,
          type = "histogram",
          nbinsx = 5,
          marker = list(color = "salmon"),
          name = input$species) %>%
    layout(title = paste("Petal.Length para", input$species))
})