Gráficos de barras

Este tutorial de R describe cómo crear un gráfico de barras utilizando el software R y el paquete ggplot2. Se puede utilizar la función geom_bar().

df <- data.frame(Empresa = 
                   c("Disney World", 
                      "Florida Hospital", 
                      "Publix Supermarkets Inc",
                      "Walmart Stores Ind",
                      "Univaersal Orlando"), 
                 Asalariados = 
                   c(51600,
                     19283,
                     14995,
                     14995,
                     12000))
knitr::kable(df)

View(df)
head(df)

##                   Empresa Asalariados
## 1            Disney World       51600
## 2        Florida Hospital       19283
## 3 Publix Supermarkets Inc       14995
## 4      Walmart Stores Ind       14995
## 5      Univaersal Orlando       12000

library(ggplot2)
ggplot(data=df, aes(x=Empresa, y=Asalariados)) + geom_bar(stat="identity")

ggplot(data=df, aes(x=Empresa, y=Asalariados)) + geom_bar(stat="identity") + coord_flip()

ggplot(data=df, aes(x=Empresa, y=Asalariados)) + geom_bar(stat="identity", width=0.5)

ggplot(data=df, aes(x=Empresa, y=Asalariados)) + geom_bar(stat="identity", color="blue", fill="white")

p <- ggplot(data=df, aes(x=Empresa, y=Asalariados)) + geom_bar(stat="identity", fill="steelblue") + theme_minimal()
p

library(ggplot2)
library(dplyr)

specie <- c(rep("sorgho" , 3), 
            rep("poacee" , 3), 
            rep("banana" , 3), 
            rep("triticum" , 3))
condition <- rep(c("normal" , "stress" , "Nitrogen") , 4)
value <- abs(rnorm(12 , 0 , 15))
data <- data.frame(specie, condition, value)

knitr::kable(data)

head(data)

##   specie condition    value
## 1 sorgho    normal 12.30164
## 2 sorgho    stress 12.41942
## 3 sorgho  Nitrogen 13.84550
## 4 poacee    normal 21.35051
## 5 poacee    stress 17.25224
## 6 poacee  Nitrogen 11.30622

ggplot(data, aes(fill=condition, y=value, x=specie)) +
  geom_bar(position="dodge", stat="identity")

ggplot(data, aes(fill=condition, y=value, x=specie)) + 
  geom_bar(position="stack", stat="identity")

df <- data.frame(Especialidad = c(rep("Finanzas", 2), 
                                  rep("Marketing", 2), 
                                  rep("Contabilidad", 2)),
                 Año = rep(c(2000, 2005), 3),
                 Matriculados = c(160, 250,
                                  140, 200,
                                  100, 150))
knitr::kable(df)

head(df, n = -4L)

##   Especialidad  Año Matriculados
## 1     Finanzas 2000          160
## 2     Finanzas 2005          250

head(df, n = 4L)

##   Especialidad  Año Matriculados
## 1     Finanzas 2000          160
## 2     Finanzas 2005          250
## 3    Marketing 2000          140
## 4    Marketing 2005          200

ggplot(df, aes(fill = Especialidad, y = Matriculados, x = Año)) +
  geom_bar(position="dodge", stat="identity", color = "black") +
  scale_fill_manual(values = c("red", "yellow","green"))

ggplot(df, aes(fill = Especialidad, y = Matriculados, x = Año)) +
  geom_bar(position="stack", stat="identity", color = "black") +
  scale_fill_manual(values = c("red", "yellow","green"))

Gráfico circular

• ggplot2 no ofrece ningún geom específica para construir diagramas circulares (piecharts). El truco es el siguiente: El marco de datos de entrada tiene 2 columnas: los nombres de los grupos (group here) y su valor (value here), se construye un gráfico de barras apilado con una sola barra utilizando la función geom_bar(), luego se hace circular con coord_polar()

• Dataset Ejemplo 14: El gerente de una universidad pidió una desagregación de los gastos de viaje de los profesores que asistían a diversas reuniones profesionales. Se observó que el 31 por ciento de los gastos estaba representado por los costes de transporte, el 25 por ciento por los costes de alojamiento, el 12 por ciento por los gastos de alimentación, el 20 por ciento por los gastos de matrícula y el resto por costes varios. Represente gráficamente estos datos.

library(ggplot2)
library(dplyr)

data <- data.frame(group = 
                     c("Transporte", 
                       "Alojamiento", 
                       "Alimentaión", 
                       "Gastos de matricula",
                       "Varios"), 
                   value = c(31, 25, 12, 20, 12))
knitr::kable(data)

ggplot(data, aes(x = "", y = value, fill=group)) +
  geom_bar(stat = "identity", width = 1) +
  coord_polar("y", start = 0)

• Renombrar etiquetas: La parte complicada es calcular la posición \(y\) de las etiquetas utilizando la transformación coord_polar. Para esto procedemos de la siguiente manera

library(ggplot2)
library(dplyr)

data <- data %>% 
  arrange(desc(group)) %>%
  mutate(prop = value / sum(data$value) *100) %>%
  mutate(ypos = cumsum(prop)- 0.5*prop )

require(scales)
ggplot(data, aes(x="", y = prop, fill=group)) +
  geom_bar(stat="identity", width=1, color="white") +
  coord_polar("y", start=0) +
  theme_void() + 
  theme(legend.position="none") +
  
  geom_text(aes(y = ypos, label = percent(value/100)), color = "white", size=6) +
  scale_fill_brewer(palette="Set1")

Diagrama de pareto

• El diagrama de Pareto es una combinación de un diagrama de barras y un diagrama de líneas que se utiliza para la visualización. En los gráficos de Pareto, el eje vertical derecho se utiliza para la frecuencia acumulada, mientras que el eje vertical izquierdo representa la frecuencia.

• Primero creamos un dataset y cargamos las librerias necesarias. Ejemplo 14: Identificar las principales causas de los problemas e intentar corregirlas rápidamente con un coste mínimo a menudo.

df <- data.frame(Error =
                   c("Códigos de procedimientos y diagnósticos",
                     "Información de proveedor", 
                     "Información del paciente",
                     "Tablas de precios",
                     "Solicitudes de contratos",
                     "Ajustes de los proveedores",
                     "Otros"), 
                 Frecuencia = c(40, 9, 6, 17, 37, 7, 4))
knitr::kable(df)

head(df)

##                                      Error Frecuencia
## 1 Códigos de procedimientos y diagnósticos         40
## 2                 Información de proveedor          9
## 3                 Información del paciente          6
## 4                        Tablas de precios         17
## 5                 Solicitudes de contratos         37
## 6               Ajustes de los proveedores          7

• El siguiente llamado nos entregará el gráfico de pareto para el dataframe df

library(qcc)

Frecuencia <- df$Frecuencia
names(Frecuencia) <- df$Error 

pareto.chart(Frecuencia, 
             ylab="Frecuencia",
             col = heat.colors(length(Frecuencia)),
             cumperc = seq(0, 100, by = 10),
             ylab2 = "Porcentaje acumulado",
             main = "Grafico de Pareto para Errores"
)

##                                           
## Pareto chart analysis for Frecuencia
##                                             Frequency  Cum.Freq. Percentage
##   Códigos de procedimientos y diagnósticos  40.000000  40.000000  33.333333
##   Solicitudes de contratos                  37.000000  77.000000  30.833333
##   Tablas de precios                         17.000000  94.000000  14.166667
##   Información de proveedor                   9.000000 103.000000   7.500000
##   Ajustes de los proveedores                 7.000000 110.000000   5.833333
##   Información del paciente                   6.000000 116.000000   5.000000
##   Otros                                      4.000000 120.000000   3.333333
##                                           
## Pareto chart analysis for Frecuencia
##                                            Cum.Percent.
##   Códigos de procedimientos y diagnósticos    33.333333
##   Solicitudes de contratos                    64.166667
##   Tablas de precios                           78.333333
##   Información de proveedor                    85.833333
##   Ajustes de los proveedores                  91.666667
##   Información del paciente                    96.666667
##   Otros                                      100.000000

Diagrama de tallo y hojas

• La sintaxis básica para dibujar el gráfico de tallos y hojas en R es la siguiente:

stem(x, scale = 1, width = 80, atom = 1e-08)

• La siguiente es la lista de argumentos soportados por el gráfico de tallo y hoja en el lenguaje de programación R:
  • x: Especifique los datos sobre los que desea dibujar el gráfico tallo y hojas. En este caso, debe utilizar el vector numérico, o una lista que contenga el vector numérico.
  • scale: Especifique la escala que desea utilizar para su gráfico.
  • width: Es opcional, pero puede usarla para especificar la anchura deseada de la parcela. Por defecto, es 80.
  • atom: Es la tolerancia.
• Consideremos el Ejemplo 17: de los apuntes de clase relacionado con la colección de 25 calificaciones en un examen de álgebra. Nótese que el resultado es similar al de los apuntes, las única diferencia es que stem() organiza las hojas de menor a mayor

df <- data.frame(calificaciones = c(78, 67, 65, 87, 75, 65, 71, 54, 94, 64, 84, 82, 81, 
                                    68, 85, 76, 89, 98, 59, 57, 79, 65, 59, 80, 67))

head(df)

##   calificaciones
## 1             78
## 2             67
## 3             65
## 4             87
## 5             75
## 6             65

knitr::kable(head(df))

• Creación del gráfico de tallo y hoja en R

stem(df$calificaciones)

## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##   5 | 4799
##   6 | 4555778
##   7 | 15689
##   8 | 0124579
##   9 | 48

• Al observar el d. de t. h. anterior podemos concluir que:
  • La calificación más alta es 98.
  • La menor es 54.
  • Las calificaciones varían de 54 a 98.
  • El tallo 9 tiene menos hojas.
  • Los tallos 6 y 8 contienen más hojas, siete en cada uno.
  • El número total de hojas representa el tamaño de la muestra.
• En el siguiente ejemplo, mostramos cómo hacer un gráfico de tallo y hoja en R utilizando el conjunto de datos ChickWeight, que proporciona R Studio. Mantenemos por defecto cada uno de los argumentos

head(ChickWeight)

##   weight Time Chick Diet
## 1     42    0     1    1
## 2     51    2     1    1
## 3     59    4     1    1
## 4     64    6     1    1
## 5     76    8     1    1
## 6     93   10     1    1

knitr::kable(head(ChickWeight))

• Creación del gráfico de tallo y hoja en R

stem(ChickWeight$weight)

## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##    2 | 599999999
##    4 | 00000111111111111111111112222222222222223333456678888888899999999999+38
##    6 | 00111111122222222333334444455555666677777888888900111111222222333334+8
##    8 | 00112223344444455555566777788999990001223333566666788888889
##   10 | 0000111122233333334566667778889901122223445555667789
##   12 | 00002223333344445555667788890113444555566788889
##   14 | 11123444455556666677788890011234444555666777777789
##   16 | 00002233334444466788990000134445555789
##   18 | 12244444555677782225677778889999
##   20 | 0123444555557900245578
##   22 | 0012357701123344556788
##   24 | 08001699
##   26 | 12344569259
##   28 | 01780145
##   30 | 355798
##   32 | 12712
##   34 | 1
##   36 | 13

• Se puede observar que la salida no es igual a los datos del ejemplo. Esto se debe a que los tallos están agrupados (el primer tallo es para 2 y 3, el segundo para 4 y 5, etc.). Para resolver este problema, se debe cambiar la altura del diagrama con el argumento scale de la siguiente manera:

stem(ChickWeight$weight, scale = 2)

## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##    3 | 599999999
##    4 | 00000111111111111111111112222222222222223333456678888888899999999999
##    5 | 000000111111112222333334445555566667778888899999
##    6 | 001111111222222223333344444555556666777778888889
##    7 | 0011111122222233333444444446667778889999
##    8 | 0011222334444445555556677778899999
##    9 | 0001223333566666788888889
##   10 | 00001111222333333345666677788899
##   11 | 01122223445555667789
##   12 | 0000222333334444555566778889
##   13 | 0113444555566788889
##   14 | 1112344445555666667778889
##   15 | 0011234444555666777777789
##   16 | 0000223333444446678899
##   17 | 0000134445555789
##   18 | 1224444455567778
##   19 | 2225677778889999
##   20 | 01234445555579
##   21 | 00245578
##   22 | 00123577
##   23 | 01123344556788
##   24 | 08
##   25 | 001699
##   26 | 12344569
##   27 | 259
##   28 | 0178
##   29 | 0145
##   30 | 35579
##   31 | 8
##   32 | 127
##   33 | 12
##   34 | 1
##   35 | 
##   36 | 1
##   37 | 3

Gráfico de series temporales

• Considere la siguiente serie de tiempo correspondiente a la evolución de la temperatura global en superficie con respecto a las temperaturas medias de 1951-1980, fuente NASA/GISS.

• La función geom_line() de ggplot2 se encarga de conectar las observaciones en orden de la variable en el eje \(x\). La siguientes instrucciones son usada para representar la gráfica de serie de tiempo, primero cargamos los datos usando read_csv()

library(readr)
library(knitr)  

df <- read_csv("https://raw.githubusercontent.com/lihkir/AnalisisEstadisticoUN/main/Data/annual_csv.csv")
knitr::kable(head(df))

head(df, n = 10L)

## # A tibble: 10 x 3
##    Source   Year  Mean
##    <chr>   <dbl> <dbl>
##  1 GCAG     2016 0.936
##  2 GISTEMP  2016 0.99 
##  3 GCAG     2015 0.900
##  4 GISTEMP  2015 0.87 
##  5 GCAG     2014 0.741
##  6 GISTEMP  2014 0.74 
##  7 GCAG     2013 0.668
##  8 GISTEMP  2013 0.65 
##  9 GCAG     2012 0.624
## 10 GISTEMP  2012 0.63

library(ggplot2)
library(dplyr)

ggplot(df, aes(x = Year, y = Mean)) +
  geom_line(color="#69b3a2", size = 1) +
  xlab("Year")

• Cargue de datos para Informe No 1

library(readr)
library(knitr)  

df <- read.csv("https://raw.githubusercontent.com/lihkir/AnalisisEstadisticoUN/main/Data/catalog.csv")

knitr::kable(head(df))

Tablas de frecuencia

data <- c(2, 3, 0, 1, 5,
          3, 2, 3, 0, 0,
          2, 1, 2, 1, 0,
          2, 1, 1, 1, 3,
          4, 0, 0, 2, 1)
knitr::kable(head(data))

library(questionr)

table <- questionr::freq(data, cum = TRUE, sort = "dec", total = TRUE)
knitr::kable(table)

str(table) 

## Classes 'freqtab' and 'data.frame':  7 obs. of  5 variables:
##  $ n      : num  7 6 6 4 1 1 25
##  $ %      : num  28 24 24 16 4 4 100
##  $ val%   : num  28 24 24 16 4 4 100
##  $ %cum   : num  28 52 76 92 96 100 100
##  $ val%cum: num  28 52 76 92 96 100 100

x <- row.names(table)
y <- table$n
names <- x[1:(length(x)-1)]
freqs <- y[1:(length(y)-1)]

df <- data.frame(x = names, y = freqs)
knitr::kable(df)

library(ggplot2)

ggplot(data=df, aes(x=x, y=y)) + 
  geom_bar(stat="identity", color="blue", fill="yellow") +
  xlab("Número de asistencias") +
  ylab("Frecuencia")

edades <- c(22, 19, 16, 13, 18, 15, 20, 14, 15, 16,
          15, 16, 20, 13, 15, 18, 15, 13, 18, 15)
knitr::kable(head(edades))

n_sturges = 1 + log(length(edades))/log(2)
n_sturgesc = ceiling(n_sturges)
n_sturgesf = floor(n_sturges)

n_clases = 0
if (n_sturgesc%%2 == 0) {
  n_clases = n_sturgesf
} else {
  n_clases = n_sturgesc
}
R = max(edades) - min(edades)
w = ceiling(R/n_clases)

bins <- seq(min(edades), max(edades) + w, by = w)
bins

## [1] 13 15 17 19 21 23

Edades <- cut(edades, bins)
Freq_table <- transform(table(Edades), Rel_Freq=prop.table(Freq), Cum_Freq=cumsum(Freq))
knitr::kable(Freq_table)

str(Freq_table)

## 'data.frame':    5 obs. of  4 variables:
##  $ Edades  : Factor w/ 5 levels "(13,15]","(15,17]",..: 1 2 3 4 5
##  $ Freq    : int  7 3 4 2 1
##  $ Rel_Freq: num  0.4118 0.1765 0.2353 0.1176 0.0588
##  $ Cum_Freq: int  7 10 14 16 17

df <- data.frame(x = Freq_table$Edades, y = Freq_table$Freq)
knitr::kable(df)

library(ggplot2)

ggplot(data=df, aes(x=x, y=y)) +
  geom_bar(stat="identity", color="blue", fill="green") +
  xlab("Rango de Edades") +
  ylab("Frecuencia")

Estadísticos

• Una función multiuso muy útil en R es summary(X), donde X puede ser uno de cualquier número de objetos, incluyendo conjuntos de datos, variables y modelos lineales, por nombrar algunos. Cuando se utiliza, el comando proporciona datos de resumen relacionados con el objeto individual que se introdujo en él. Así, la función de resumen tiene diferentes resultados dependiendo del tipo de objeto que tome como argumento. Además de ser ampliamente aplicable, este método es valioso porque a menudo proporciona exactamente lo que se necesita en términos de estadísticas de resumen.

• Considere el Ejemplo 65 de las notas de clase, donde se entrega la siguiente tabla con información del precio del cerdo y queso en las capitales más importantes del mundo

df <- data.frame(Capital = c("Berna", "Bonn", "Brasilia", "BuenosAires",
                             "Camberra", "Londres", "Madrid", "México", "Ottawa",
                             "Paris", "Pretoria", "Roma", "Estocolmo", "Tokio",
                             "Washington"),
                 Cerdo = c(6.61, 2.38, 1.27, 1.36, 2.06, 1.56, 2.33, 1.08, 1.99,
                           2.47, 1.95, 2.46, 5.35, 4.19, 3.29),
                 Queso = c(4.00, 2.74, 1.08, 2.03, 2.60, 1.81, 3.15, 2.29, 3.98,
                           2.37, 1.76, 2.96, 2.54, 2.38, 2.69))
knitr::kable(df)

• Usando la función summary() podemos obtener estadísticos de interes y valores de posición:

summary(df$Cerdo)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.080   1.755   2.330   2.690   2.880   6.610

• Como se puede notar la función summary() no nos entrega todos los estadísticos de interés, para solucionar esto podemos hacer uso de la librería, pastecs y la función stat.desc(), como se muestra a continuación

library(pastecs)
stat.desc(df)

##          Capital      Cerdo      Queso
## nbr.val       NA 15.0000000 15.0000000
## nbr.null      NA  0.0000000  0.0000000
## nbr.na        NA  0.0000000  0.0000000
## min           NA  1.0800000  1.0800000
## max           NA  6.6100000  4.0000000
## range         NA  5.5300000  2.9200000
## sum           NA 40.3500000 38.3800000
## median        NA  2.3300000  2.5400000
## mean          NA  2.6900000  2.5586667
## SE.mean       NA  0.4051325  0.2007673
## CI.mean       NA  0.8689229  0.4306029
## var           NA  2.4619857  0.6046124
## std.dev       NA  1.5690716  0.7775682
## coef.var      NA  0.5832980  0.3038959

Caja y extensión

• Los gráficos de caja (box plots), también conocidos como diagramas de cajas y bigotes, son una representación gráfica que permite resumir las características principales de los datos (posición, dispersión, asimetría, …) e identificar la presencia de valores atípicos. En esta sección revisaremos cómo hacer box plots en R base y en ggplot2.

• Considere el Ejemplo 73 de los apuntes de clase, donde se entrega un conjunto de datos distribuidos en una tabla de frecuencias. Utilizaremos las funciones boxplot() y geom_boxplot() para construir gráficos de caja. Cualquiera de las dos opciones son validas y se pueden elegir de acuerdo a la preferencia

• Utilizando boxplot() R base

data <- c(47, 52, 52, 57, 58, 58, 60, 65, 66, 66, 71, 71, 72, 73, 96)
boxplot(data, horizontal=TRUE, col='steelblue')

• Usando geom_boxplot() de la librería ggplot2

library(tidyverse)
library(hrbrthemes)
library(viridis)

df <- data.frame(data)
df %>% ggplot(aes(x = "", y = data)) +
  geom_boxplot(color="red", fill="orange", alpha=0.5) +
  theme_ipsum() +
  theme(legend.position="none", plot.title = element_text(size=11)) +
  ggtitle("Basic boxplot") +
  coord_flip() +
  xlab("") +
  ylab("")

Distribución binomial

La distribución binomial es una distribución de probabilidad discreta. Describe el resultado de \(n\) ensayos independientes en un experimento. Se supone que cada ensayo tiene sólo dos resultados, éxito o fracaso. Si la probabilidad de un ensayo con éxito es \(p\), entonces la probabilidad de tener \(x\) resultados con éxito en un experimento de \(n\) ensayos independientes es la siguiente

\[P(X=x)=\displaystyle{\binom{n}{x}}p^{x}(1-p)^{n-x}\]

• Ejemplo

Supongamos que hay doce preguntas de opción múltiple en un examen de la clase de inglés. Cada pregunta tiene cinco respuestas posibles, y sólo una de ellas es correcta. Encuentra la probabilidad de tener cuatro o menos respuestas correctas si un estudiante intenta responder a cada pregunta al azar.

• Solución

Como sólo una de las cinco respuestas posibles es correcta, la probabilidad de responder correctamente a una pregunta por azar es \(1/5=0.2\). Podemos encontrar la probabilidad de tener exactamente \(4\) respuestas correctas por intentos aleatorios de la siguiente manera.

dbinom(4, size = 12, prob = 0.2) 

## [1] 0.1328756

Para encontrar la probabilidad de tener cuatro o menos respuestas correctas mediante intentos aleatorios, aplicamos la función dbinom con \(x = 0,\dots,4\).

sum(dbinom(0:4, size = 12, prob = 0.2))

## [1] 0.9274445

Como alternativa, podemos utilizar la función de probabilidad acumulada para la distribución binomial pbinom

pbinom(4, size=12, prob=0.2)

## [1] 0.9274445

La probabilidad de que cuatro o menos preguntas sean contestadas correctamente al azar en un cuestionario de opción múltiple de doce preguntas es del \(92.7\%\).

Distribución normal

• La distribución normal está definida por la siguiente función de densidad de probabilidad, donde \(\mu\) es la media de la población y \(\sigma^2\) es la varianza

\[f(x; \mu, \sigma)=\displaystyle{\frac{1}{\sqrt{2\pi}\sigma}}e^{-(x-\mu)^{2}/(2\sigma^{2})},\quad -\infty<x<\infty\]

• Si una variable aleatoria \(X\) sigue la distribución normal, entonces escribimos: \(X\sim N(\mu, \sigma^{2})\). En particular, la distribución normal con \(\mu = 0\) y \(\sigma = 1\) se denomina distribución normal estándar, y se denota como \(N(0, 1)\)

• La distribución normal es importante debido al Teorema del Límite Central, que establece que la población de todas las medias de muestras aleatorias de tamaño n de una población con media \(\mu\) y varianza \(\sigma^2\) se aproxima a una distribución normal con media \(\mu\) y \(\sigma^{2}/n\) cuando \(n\) se aproxima al infinito}

• Si \(X\) es una variable aleatoria normalmente distribuida, con media = \(\mu\) y desviación estándar = \(\sigma\), entonces

\[ \begin{align*} P(X<x_{\max})&=\text{pnorm}(x_{\max},\, \text{mean}=\mu,\, \text{sd}=\sigma,\, \text{lower.tail=TRUE})\\ P(X>x_{\min})&=\text{pnorm}(x_{\min},\, \text{mean}=\mu,\, \text{sd}=\sigma,\, \text{lower.tail=FALSE})\\ P(x_{\min}<X<x_{\max})&=\text{pnorm}(x_{\max},\, \text{mean}=\mu,\, \text{sd}=\sigma,\, \text{lower.tail=TRUE})\\ &\,-\text{pnorm}(x_{\min},\, \text{mean}=\mu,\, \text{sd}=\sigma,\, \text{lower.tail=TRUE}) \end{align*} \]

• Ejemplo

Supongamos que el CI se distribuye normalmente con una media de 100 y una desviación estándar de 15

1. ¿Qué porcentaje de personas tiene un coeficiente intelectual inferior a 125?

pnorm(125, mean = 100, sd = 15, lower.tail=TRUE)

## [1] 0.9522096

2. ¿Qué porcentaje de personas tiene un coeficiente intelectual superior a 110?

pnorm(110, mean = 100, sd = 15, lower.tail=FALSE)

## [1] 0.2524925

3. ¿Qué porcentaje de personas tienen un coeficiente intelectual entre 110 y 125?

pnorm(125, mean = 100, sd = 15, lower.tail=TRUE) - pnorm(110, mean = 100, sd = 15, lower.tail=TRUE)

## [1] 0.2047022

• Obtención de percentiles de una distribución normal con media \(\mu\) y desviación típica \(\sigma\):

• Suponga que quiere encontrar el valor \(x\) que separa el \(k\%\) inferior de los valores de una distribución con media \(\mu\) y desviación estándar \(\sigma\). Denotamos este valor en el texto como \(P_{k}\)

\[P_{k}=\textrm{qnorm}(k (\text{formato decimal}),\, \text{mean}=\mu,\, \text{sd}=\sigma,\, \text{lower.tail=TRUE})\]

• Ejemplo

Supongamos que el CI se distribuye normalmente con una media de 100 y una desviación estándar de 15

1. ¿Qué coeficiente intelectual separa al 25% inferior de los demás? (Encuentra el P25.)

qnorm(.25, mean = 100, sd = 15, lower.tail=TRUE)

## [1] 89.88265

2. ¿Qué coeficiente intelectual separa al 10% superior de los demás? (Encuentra el P90.)

qnorm(.90, mean = 100, sd = 15, lower.tail=TRUE)

## [1] 119.2233

Teorema del límite central

Sea \(X_{1}, X_{2},\dots, X_{n}\) una muestra aleatoria de una distribución con media \(\mu\) y varianza \(\sigma^{2}\). Entonces si \(n\) es suficientemente grande, \(\overline{X}\) tiene aproximadamente una distribución normal con \(\mu_{\overline{X}}=\mu\) y \(\sigma_{\overline{X}}^{2}=\sigma^{2}/n\), y \(T_{\text{o}}\) también tiene aproximadamente una distribución normal con \(\mu_{T_{\text{o}}}=n\mu, \sigma_{T_{\text{o}}}^{2}=n\sigma^{2}\). Mientras más grande es el valor de \(n\), mejor es la aproximación.

• Simulación (TLC): La repetición de una simulación es la clave para obtener información de una simulación, ya que podemos utilizar fácilmente los resúmenes gráficos de los números simulados para el mecanismo que los produce

mu <- 100; sigma <- 16
M <- 4; n <- 16 # cambiar M para más replicas
res <- numeric(M)
for (i in 1:M) {
    res[i] <- mean(rnorm(n, mean=mu, sd=sigma))
}
res

## [1]  97.87882  99.32994 101.36965 103.20213

• Comenzamos haciendo una función para encontrar la puntuación \(z\) de una muestra de \(\overline{x}\) donde centramos por la media de la población y escalamos por la desviación estándar de la población:

zstat <- function(x, mu, sigma) {
  (mean(x) - mu) / (sigma/sqrt(length(x)))
}

• Distribución exponencial: Se dice que \(X\) tiene una distribución exponencial con parámetro \(\lambda(\lambda>0)\) si la función de densidad de probabilidad de \(X\) es

\[f(x; \lambda)= \begin{cases} \lambda e^{-\lambda x}, & x\geq 0\\ 0, & \text{otro caso} \end{cases}\]

• La función dexp(x, rate) nos permitirá visualizar la función de densidad de probabilidad asociada a la distribución exponencial, esta función requiere como argumento de entrada el valor de \(\lambda=\)rate y el valor de \(x\) observado

curve(dexp(x, rate = 1.0), from=0, to=10, col='blue')
curve(dexp(x, rate = 1.5), from=0, to=10, col='red', add=TRUE)
curve(dexp(x, rate = 2.0), from=0, to=10, col='purple', add=TRUE)

legend("topright", 
       legend=c("rate=1.0", "rate=1.5", "rate=2.0"),
       col = c("blue", "red", "purple"), 
       lty=1, cex=1.2)

• El código para generar la distribución exponencial aleatoria en R es rexp(n, rate) donde n se refiere al tamaño de la muestra y rate es el parámetro de la tasa \(\lambda\). La media de la distribución exponencial es \(\mu=1/\lambda\) y la desviación estándar es también \(1/\lambda\). En nuestro ejercicio, \(\lambda\) se fija en 2.0 para todas las simulaciones. Usamos el gráfico de histograma y QQ-plot para estudiar si existe normalidad

par(mfrow=c(1,2))
hist(rexp(10000, rate = 2.0), main=paste("Dist. Exp.:  n =", 10000), col = "gold")
qqnorm(rexp(10000, rate = 2.0), main=paste("Dist. Exp.; n =", 10000), col = "gold")
qqline(rexp(10000, rate = 2.0))
grid(lwd = 2)

• La funciónrexp() nos permite obtener \(n\) observaciones de una distribución exponencial usando la sintaxis rexp(n, rate)

M <- 2000; n <- 7
rate <- 2; mu <- sigma <- 1/rate
res <- replicate(M, {
    x <- rexp(n, rate=rate)
    zstat(x, mu, sigma)
})

par(mfrow=c(1,2))
hist(res, main=paste("Dist. Exp.: M =", M, "; n =", n), col = "blue")
qqnorm(res, main=paste("Dist. Exp.; n =", n), col = "blue")
qqline(res)
grid(lwd = 2)

• El gráfico muestra una curva ascendente indicando una distribución de cola derecha. El fenómeno del límite central no ha eliminado la asimetría inicial de la población cuando \(n=7\)

• Repetimos con \(n = 150\). El gráfico de histograma y QQ-plot muestra ahora que la distribución muestral es ahora aproximadamente normal.

M <- 2000; n <- 150
rate <- 2; mu <- sigma <- 1/rate
res <- replicate(M, {
    x <- rexp(n, rate=rate)
    zstat(x, mu, sigma)
})

par(mfrow=c(1,2))
hist(res, main=paste("Dist. Exp.: M =", M, "; n =", n), col = "red")
qqnorm(res, main=paste("Dist. Exp.; n =", n), col = "red")
qqline(res)
grid(lwd = 2)

Modelado de la tendencia de la serie y componentes estacionales

• En primer lugar, modelaremos la tendencia de la serie mediante la regresión de la serie con la variable de tendencia trend, en la partición de entrenamiento

md_trend <- lm(y ~ trend, data = train)

• Utilizaremos la función de summary para revisar los detalles del modelo:

summary(md_trend)

## 
## Call:
## lm(formula = y ~ trend, data = train)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -547.2 -307.4 -153.2  333.1 1052.6 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1751.0074    52.6435   33.26  < 2e-16 ***
## trend          2.4489     0.4021    6.09 4.86e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 394.4 on 224 degrees of freedom
## Multiple R-squared:  0.1421, Adjusted R-squared:  0.1382 
## F-statistic: 37.09 on 1 and 224 DF,  p-value: 4.861e-09

• Como puede observar en el resultado de la regresión anterior, la cuarta columna representa el valor \(p\) de cada uno de los coeficientes del modelo. El valor \(p\) proporciona la probabilidad de que rechacemos la hipótesis nula dado que es verdadera, o el error de tipo I. Por lo tanto, para el valor \(p\) menor que, el valor umbral, rechazaremos la hipótesis nula con un nivel de significación de \(\alpha\), donde los valores típicos de son 0.1, 0.05, 0.01, etc.

• Utilizaremos md_trend para predecir los valores ajustados y previstos del modelo de tendencia que entrenamos antes

train$yhat <- predict(md_trend, newdata = train)
head(train)

##           ds      y trend seasonal     yhat
## 1 2000-01-01 2510.5     1      Jan 1753.456
## 2 2000-02-01 2330.7     2      Feb 1755.905
## 3 2000-03-01 2050.6     3      Mar 1758.354
## 4 2000-04-01 1783.3     4      Apr 1760.803
## 5 2000-05-01 1632.9     5      May 1763.252
## 6 2000-06-01 1513.1     6      Jun 1765.701

test$yhat <- predict(md_trend, newdata = test)
head(test)

##             ds      y trend seasonal     yhat
## 227 2018-11-01 2709.9   227      Nov 2306.917
## 228 2018-12-01 2993.1   228      Dec 2309.366
## 229 2019-01-01 3399.9   229      Jan 2311.815
## 230 2019-02-01 2999.2   230      Feb 2314.264
## 231 2019-03-01 2899.9   231      Mar 2316.713
## 232 2019-04-01 2201.1   232      Apr 2319.162

• Vamos a crear una función de utilidad que traza las series y la salida del modelo, utilizando el paquete plotly

library(plotly)

plot_lm <- function(data, train, test, title = NULL){
  p <- plot_ly(data = data,
               x = ~ ds,
               y = ~ y,
               type = "scatter",
               mode = "line",
               name = "Actual") %>%
    add_lines(x = ~ train$ds,
              y = ~ train$yhat,
              line = list(color = "red"),
              name = "Fitted") %>%
    add_lines(x = ~ test$ds,
              y = ~ test$yhat,
              line = list(color = "green", dash = "dot", width = 3),
              name = "Forecasted") %>%
    layout(title = title,
           xaxis = list(title = "Year"),
           yaxis = list(title = "Billion Cubic Feet"),
           legend = list(x = 0.05, y = 0.95))
  return(p)
}

• Los argumentos de la función son los siguientes

  • data: Los datos de entrada, un objeto data.frame que sigue la misma estructura que la de USgas_df (incluyendo la variable yhat)
  • train: El conjunto de entrenamiento correspondiente que se utilizó para entrenar el modelo
  • test: Igualmente, el conjunto de pruebas correspondiente que se utilizó para evaluar el modelo de predicción title: El título de la parcela, por defecto, establecido como NULL

• Establezcamos las entradas de la función plot_lm con la salida del modelo

plot_lm(data = USgas_df,
        train = train,
        test = test,
        title = "Predicting the Trend Component of the Series")

• En general, el modelo fue capaz de captar el movimiento general de la tendencia, aunque una tendencia lineal puede no captar la ruptura estructural de la tendencia que se produjo en torno a 2010. Más adelante, en este capítulo, veremos un método avanzado para captar una tendencia no lineal. Por último, pero no menos importante, para el análisis de comparación, queremos medir la tasa de error del modelo tanto en en los conjuntos de entrenamiento y de prueba. La norma usada se define como

\[ \textsf{MAPE} = \frac{1}{M}\sum_{i=1}^{n}\left|\frac{A_{i}-F_{i}}{A_{i}}\right| \]

mape_trend <- c(mean(abs(train$y - train$yhat)/train$y), mean(abs(test$y - test$yhat)/test$y))
mape_trend

## [1] 0.1644088 0.1299951

• El proceso de modelamiento y predicción del componente estacional, el cual sigue el mismo proceso que aplicamos con la tendencia, haciendo una regresión de la serie con la variable estacional seasonal que creamos antes

md_seasonal <- lm(y ~ seasonal, data = train)
summary(md_seasonal)

## 
## Call:
## lm(formula = y ~ seasonal, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -608.98 -162.34  -50.77  148.40  566.89 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2774.28      49.75  55.759  < 2e-16 ***
## seasonalFeb  -297.92      70.36  -4.234 3.41e-05 ***
## seasonalMar  -479.10      70.36  -6.809 9.77e-11 ***
## seasonalApr  -905.28      70.36 -12.866  < 2e-16 ***
## seasonalMay -1088.42      70.36 -15.468  < 2e-16 ***
## seasonalJun -1105.49      70.36 -15.711  < 2e-16 ***
## seasonalJul  -939.35      70.36 -13.350  < 2e-16 ***
## seasonalAug  -914.12      70.36 -12.991  < 2e-16 ***
## seasonalSep -1114.74      70.36 -15.843  < 2e-16 ***
## seasonalOct -1022.21      70.36 -14.527  < 2e-16 ***
## seasonalNov  -797.53      71.33 -11.180  < 2e-16 ***
## seasonalDec  -256.67      71.33  -3.598 0.000398 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 216.9 on 214 degrees of freedom
## Multiple R-squared:  0.7521, Adjusted R-squared:  0.7394 
## F-statistic: 59.04 on 11 and 214 DF,  p-value: < 2.2e-16

• Dado que regresamos la variable dependiente con una variable categórica, el modelo de regresión crea coeficientes para 11 de las 12 categorías, que son las que están incrustadas con los valores de la pendiente. Como se puede ver en el resumen de la regresión del modelo estacional, todos los coeficientes del modelo son estadísticamente significativos. Además, se puede observar que el R-cuadrado ajustado de el modelo estacional es algo mayor con respecto al modelo de tendencia (0.78 frente a 0.1). Antes de trazar el modelo ajustado y los valores de previsión con la función plot_lm, vamos a actualizar los valores de yhat con la función predict

train$yhat <- predict(md_seasonal, newdata = train)
head(train)

##           ds      y trend seasonal     yhat
## 1 2000-01-01 2510.5     1      Jan 2774.279
## 2 2000-02-01 2330.7     2      Feb 2476.358
## 3 2000-03-01 2050.6     3      Mar 2295.179
## 4 2000-04-01 1783.3     4      Apr 1869.000
## 5 2000-05-01 1632.9     5      May 1685.863
## 6 2000-06-01 1513.1     6      Jun 1668.784

test$yhat <- predict(md_seasonal, newdata = test)
head(test)

##             ds      y trend seasonal     yhat
## 227 2018-11-01 2709.9   227      Nov 1976.750
## 228 2018-12-01 2993.1   228      Dec 2517.606
## 229 2019-01-01 3399.9   229      Jan 2774.279
## 230 2019-02-01 2999.2   230      Feb 2476.358
## 231 2019-03-01 2899.9   231      Mar 2295.179
## 232 2019-04-01 2201.1   232      Apr 1869.000

• Ahora podemos utilizar la función plot_lm para visualizar el modelo ajustado y los valores predecidos

plot_lm(data = USgas_df,
        train = train,
        test = test,
        title = "Predicting the Seasonal Component of the Series")

• Como se puede ver en el gráfico anterior, el modelo está haciendo un trabajo bastante bueno para capturar la estructura del patrón estacional de la serie. Sin embargo, se puede observar que la tendencia de la serie no presenta una buen predicción. Antes de añadir tanto la tendencia como los componentes estacionales, puntuaremos el rendimiento del modelo:

mape_seasonal <- c(mean(abs(train$y - train$yhat) / train$y), mean(abs(test$y - test$yhat) / test$y))
mape_seasonal

## [1] 0.08628973 0.21502100

• El alto índice de error en el conjunto de pruebas está relacionado con el componente de tendencia que no se incluyó en el modelo. El siguiente paso es unir los dos componentes en un modelo y pronosticar los valores característicos de la serie

md1 <- lm(y ~ seasonal + trend, data = train)
summary(md1)

## 
## Call:
## lm(formula = y ~ seasonal + trend, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -514.73  -77.17  -17.70   85.80  336.95 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2488.8994    32.6011  76.344  < 2e-16 ***
## seasonalFeb  -300.5392    41.4864  -7.244 7.84e-12 ***
## seasonalMar  -484.3363    41.4870 -11.674  < 2e-16 ***
## seasonalApr  -913.1334    41.4880 -22.010  < 2e-16 ***
## seasonalMay -1098.8884    41.4895 -26.486  < 2e-16 ***
## seasonalJun -1118.5855    41.4913 -26.960  < 2e-16 ***
## seasonalJul  -955.0563    41.4936 -23.017  < 2e-16 ***
## seasonalAug  -932.4482    41.4962 -22.471  < 2e-16 ***
## seasonalSep -1135.6874    41.4993 -27.366  < 2e-16 ***
## seasonalOct -1045.7687    41.5028 -25.198  < 2e-16 ***
## seasonalNov  -808.0016    42.0617 -19.210  < 2e-16 ***
## seasonalDec  -269.7642    42.0635  -6.413 9.05e-10 ***
## trend           2.6182     0.1305  20.065  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 127.9 on 213 degrees of freedom
## Multiple R-squared:  0.9142, Adjusted R-squared:  0.9094 
## F-statistic: 189.2 on 12 and 213 DF,  p-value: < 2.2e-16

• La regresión de la serie \(y\) con la tendencia(trend) y los componentes estacionales (seasonal) juntos proporciona un aumento adicional del R-cuadrado ajustado del modelo, que pasa de 0.78 a 091. Esto puede verse en el gráfico de los resultados del modelo

train$yhat <- predict(md1, newdata = train)
test$yhat <- predict(md1, newdata = test)

plot_lm(data = USgas_df,
        train = train,
        test = test,
        title = "Predicting the Trend and Seasonal Components of the Series")

• Midamos la puntuación MAPE del modelo en las particiones de entrenamiento y de prueba:

mape_md1 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md1

## [1] 0.04774945 0.09143290

• La regresión de las series con la tendencia y los componentes estacionales proporciona una mejora significativamente tanto la calidad del ajuste del modelo como la precisión del mismo. Sin embargo, al observar el gráfico del ajuste del modelo y la predicción, se puede observar que la tendencia del modelo es demasiado lineal, y no tiene la ruptura estructural de la tendencia de la serie.

• Este es el punto en el que la adición de un componente polinómico para el modelo podría proporcionar una mejora adicional de la precisión del modelo. Una técnica sencilla para captar una tendencia no lineal consiste en añadir un componente polinómico a la tendencia de la serie para capturar la curvatura de la tendencia en el tiempo. Utilizaremos el argumento I que nos permite aplicar operaciones matemáticas sobre cualquiera de los objetos de entrada. Por lo tanto, utilizaremos este argumento para añadir un segundo grado del polinomio para la entrada de la tendencia

md2 <- lm(y ~ seasonal + trend + I(trend^2), data = train)
summary(md2)

## 
## Call:
## lm(formula = y ~ seasonal + trend + I(trend^2), data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -468.47  -54.66   -2.21   63.11  294.32 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.635e+03  3.224e+01  81.738  < 2e-16 ***
## seasonalFeb -3.004e+02  3.540e+01  -8.487 3.69e-15 ***
## seasonalMar -4.841e+02  3.540e+01 -13.676  < 2e-16 ***
## seasonalApr -9.128e+02  3.540e+01 -25.787  < 2e-16 ***
## seasonalMay -1.099e+03  3.540e+01 -31.033  < 2e-16 ***
## seasonalJun -1.118e+03  3.540e+01 -31.588  < 2e-16 ***
## seasonalJul -9.547e+02  3.540e+01 -26.968  < 2e-16 ***
## seasonalAug -9.322e+02  3.541e+01 -26.329  < 2e-16 ***
## seasonalSep -1.136e+03  3.541e+01 -32.070  < 2e-16 ***
## seasonalOct -1.046e+03  3.541e+01 -29.532  < 2e-16 ***
## seasonalNov -8.001e+02  3.590e+01 -22.286  < 2e-16 ***
## seasonalDec -2.618e+02  3.590e+01  -7.293 5.95e-12 ***
## trend       -1.270e+00  4.472e-01  -2.840  0.00494 ** 
## I(trend^2)   1.713e-02  1.908e-03   8.977  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 109.1 on 212 degrees of freedom
## Multiple R-squared:  0.9379, Adjusted R-squared:  0.9341 
## F-statistic: 246.1 on 13 and 212 DF,  p-value: < 2.2e-16

• La adición del polinomio de segundo grado al modelo de regresión no supuso una mejora significativa de la bondad de ajuste del modelo. En el otro modelo, como se puede ver en el siguiente gráfico de salida, este simple cambio en la estructura del modelo nos permite capturar la ruptura estructural de la tendencia en el tiempo

train$yhat <- predict(md2, newdata = train)
test$yhat <- predict(md2, newdata = test)

plot_lm(data = USgas_df,
        train = train,
        test = test,
        title = "Predicting the Trend (Polynomial) and Seasonal Components of the Series")

• Como podemos ver en el modelo que sigue a la puntuación MAPE, la precisión del modelo mejoró significativamente al añadir la tendencia polinómica al modelo de regresión, ya que el error en el conjunto de pruebas se redujo del 9.2% al 4.5%

mape_md2 <- c(mean(abs(train$y - train$yhat) / train$y),
mean(abs(test$y - test$yhat) / test$y))
mape_md2

## [1] 0.03688770 0.04212618

library(TSstudio)
data(USgas)
acf(USgas, lag.max = 60)

Valores propios / Varianzas

• Los valores propios miden la cantidad de variación de cada componente principal. Los valores propios son grandes para los primeros PC y pequeños para los siguientes PC. Examinamos los valores propios para determinar el número de componentes principales que hay que considerar. Los valores propios y la proporción de varianza (es decir, la información) retenida por los componentes principales PC pueden extraerse utilizando la función get_eigenvalue()

library("factoextra")

eig.val <- get_eigenvalue(res.pca)
eig.val

##        eigenvalue variance.percent cumulative.variance.percent
## Dim.1   4.1242133        41.242133                    41.24213
## Dim.2   1.8385309        18.385309                    59.62744
## Dim.3   1.2391403        12.391403                    72.01885
## Dim.4   0.8194402         8.194402                    80.21325
## Dim.5   0.7015528         7.015528                    87.22878
## Dim.6   0.4228828         4.228828                    91.45760
## Dim.7   0.3025817         3.025817                    94.48342
## Dim.8   0.2744700         2.744700                    97.22812
## Dim.9   0.1552169         1.552169                    98.78029
## Dim.10  0.1219710         1.219710                   100.00000

• La suma de todos los valores propios da una varianza total de 10. La proporción de la variación explicada por cada valor propio se indica en la segunda columna. Por ejemplo, 4.124 dividido por 10 es igual a 0.4124, es decir, alrededor del 41.24% de la variación se explica por este primer valor propio.

• En nuestro análisis, los tres primeros componentes principales explican el 72% de la variación. Este es un porcentaje aceptablemente grande. Un método alternativo para determinar el número de componentes principales es observar un Scree Plot usando fviz_eig(), que es el gráfico de los valores propios ordenados de mayor a menor.

fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50))

• Un método sencillo para extraer los resultados, para las variables, de la salida de un PCA es utilizar la función get_pca_var(). Esta función proporciona una lista de matrices que contiene todos los resultados de las variables activas (coordenadas, correlación entre variables y ejes, coseno cuadrado y contribuciones)

var <- get_pca_var(res.pca)
var

## Principal Component Analysis Results for variables
##  ===================================================
##   Name       Description                                    
## 1 "$coord"   "Coordinates for the variables"                
## 2 "$cor"     "Correlations between variables and dimensions"
## 3 "$cos2"    "Cos2 for the variables"                       
## 4 "$contrib" "contributions of the variables"

• Los componentes de get_pca_var() se pueden utilizar en el gráfico de variables como sigue
  • var$coord: coordenadas de las variables para crear un gráfico de dispersión
  • var$cos2: representa la calidad de la representación de las variables en el mapa de factores. Se calcula como el cuadrado de las coordenadas: var.cos2 = var.coord * var.coord.
  • var$contrib: contiene las contribuciones (en porcentaje) de las variables a los componentes del principio. La contribución de una variable (var) a un componente principal determinado es (en porcentaje) : (var.cos2 * 100) / (cos2 total del componente).

# Coordinates
var$coord

##                     Dim.1       Dim.2       Dim.3       Dim.4      Dim.5
## X100m        -0.850625692 -0.17939806  0.30155643  0.03357320 -0.1944440
## Long.jump     0.794180641  0.28085695 -0.19054653 -0.11538956  0.2331567
## Shot.put      0.733912733  0.08540412  0.51759781  0.12846837 -0.2488129
## High.jump     0.610083985 -0.46521415  0.33008517  0.14455012  0.4027002
## X400m        -0.701603377  0.29017826  0.28353292  0.43082552  0.1039085
## X110m.hurdle -0.764125197 -0.02474081  0.44888733 -0.01689589  0.2242200
## Discus        0.743209016  0.04966086  0.17652518  0.39500915 -0.4082391
## Pole.vault   -0.217268042  0.80745110  0.09405773 -0.33898477 -0.2216853
## Javeline      0.428226639  0.38610928  0.60412432 -0.33173454  0.1978128
## X1500m        0.004278487  0.78448019 -0.21947068  0.44800961  0.2632527

# Cos2: quality on the factore map
var$cos2

##                     Dim.1        Dim.2       Dim.3        Dim.4      Dim.5
## X100m        7.235641e-01 0.0321836641 0.090936280 0.0011271597 0.03780845
## Long.jump    6.307229e-01 0.0788806285 0.036307981 0.0133147506 0.05436203
## Shot.put     5.386279e-01 0.0072938636 0.267907488 0.0165041211 0.06190783
## High.jump    3.722025e-01 0.2164242070 0.108956221 0.0208947375 0.16216747
## X400m        4.922473e-01 0.0842034209 0.080390914 0.1856106269 0.01079698
## X110m.hurdle 5.838873e-01 0.0006121077 0.201499837 0.0002854712 0.05027463
## Discus       5.523596e-01 0.0024662013 0.031161138 0.1560322304 0.16665918
## Pole.vault   4.720540e-02 0.6519772763 0.008846856 0.1149106765 0.04914437
## Javeline     1.833781e-01 0.1490803723 0.364966189 0.1100478063 0.03912992
## X1500m       1.830545e-05 0.6154091638 0.048167378 0.2007126089 0.06930197

# Contributions to the principal components
var$contrib

##                     Dim.1      Dim.2      Dim.3       Dim.4     Dim.5
## X100m        1.754429e+01  1.7505098  7.3386590  0.13755240  5.389252
## Long.jump    1.529317e+01  4.2904162  2.9300944  1.62485936  7.748815
## Shot.put     1.306014e+01  0.3967224 21.6204325  2.01407269  8.824401
## High.jump    9.024811e+00 11.7715838  8.7928883  2.54987951 23.115504
## X400m        1.193554e+01  4.5799296  6.4876363 22.65090599  1.539012
## X110m.hurdle 1.415754e+01  0.0332933 16.2612611  0.03483735  7.166193
## Discus       1.339309e+01  0.1341398  2.5147385 19.04132022 23.755756
## Pole.vault   1.144592e+00 35.4618611  0.7139512 14.02307063  7.005084
## Javeline     4.446377e+00  8.1086683 29.4531777 13.42963254  5.577615
## X1500m       4.438531e-04 33.4728757  3.8871610 24.49386930  9.878367

Círculo de correlación

• La correlación entre una variable y una componente principal (PC) se utiliza como las coordenadas de la variable en el PC. La representación de las variables difiere del trazado de las observaciones: Las observaciones se representan por sus proyecciones, pero las variables están representadas por sus correlaciones.

var$coord

##                     Dim.1       Dim.2       Dim.3       Dim.4      Dim.5
## X100m        -0.850625692 -0.17939806  0.30155643  0.03357320 -0.1944440
## Long.jump     0.794180641  0.28085695 -0.19054653 -0.11538956  0.2331567
## Shot.put      0.733912733  0.08540412  0.51759781  0.12846837 -0.2488129
## High.jump     0.610083985 -0.46521415  0.33008517  0.14455012  0.4027002
## X400m        -0.701603377  0.29017826  0.28353292  0.43082552  0.1039085
## X110m.hurdle -0.764125197 -0.02474081  0.44888733 -0.01689589  0.2242200
## Discus        0.743209016  0.04966086  0.17652518  0.39500915 -0.4082391
## Pole.vault   -0.217268042  0.80745110  0.09405773 -0.33898477 -0.2216853
## Javeline      0.428226639  0.38610928  0.60412432 -0.33173454  0.1978128
## X1500m        0.004278487  0.78448019 -0.21947068  0.44800961  0.2632527

fviz_pca_var(res.pca, col.var = "black")

Erosión del suelo

df_soil <- read.csv("cmz-48-2030-yields.csv")
View(distinct(df_soil))

df_soil <- select(df_soil, 
                  areasymbol, 
                  soil_erosion = soil.erosion, 
                  watereros, winderos, sci, sciom, scifo, 
                  slope = Slope)
View(df_soil)
head(df_soil)

##   areasymbol soil_erosion watereros winderos  sci sciom scifo slope
## 1      OK009         0.63      0.63     0.00 0.43 -0.21  0.92   0.5
## 2      OK009         0.19      0.19     0.00 0.70  0.37  0.93   0.5
## 3      OK009         1.85      0.74     1.11 0.06 -0.13  0.15   0.5
## 4      OK009         0.81      0.36     0.46 0.33  0.32  0.16   0.5
## 5      OK009         1.56      1.56     0.00 0.30 -0.37  0.92   2.0
## 6      OK009         0.39      0.39     0.00 0.58  0.10  0.93   2.0

df_areasymbol <- group_by(df_soil, areasymbol)
View(df_areasymbol)
head(df_areasymbol)

## # A tibble: 6 x 8
## # Groups:   areasymbol [1]
##   areasymbol soil_erosion watereros winderos   sci sciom scifo slope
##   <chr>             <dbl>     <dbl>    <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 OK009              0.63      0.63     0     0.43 -0.21  0.92   0.5
## 2 OK009              0.19      0.19     0     0.7   0.37  0.93   0.5
## 3 OK009              1.85      0.74     1.11  0.06 -0.13  0.15   0.5
## 4 OK009              0.81      0.36     0.46  0.33  0.32  0.16   0.5
## 5 OK009              1.56      1.56     0     0.3  -0.37  0.92   2  
## 6 OK009              0.39      0.39     0     0.58  0.1   0.93   2

df_pca <- summarize(df_areasymbol, 
          soil_erosion = mean(soil_erosion, na.rm = TRUE),
          watereros = mean(watereros, na.rm = TRUE),
          winderos = mean(winderos, na.rm = TRUE),
          sci = mean(sci, na.rm = TRUE),
          sciom = mean(sciom, na.rm = TRUE),
          scifo = mean(scifo, na.rm = TRUE),
          slope = mean(slope, na.rm = TRUE))

View(df_pca)
str(df_pca)

## tibble [50 x 8] (S3: tbl_df/tbl/data.frame)
##  $ areasymbol  : chr [1:50] "OK009" "OK015" "OK031" "OK033" ...
##  $ soil_erosion: num [1:50] 7.42 7.55 21.03 19.72 12.08 ...
##  $ watereros   : num [1:50] 3.27 5.53 6.07 5.49 4.37 ...
##  $ winderos    : num [1:50] 4.15 2.01 14.96 14.23 7.71 ...
##  $ sci         : num [1:50] -0.191 -0.249 -1.324 -1.231 -0.59 ...
##  $ sciom       : num [1:50] -0.159 -0.276 -0.311 -0.337 -0.238 ...
##  $ scifo       : num [1:50] 0.641 0.641 0.641 0.641 0.641 ...
##  $ slope       : num [1:50] 4.36 5.27 4.57 3.15 4.09 ...

df_pca <- df_pca[, 2:ncol(df_pca)]
df_pca

## # A tibble: 50 x 7
##    soil_erosion watereros winderos     sci    sciom scifo slope
##           <dbl>     <dbl>    <dbl>   <dbl>    <dbl> <dbl> <dbl>
##  1         7.42      3.27     4.15 -0.191  -0.159   0.641  4.36
##  2         7.55      5.53     2.01 -0.249  -0.276   0.641  5.27
##  3        21.0       6.07    15.0  -1.32   -0.311   0.641  4.57
##  4        19.7       5.49    14.2  -1.23   -0.337   0.641  3.15
##  5        12.1       4.37     7.71 -0.590  -0.238   0.641  4.09
##  6         6.46      2.62     3.84 -0.103  -0.126   0.641  3.55
##  7         5.72      2.13     3.58 -0.0124 -0.0470  0.641  3.72
##  8         3.24      1.53     1.71  0.205   0.00862 0.641  2.85
##  9         6.81      4.18     2.63 -0.147  -0.170   0.641  4.24
## 10         5.81      2.76     3.05 -0.0643 -0.159   0.641  2.84
## # ... with 40 more rows

library("FactoMineR")

res_pca <- PCA(df_pca, graph = FALSE)
print(res_pca)

## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 7 variables
## *The results are available in the following objects:
## 
##    name               description                          
## 1  "$eig"             "eigenvalues"                        
## 2  "$var"             "results for the variables"          
## 3  "$var$coord"       "coord. for the variables"           
## 4  "$var$cor"         "correlations variables - dimensions"
## 5  "$var$cos2"        "cos2 for the variables"             
## 6  "$var$contrib"     "contributions of the variables"     
## 7  "$ind"             "results for the individuals"        
## 8  "$ind$coord"       "coord. for the individuals"         
## 9  "$ind$cos2"        "cos2 for the individuals"           
## 10 "$ind$contrib"     "contributions of the individuals"   
## 11 "$call"            "summary statistics"                 
## 12 "$call$centre"     "mean of the variables"              
## 13 "$call$ecart.type" "standard error of the variables"    
## 14 "$call$row.w"      "weights for the individuals"        
## 15 "$call$col.w"      "weights for the variables"

library("factoextra")

eig_val <- get_eigenvalue(res_pca)
eig_val

##         eigenvalue variance.percent cumulative.variance.percent
## Dim.1 4.065112e+00     5.807303e+01                    58.07303
## Dim.2 1.000000e+00     1.428571e+01                    72.35875
## Dim.3 9.087609e-01     1.298230e+01                    85.34104
## Dim.4 8.407550e-01     1.201079e+01                    97.35183
## Dim.5 1.853716e-01     2.648166e+00                   100.00000
## Dim.6 2.648985e-07     3.784264e-06                   100.00000
## Dim.7 8.820093e-10     1.260013e-08                   100.00000

fviz_eig(res_pca, addlabels = TRUE, ylim = c(0, 60))

var_res <- get_pca_var(res_pca)
var_res

## Principal Component Analysis Results for variables
##  ===================================================
##   Name       Description                                    
## 1 "$coord"   "Coordinates for the variables"                
## 2 "$cor"     "Correlations between variables and dimensions"
## 3 "$cos2"    "Cos2 for the variables"                       
## 4 "$contrib" "contributions of the variables"

library("corrplot")

corrplot(var_res$cos2, is.corr=FALSE)

Indicadores económicos

countries <- read.csv('countries_of_the_world.csv', na.string = c("", "NA"))
head(countries)

##           Country                              Region Population Area..sq..mi..
## 1    Afghanistan        ASIA (EX. NEAR EAST)            31056997         647500
## 2        Albania  EASTERN EUROPE                         3581655          28748
## 3        Algeria  NORTHERN AFRICA                       32930091        2381740
## 4 American Samoa  OCEANIA                                  57794            199
## 5        Andorra  WESTERN EUROPE                           71201            468
## 6         Angola  SUB-SAHARAN AFRICA                    12127071        1246700
##   Pop..Density..per.sq..mi.. Coastline..coast.area.ratio. Net.migration
## 1                       48,0                         0,00         23,06
## 2                      124,6                         1,26         -4,93
## 3                       13,8                         0,04         -0,39
## 4                      290,4                        58,29        -20,71
## 5                      152,1                         0,00           6,6
## 6                        9,7                         0,13             0
##   Infant.mortality..per.1000.births. GDP....per.capita. Literacy....
## 1                             163,07                700         36,0
## 2                              21,52               4500         86,5
## 3                                 31               6000         70,0
## 4                               9,27               8000         97,0
## 5                               4,05              19000        100,0
## 6                             191,19               1900         42,0
##   Phones..per.1000. Arable.... Crops.... Other.... Climate Birthrate Deathrate
## 1               3,2      12,13      0,22     87,65       1      46,6     20,34
## 2              71,2      21,09      4,42     74,49       3     15,11      5,22
## 3              78,1       3,22      0,25     96,53       1     17,14      4,61
## 4             259,5         10        15        75       2     22,46      3,27
## 5             497,2       2,22         0     97,78       3      8,71      6,25
## 6               7,8       2,41      0,24     97,35    <NA>     45,11      24,2
##   Agriculture Industry Service
## 1        0,38     0,24    0,38
## 2       0,232    0,188   0,579
## 3       0,101      0,6   0,298
## 4        <NA>     <NA>    <NA>
## 5        <NA>     <NA>    <NA>
## 6       0,096    0,658   0,246

any(is.na(countries))

## [1] TRUE

sum(is.na(countries))

## [1] 110

library(Amelia)

missmap(countries, legend = TRUE, col = c("red", "dodgerblue"))

str(countries)

## 'data.frame':    227 obs. of  20 variables:
##  $ Country                           : chr  "Afghanistan " "Albania " "Algeria " "American Samoa " ...
##  $ Region                            : chr  "ASIA (EX. NEAR EAST)         " "EASTERN EUROPE                     " "NORTHERN AFRICA                    " "OCEANIA                            " ...
##  $ Population                        : int  31056997 3581655 32930091 57794 71201 12127071 13477 69108 39921833 2976372 ...
##  $ Area..sq..mi..                    : int  647500 28748 2381740 199 468 1246700 102 443 2766890 29800 ...
##  $ Pop..Density..per.sq..mi..        : chr  "48,0" "124,6" "13,8" "290,4" ...
##  $ Coastline..coast.area.ratio.      : chr  "0,00" "1,26" "0,04" "58,29" ...
##  $ Net.migration                     : chr  "23,06" "-4,93" "-0,39" "-20,71" ...
##  $ Infant.mortality..per.1000.births.: chr  "163,07" "21,52" "31" "9,27" ...
##  $ GDP....per.capita.                : int  700 4500 6000 8000 19000 1900 8600 11000 11200 3500 ...
##  $ Literacy....                      : chr  "36,0" "86,5" "70,0" "97,0" ...
##  $ Phones..per.1000.                 : chr  "3,2" "71,2" "78,1" "259,5" ...
##  $ Arable....                        : chr  "12,13" "21,09" "3,22" "10" ...
##  $ Crops....                         : chr  "0,22" "4,42" "0,25" "15" ...
##  $ Other....                         : chr  "87,65" "74,49" "96,53" "75" ...
##  $ Climate                           : chr  "1" "3" "1" "2" ...
##  $ Birthrate                         : chr  "46,6" "15,11" "17,14" "22,46" ...
##  $ Deathrate                         : chr  "20,34" "5,22" "4,61" "3,27" ...
##  $ Agriculture                       : chr  "0,38" "0,232" "0,101" NA ...
##  $ Industry                          : chr  "0,24" "0,188" "0,6" NA ...
##  $ Service                           : chr  "0,38" "0,579" "0,298" NA ...

head(as.character(countries$Agriculture))

## [1] "0,38"  "0,232" "0,101" NA      NA      "0,096"

for (i in 3:length(names(countries))){
    countries[,i] <- as.numeric(gsub(",",'.',(sapply(countries[,i], as.character))))
}

str(countries)

## 'data.frame':    227 obs. of  20 variables:
##  $ Country                           : chr  "Afghanistan " "Albania " "Algeria " "American Samoa " ...
##  $ Region                            : chr  "ASIA (EX. NEAR EAST)         " "EASTERN EUROPE                     " "NORTHERN AFRICA                    " "OCEANIA                            " ...
##  $ Population                        : num  31056997 3581655 32930091 57794 71201 ...
##  $ Area..sq..mi..                    : num  647500 28748 2381740 199 468 ...
##  $ Pop..Density..per.sq..mi..        : num  48 124.6 13.8 290.4 152.1 ...
##  $ Coastline..coast.area.ratio.      : num  0 1.26 0.04 58.29 0 ...
##  $ Net.migration                     : num  23.06 -4.93 -0.39 -20.71 6.6 ...
##  $ Infant.mortality..per.1000.births.: num  163.07 21.52 31 9.27 4.05 ...
##  $ GDP....per.capita.                : num  700 4500 6000 8000 19000 1900 8600 11000 11200 3500 ...
##  $ Literacy....                      : num  36 86.5 70 97 100 42 95 89 97.1 98.6 ...
##  $ Phones..per.1000.                 : num  3.2 71.2 78.1 259.5 497.2 ...
##  $ Arable....                        : num  12.13 21.09 3.22 10 2.22 ...
##  $ Crops....                         : num  0.22 4.42 0.25 15 0 0.24 0 4.55 0.48 2.3 ...
##  $ Other....                         : num  87.7 74.5 96.5 75 97.8 ...
##  $ Climate                           : num  1 3 1 2 3 NA 2 2 3 4 ...
##  $ Birthrate                         : num  46.6 15.11 17.14 22.46 8.71 ...
##  $ Deathrate                         : num  20.34 5.22 4.61 3.27 6.25 ...
##  $ Agriculture                       : num  0.38 0.232 0.101 NA NA 0.096 0.04 0.038 0.095 0.239 ...
##  $ Industry                          : num  0.24 0.188 0.6 NA NA 0.658 0.18 0.22 0.358 0.343 ...
##  $ Service                           : num  0.38 0.579 0.298 NA NA 0.246 0.78 0.743 0.547 0.418 ...

countries$Climate = ifelse(is.na(countries$Climate), 'Unknown', countries$Climate)
countries$Climate = as.factor(countries$Climate)

num_cols = c(3:20)
num_cols = num_cols[num_cols != 15] 

for (index in num_cols)
{countries[,index] = ifelse(is.na(countries[,index]),ave(countries[,index], 
                    FUN =function(x) mean(x, na.rm=TRUE)), countries[,index]) }

missmap(countries, legend = TRUE, col = c("red", "dodgerblue"))