Normal cemento Bogota

R Markdown

library(readxl)
Base_de_datos <- read_excel("C:/Users/USUARIO/Desktop/Base de datos.xlsx")
View(Base_de_datos)

c <- Base_de_datos

Visualizar una distribucion normal

head(c)

## # A tibble: 6 × 4
##   Periodo Ciudad `Total area licenciada metros cuadrados` Despachos de cemento…¹
##     <dbl> <chr>                                     <dbl>                  <dbl>
## 1    2005 Bogota                                   333345                124122 
## 2    2005 Bogota                                   291808                117670.
## 3    2005 Bogota                                   328718                125706.
## 4    2005 Bogota                                   234104                133875 
## 5    2005 Bogota                                   379570                119691 
## 6    2005 Bogota                                   347630                132501.
## # ℹ abbreviated name: ¹​`Despachos de cemento gris toneladas`

Generamos el histograma de la variable “DESPACHO DE CEMENTO GRIS (TON)” y superponemos la curva de densidad de una v.a Normal con media y desviacion típica iguales a la de la varible “DESPACHO DE CEMENTO GRIS (TON)”:

Base de datos, Queremos ver si se distribuye normalamente, para poder predecir resultados de la variable de estudio.

Analisis por medio de histogramas.

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.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

 ggplot(data = c, mapping = aes(x =c$`Despachos de cemento gris toneladas`)) + 
  
  geom_histogram(aes(y =..density..), color="black", fill="bisque3", position = 'identity', bins = 45) +
  
  stat_function(
    fun = dnorm, 
    args = list(mean = mean(c$`Despachos de cemento gris toneladas`), sd = sd(c$`Despachos de cemento gris toneladas`)), 
    lwd = 1, 
    col = 'red2' ) +
  
  scale_y_continuous( n.breaks = 15 )+
  scale_x_continuous(n.breaks = 15)+
  labs(x = "DESPACHOS DE CEMENTO GRIS (TON)", y = "Frecuencia Relativa")

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Use of `` c$`Despachos de cemento gris toneladas` `` is discouraged.
## ℹ Use `Despachos de cemento gris toneladas` instead.
## Use of `` c$`Despachos de cemento gris toneladas` `` is discouraged.
## ℹ Use `Despachos de cemento gris toneladas` instead.

Segun mi criterio, el perfil del histograma y la curva de densidad normal no son aproximadamente iguales, hay demasiada diferencia como para considerar que X=DESPACHOS DE CEMENTO GRIS (TON) tiene una distribucion significativamente normal. Por lo que no habria evidencia estadistica significativa como para aceptar que X, de muestra tiene una distribucion normal.

s=sd(c$`Despachos de cemento gris toneladas`)
n=length(c$`Despachos de cemento gris toneladas`)
max=max(c$`Despachos de cemento gris toneladas`)
min=min(c$`Despachos de cemento gris toneladas`)

#Amplitud teorica
At = 3.5*s*n^-(1/3) 

#Numero de intervalos
lambda = ceiling((max-min)/At) 

#Amplitud de los intervalos
A = ceiling((max-min)/lambda) 

#Primer extremo
L1 = min

#Vector con los extremos de los intervalos
L = L1 + 0:lambda * A  

library(tidyverse)

ggplot(data = c, mapping = aes(x = c$`Despachos de cemento gris toneladas`)) + 
  
  geom_histogram(aes(y =..density..), color="black", fill="bisque3", position = 'identity', breaks=L) +
  
  stat_function(
    fun = dnorm, 
    args = list(mean = mean(c$`Despachos de cemento gris toneladas`), sd = sd(c$`Despachos de cemento gris toneladas`)), 
    lwd = 1.5, 
    col = 'red2' ) +
  
  scale_y_continuous( n.breaks = 15 )+
  scale_x_continuous(n.breaks = 15)+
  labs(x = "DESPACHOS DE CEMENTO GRIS (TON)", y = "Frecuencia Relativa")

## Warning: Use of `` c$`Despachos de cemento gris toneladas` `` is discouraged.
## ℹ Use `Despachos de cemento gris toneladas` instead.
## Use of `` c$`Despachos de cemento gris toneladas` `` is discouraged.
## ℹ Use `Despachos de cemento gris toneladas` instead.

reducir el numero de intervalos el perfil del histograma se aproxima mas a la curva de densidad normal.

Veamos ahora cual sería el resultado si simulamos valores de un v.a. Normal con media y desviación típica igual a la media y desviación tipica de la variable “DESPACHOS DE CEMENTO GRIS (TON)”, y representamos su histograma, y ademas superponemos la curva normal como hicimos antes.

set.seed(123)

df<- tibble(Variable_Normal_Simulada = rnorm(n=200000, mean=mean(c$`Despachos de cemento gris toneladas`), sd=sd(c$`Despachos de cemento gris toneladas`)) )

knitr::kable(head(df) , align="c")

Variable_Normal_Simulada
115062.2
121144.9
154088.7
126682.3
127764.7
156968.1

##Gráfico Q-Q-Norm en R Investigan….

ggplot(c, aes(sample =  c$`Despachos de cemento gris toneladas`) ) +

  stat_qq( color="red2") + stat_qq_line(color="black", size=0.6) +
    
    scale_y_continuous( n.breaks = 15 )+
    scale_x_continuous(n.breaks =20)+
    xlab("Cuantiles teoricos N(0,1)")+
    ylab("Cuantiles variable peso Media ")  

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Use of `` c$`Despachos de cemento gris toneladas` `` is discouraged.
## ℹ Use `Despachos de cemento gris toneladas` instead.
## Use of `` c$`Despachos de cemento gris toneladas` `` is discouraged.
## ℹ Use `Despachos de cemento gris toneladas` instead.

El analisis de Normalidad basado en los graficos Q-Q-Norm se basa en esta última conclusión. El criterio antes expuesto tiene su fundamento teorico en este último hecho.

Como consideramos que el entendimiento de los graficos Q-Q-Norm es algo complicado, al menos en un principio, se va a desarrollar otro grafico mas intuitivo, a mi juicio, que tambien puede usarse para realizar un análisis de normalidad.

##Grafico de Cuantiles Comparados

cuantiles_normales <- quantile(rnorm(200000 , mean=mean(c$`Despachos de cemento gris toneladas`), sd=sd(c$`Despachos de cemento gris toneladas`)) , seq(0, 1 , by=0.05) )

cuantiles_variable <- quantile( c$`Despachos de cemento gris toneladas` ,  seq(0, 1 , by=0.05))

orden_cuantiles <- c(seq(0, 1 , by=0.05), seq(0, 1 , by=0.05) )
cuantiles <- c(cuantiles_normales , cuantiles_variable)

dff<-tibble( orden_cuantiles , cuantiles, tipo=c(rep("Normal", length(cuantiles_normales)) , rep("Variable", length(cuantiles_variable))) )

knitr::kable(head(df) , align="c")

Variable_Normal_Simulada
115062.2
121144.9
154088.7
126682.3
127764.7
156968.1

ggplot(data = dff , aes(x = orden_cuantiles, y = cuantiles, color = tipo))+
  
  geom_point( ) +
  
  scale_y_continuous( n.breaks = 15 )+
    scale_x_continuous(n.breaks =15)+
    xlab("Orden del cuantil")+
    ylab("Cuantiles") +
  scale_color_manual(breaks = c("Normal", "Variable"),
                        values=c("red", "black" ))

Como puede verse hay bastante similitud entre los cuantiles de la variable de interés y los de la variable normal simulada con media y desviacion tipica igual a la de la variable de interés.

A=seq(0, 1 , by=0.005) (Se representarán muchos mas cuantiles en este caso)

cuantiles_normales <- quantile(rnorm(200000 , mean=mean(c$`Despachos de cemento gris toneladas`), sd=sd(c$`Despachos de cemento gris toneladas`)) , seq(0, 1 , by=0.05) )

cuantiles_variable <- quantile(c$`Despachos de cemento gris toneladas` ,  seq(0, 1 , by=0.05))

orden_cuantiles <- c(seq(0, 1 , by=0.05), seq(0, 1 , by=0.05) )

cuantiles <- c(cuantiles_normales , cuantiles_variable)

df<-tibble( orden_cuantiles , cuantiles, tipo=c(rep("Normal", length(cuantiles_normales)) , rep("Variable", length(cuantiles_variable))) )

knitr::kable(head(df) , align="c")

orden_cuantiles	cuantiles	tipo
0.00	43720.93	Normal
0.05	95240.73	Normal
0.10	101859.49	Normal
0.15	106378.25	Normal
0.20	109901.39	Normal
0.25	112965.33	Normal

ggplot(data = df , aes(x = orden_cuantiles, y = cuantiles, color = tipo))+
  
  geom_point( ) +
  
  scale_y_continuous( n.breaks = 15)+
    scale_x_continuous(n.breaks =15)+
    xlab("Orden del cuantil")+
    ylab("Cuantiles") +
  scale_color_manual(breaks = c("Normal", "Variable"),
                        values=c("red", "black" ))

Como puede verse hay bastante similitud entre los cuantiles de la variable de interés y los de la variable normal simulada con media y desviacion tipica igual a la de la variable de interés. Por lo que segun este metodo podria considerarse que la distribucion de peso es significativamente Normal, y con ello aceptarlo a nivel poblacional.

##Gráfico de Comparación de Frecuencias Relativas en R

Tabla_Frecuencias  <- function(X, V1 ,p1 ){

tabla<- matrix(NA, nrow=p1+1, ncol= 4 )
X<- cbind(X, V1)

for(i in 0:p1){

  tabla[i+1, 1 ] <-  dim( ( X  %>% filter(V1==i ) ) )[1]  
}

for(i in 0:p1){
   
  tabla[i+1, 2 ] <- ( dim( ( X  %>% filter(V1==i ) ) )[1] / dim(X)[1] )
}

for(i in 0:p1){

  tabla[i+1, 3] <- cumsum(tabla[1:(i+1) ,1])[i+1]
}

for(i in 0:p1){
   
  tabla[i+1, 4] <- cumsum(tabla[1:(i+1) ,2])[i+1]
}

 rownames(tabla)<-0:p1
 colnames(tabla)<-c("F.Abs","F.Rel", "F.Abs.Cum","F.Rel.Cum")
 
 tabla <- as.data.frame(tabla)
 
 return(tabla) 
}

DESPACHO <- c$`Despachos de cemento gris toneladas`

s=sd(DESPACHO)
n=length(DESPACHO)
max=max(DESPACHO)
min=min(DESPACHO)

#Numero de intervalos --> 45

#Amplitud de los intervalos
A = ceiling((max-min)/45) 

#Primer extremo
L1 = min - 5

#Vector con los extremos de los intervalos
( L = L1 + 0:45 * A  )

##  [1]   9728.558  13320.558  16912.558  20504.558  24096.558  27688.558
##  [7]  31280.558  34872.558  38464.558  42056.558  45648.558  49240.558
## [13]  52832.558  56424.558  60016.558  63608.558  67200.558  70792.558
## [19]  74384.558  77976.558  81568.558  85160.558  88752.558  92344.558
## [25]  95936.558  99528.558 103120.558 106712.558 110304.558 113896.558
## [31] 117488.558 121080.558 124672.558 128264.558 131856.558 135448.558
## [37] 139040.558 142632.558 146224.558 149816.558 153408.558 157000.558
## [43] 160592.558 164184.558 167776.558 171368.558

#Calculamos las frecuencias relativas
set.seed(666)

variable_normal <- rnorm(200000, mean = mean(DESPACHO), sd = sd(DESPACHO))

a<-(cut(variable_normal, breaks = L , include.lowest = T,  labels=F))

variable_normal_categorizada <- a-1

df<- data.frame(variable_normal_categorizada)

Frecuencias_Relativas_Normal <- Tabla_Frecuencias(df, variable_normal_categorizada, p1=44)$F.Rel

df2<- data.frame(Frecuencias_Relativas_Normal, intervalo=1:45)

##########################################################################################################

b<-(cut(DESPACHO, breaks = L , include.lowest = T,  labels=F))

DESPACHO_categorizada <- b-1

df3<- data.frame(DESPACHO_categorizada)

Frecuencias_Relativas_DESPACHO <- Tabla_Frecuencias(df3, DESPACHO_categorizada, p1=44)$F.Rel

##########################################################################################################

df_completo <- data.frame( Frecuencias_Relativas=c(Frecuencias_Relativas_DESPACHO, Frecuencias_Relativas_Normal), intervalo=rep(1:45, 2), tipo=c(rep("Variable de interes",45),rep("Variable Normal",45)))

library(tidyverse)

ggplot(data = df_completo, aes(y=Frecuencias_Relativas, x=intervalo, color=tipo)) +

  geom_point()+
  
  geom_line()+

  scale_y_continuous( n.breaks = 15 )+
  scale_x_continuous(n.breaks = 15)+
  labs(x = "Intervalos", y = "Frecuencia Relativa")

En el grafico se pueden ver las frecuencias relativas de la variable cuantitativa de interes DESPACHO y la variable Normal simulada , en los mismos intervalos. Esto permite analizar visualmente la disparidad o semejanza de las frecuencias relativas en el caso de la variable original frente al caso de que la variable tuviese distribucion normal.

Queda a ojo del analista decidir si las frecuencias relativas de la variable de interes y la variable normal estan suficientemente proximas como para considerarlas aproximadamente iguales y con ello considerar que la variable de interes tiene una distribucion significativamente normal.

##Métodos Analíticos:

Coeficiente de Asimetria

Asimetria en R

library(moments)

skewness(c$`Despachos de cemento gris toneladas`)

## [1] -1.248709

[0, 1.5]

Podemos concluir que hay evidencia estadisitica de que la distribucion de frecuencias de la variable de interes es una distribucion aproximadamente Normal con media y varianza igual a la de la variable de interes.

Si simulamos una variable normal se obtiene el siguiente resultado:

library(moments)

set.seed(666)

skewness(rnorm(n=200000, mean=mean(c$`Despachos de cemento gris toneladas`), sd=sd(c$`Despachos de cemento gris toneladas`))) 

## [1] -0.006374078

Como era de esperar el valor obtenido es prácticamente 0, lo que indica que la distribucion de frecuencias de la variable normal simulada es efectivamente una distribucion Normal.

##Curtosis en R

library(moments)

kurtosis(c$`Despachos de cemento gris toneladas`)

## [1] 9.575457

library(moments)

kurtosis(rnorm(n=15000, mean=mean(c$`Despachos de cemento gris toneladas`), sd=sd(c$`Despachos de cemento gris toneladas`))) 

## [1] 2.910792

Como era de espera la curtosis de la variable normal simulada es aproximadamente 3.

#Contrastes de hipótesis

A continuacion vamos a exponer los contrastes de hipotesis mas empleados para analizar la normalidad de una variable cuantitativa.

Un contraste de hipotesis es un procedimiento estadistico para rechazar una hipotesis basandose en datos/evidencia.

#Test de Shapiro-Wilk en R

shapiro.test( x = c$`Despachos de cemento gris toneladas` )

## 
##  Shapiro-Wilk normality test
## 
## data:  c$`Despachos de cemento gris toneladas`
## W = 0.93197, p-value = 9.434e-09

Como pvalor=9.434e-09 < α=0.05, Se rechaza la hipotesis nula

para un nivel de significacion de 0.05, no hay suficiente evidencia estadistica en los datos como para rechazar la hipotesis de que la variable de interes tenga una distribución de frecuencias Normal con media y desviacion tipica igual a la de la variable de interés.

#Test de Kolmogorov-Smirnov-Lilliefors en R

library(nortest)

lillie.test(x = c$`Despachos de cemento gris toneladas` )

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  c$`Despachos de cemento gris toneladas`
## D = 0.053194, p-value = 0.1216

Como pvalor=0.1216>α=0.05, para un nivel de significacion de 0.05, no hay suficiente evidencia estadistica en los datos como para rechazar la hipotesis de que la variable de interes tenga una distribución de frecuencias Normal con media y desviacion tipica igual a la de la variable de interés.

#Comparacion de funciones de distribución en R

Detras del test de Kolmogorov-Smirnov-Lilliefors esta la comparacion entre la funcion de distribucion de la variable de interes y la de una v.a. Normal con media y desviacion tipica igual a la de la variable de interés.

En el fondo el grafico de comparación de cuantiles expuesto en la seccion de métodos gráficos es la superposicion de los graficos de las funciones de distribucion empiricas de la variable de interes y la v.a. Normal con media y desviacion tipica igual a la de interés.

Podemos obtener con R el grafico de la funcion de distribucion empirica de una variable del siguiente modo:

despacho_medio=c$`Despachos de cemento gris toneladas`

plot(ecdf(despacho_medio))

m=mean(c$`Despachos de cemento gris toneladas`)