#UNIVRSIDAD CENTRAL DEL ECUADOR
#Tema:Estadística Inferencial-Variable Ordinal
#AUTOR:Grupo 6
#CARGA DE DATOS
setwd("/cloud/project")
datos <- read.csv("DATOS.csv", header = TRUE, sep = ";" , dec = ".")
str(datos)
## 'data.frame': 10190 obs. of 17 variables:
## $ Distrito_edit : chr "1" "1" "1" "1" ...
## $ Year_edit_Fecha_del_derrame : int 2013 2013 2013 2013 2013 2013 2013 2013 2013 2013 ...
## $ Mes_edit_Fecha_del_derrame : int 6 3 4 4 6 6 3 9 10 6 ...
## $ Categoria_Instalaciones : chr "Instalacion fija" "Pozos" "Pozos" "Pozos" ...
## $ Operacion_general : chr "Produccion" "Otro" "Produccion" "Produccion" ...
## $ Categoria_Fuente : chr NA "Tanques/Almacenamiento" "Lineas/Tuberias" "Infraestructura Fija" ...
## $ Grupo_causas_probable : chr NA "Afectaciones externas" "Factores humanos" "Problemas tecnicos" ...
## $ Liberacion_petroleo_crudo_edicion : num 0 0 0 0 0 ...
## $ Edicion_recuperacion_petroleo_crudo : num NA 0 0 0 0 0 0 0 0 NA ...
## $ Volumen_liberado_Cond_Final : num 0 0 0 10 0 0 0 1 0 0 ...
## $ Liberacion_agua_de_produccion_edicion: num 6720 3780 5040 420 10920 ...
## $ Liberacion_volumen_gas : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Volumen_condensado_recuperado : num NA 0 0 1 0 0 0 0 0 NA ...
## $ Edicion_Recuperacion_agua_producida : num NA 420 4620 0 10920 ...
## $ Derrame_sobre_agua_limpio : chr "NO" "NO" "NO" "NO" ...
## $ Estado_general : chr "Observaciones tecnicas" NA NA NA ...
## $ Codigo_area : int 1 1 1 1 1 1 1 1 1 3 ...
#VARIABLE ORDINAL 1:CÓDIGO DE ÁREA ----------------------------------------------------
#Extraer la variable
Codigo_area <- datos$Codigo_area
# Filtrar valores NA, -1 y vacíos
Codigo_area <- Codigo_area[!is.na(Codigo_area) & Codigo_area != -1 & Codigo_area != ""]
# Crear la tabla de frecuencias
TDF_Codigo_area <- table(Codigo_area)
Tabla <- as.data.frame(TDF_Codigo_area)
colnames(Tabla) <- c("Codigo_area", "ni")
# Calcular la frecuencia relativa (hi)
Tabla$hi <- (Tabla$ni / sum(Tabla$ni)) * 100
Tabla$Codigo_area <- as.character(Tabla$Codigo_area)
# Agregar fila total
fila_total <- data.frame(Codigo_area = "Total",
ni = sum(Tabla$ni),
hi = sum(Tabla$hi),
stringsAsFactors = FALSE)
# Agregar fila total a la tabla
Tabla <- rbind(Tabla, fila_total)
# Mostrar tabla final
print(Tabla)
## Codigo_area ni hi
## 1 1 7285 92.34376981
## 2 2 366 4.63937128
## 3 3 78 0.98871847
## 4 4 155 1.96476106
## 5 5 3 0.03802763
## 6 7 2 0.02535176
## 7 Total 7889 100.00000000
#Graficas
# Ajustar márgenes para etiquetas largas
par(mar = c(8, 4, 3, 2))
# Gráfico de barras (solo sin la fila Total)
barplot(TDF_Codigo_area,
main = "Gráfica No.2:
Distribución de los codigos Area ",
xlab = "",
ylab = "Frecuencia Absoluta",
col = "aquamarine",
las = 3, # Rotar etiquetas en eje x
cex.names = 0.8,
cex.axis = 0.9,
cex.main = 1.2,
ylim = c(0, max(TDF_Codigo_area) * 1.2),
border = "black")
# Eliminar la fila "Total" para graficar
Tabla <- Tabla[Tabla$Codigo_area != "Total", ] # Eliminar fila "Total" para graficar
# Convertir Codigo_area a factor con orden correcto
Tabla$Codigo_area <- factor(Tabla$Codigo_area, levels = Tabla$Codigo_area) # asegurar orden
# Gráfico de barras de frecuencia relativa (%)
barplot(Tabla$hi,
names.arg = Tabla$Codigo_area,
main = "Gráfica No.2: Distribución Relativa de Codigo_area",
xlab = "Código de Área",
ylab = "Frecuencia Relativa (%)",
col = "aquamarine2",
las = 3,
cex.names = 0.8,
cex.axis = 0.9,
cex.main = 1.2,
ylim = c(0, max(Tabla$hi) * 1.2),
border = "black")
# Excluir fila "Total" y convertir a numérico
Tabla <- Tabla[Tabla$Codigo_area != "Total", ]
Tabla$Codigo_area <- as.numeric(Tabla$Codigo_area)
# Parámetros estimados para la binomial
n <- max(Tabla$Codigo_area) # Número máximo de la variable Codigo_area
media <- sum(Tabla$Codigo_area * Tabla$ni) / sum(Tabla$ni) # Media ponderada
p <- media / n # Estimación de la probabilidad
# Mostrar el valor de p
p
## [1] 0.1879833
# Frecuencias esperadas
Fe2 <- dbinom(Tabla$Codigo_area, size = n, prob = p) * sum(Tabla$ni) # Frecuencias esperadas
Fo2 <- Tabla$ni # Frecuencias observadas
# Agregar columna de frecuencias esperadas a la tabla
Tabla$esperados <- Fe2
# Mostrar la tabla con las frecuencias esperadas
print(Tabla)
## Codigo_area ni hi esperados
## 1 1 7285 92.34376981 3141.3569861
## 2 2 366 4.63937128 1818.0738411
## 3 3 78 0.98871847 561.1829548
## 4 4 155 1.96476106 97.4361133
## 5 5 3 0.03802763 9.0226510
## 6 6 2 0.02535176 0.3481265
# Gráfico Experimental vs Teórica (Binomial)
barplot(rbind(Fo2, Fe2),
beside = TRUE,
col = c("aquamarine", "aquamarine4"),
legend.text = c("Experimental", "Teórica"),
args.legend = list(x = "topright"),
main = "Gráfica No.3: Distribución Experimental vs Teórica (Binomial)",
names.arg = Tabla$Codigo_area, # Cambié Puntaje por Codigo_area
ylab = "Frecuencia")
#Test Chi-cuadrado
x2_2 <- sum(((Fo2 - Fe2)^2) / Fe2)
k <- length(Fo2) # número de clases
gl_2 <- k - 1 - 1 # grados de libertad: k - 1 - 1
vc_2 <- qchisq(0.95, df = gl_2)
cat("Chi-cuadrado calculado (x2_2):", x2_2, "\n")
## Chi-cuadrado calculado (x2_2): 7087.365
cat("Valor crítico (vc_2):", vc_2, "\n")
## Valor crítico (vc_2): 9.487729
if (x2_2 < vc_2) {
cat("Resultado: No se rechaza H0 (la binomial se ajusta bien)\n")
} else {
cat("Resultado: Se rechaza H0 (la binomial no se ajusta bien)\n")
}
## Resultado: Se rechaza H0 (la binomial no se ajusta bien)
#Correlación numérica
correlacion <- cor(Fo2, Fe2)
cat("Correlación Experimental vs Teórico:", correlacion, "\n")
## Correlación Experimental vs Teórico: 0.8635876
#Gráfico de Correlación
plot(Fe2, Fo2,
main = "Gráfica No.4:
Correlación entre Frecuencias Experimental y Teórica ",
xlab = "Frecuencia Experimental(Fe2)",
ylab = "Frecuencia Teórica (Fo2)",
pch = 19,
col = "green3",
xlim = c(min(Fe2), max(Fe2)),
ylim = c(min(Fo2), max(Fo2)))
abline(lm(Fo2 ~ Fe2), col = "green4", lwd = 2)
grid()
# Mostrar r en el gráfico
text(x = max(Fe2)*0.7,
y = max(Fo2)*0.9,
labels = paste("r =", round(correlacion, 3)),
col = "darkred", cex = 1.2)
#Modelo de Poisson
# Excluir fila "Total" y convertir a numérico
Tabla <- Tabla[Tabla$Codigo_area != "Total", ]
Tabla$Codigo_area <- as.numeric(Tabla$Codigo_area)
Tabla$ni <- as.numeric(Tabla$ni)
# Estimación de lambda para Poisson
lambda <- sum(Tabla$Codigo_area * Tabla$ni) / sum(Tabla$ni)
# Mostrar el valor de lambda estimado
cat("Lambda estimado para Poisson:", lambda, "\n")
## Lambda estimado para Poisson: 1.1279
# Rango de Codigo_area observados
x_vals <- Tabla$Codigo_area
# Calcular probabilidades Poisson
probs_pois <- dpois(x_vals, lambda = lambda)
# Frecuencias esperadas según Poisson
Fe_pois <- probs_pois * sum(Tabla$ni)
# Añadir las frecuencias esperadas a la tabla
Tabla$esperados_pois <- Fe_pois
# Mostrar la tabla con las frecuencias esperadas
print(Tabla)
## Codigo_area ni hi esperados esperados_pois
## 1 1 7285 92.34376981 3141.3569861 2880.39352
## 2 2 366 4.63937128 1818.0738411 1624.39736
## 3 3 78 0.98871847 561.1829548 610.71905
## 4 4 155 1.96476106 97.4361133 172.20744
## 5 5 3 0.03802763 9.0226510 38.84654
## 6 6 2 0.02535176 0.3481265 7.30250
# Gráfico Experimental vs Teórico (Poisson)
barplot(rbind(Tabla$ni, Tabla$esperados_pois),
beside = TRUE,
col = c("cyan", "cyan4"),
legend.text = c("Experimental", "Teórico (Poisson)"),
args.legend = list(x = "topright"),
names.arg = Tabla$Codigo_area, # Cambié Puntaje por Codigo_area
main = "Gráfica No.5: Distribución Experimental vs Teórico (Poisson)",
xlab = "Código de Área",
ylab = "Frecuencia")
#Test de Chi-cuadrado (Poisson)
Fo_pois <- Tabla$ni
Fe_pois <- Tabla$esperados_pois
x2_pois <- sum((Fo_pois - Fe_pois)^2 / Fe_pois)
k <- length(Fo_pois)
gl_pois <- k - 1 - 1 # k - 1 - parámetros estimados (λ)
vc_pois <- qchisq(0.95, df = gl_pois)
cat("Chi-cuadrado Poisson:", x2_pois, "\n")
## Chi-cuadrado Poisson: 8213.576
cat("Valor crítico (Poisson):", vc_pois, "\n")
## Valor crítico (Poisson): 9.487729
if (x2_pois < vc_pois) {
cat("Resultado: No se rechaza H0 (la Poisson se ajusta bien)\n")
} else {
cat("Resultado: Se rechaza H0 (la Poisson no se ajusta bien)\n")
}
## Resultado: Se rechaza H0 (la Poisson no se ajusta bien)
#Correlación Experimental vs Teórico (Poisson)
cor_pois <- cor(Fo_pois, Fe_pois)
cat("Correlación Experimental vs Teórico (Poisson):", round(cor_pois, 4), "\n")
## Correlación Experimental vs Teórico (Poisson): 0.8705
plot(Fe_pois, Fo_pois,
main = "Gráfica No.6: Correlación de
Frecuencia Observada vs Esperada (Poisson)",
xlab = "Frecuencia Teórico (Poisson)",
ylab = "Frecuencia Experimental",
pch = 19,
col = "deepskyblue",
xlim = c(min(Fe_pois), max(Fe_pois)),
ylim = c(min(Fo_pois), max(Fo_pois)))
abline(lm(Fo_pois ~ Fe_pois), col = "deepskyblue4", lwd = 2)
grid()
# Agregar valor de la correlación r
text(x = max(Fe_pois)*0.7,
y = max(Fo_pois)*0.9,
labels = paste("r =", round(cor(Fo_pois, Fe_pois), 3)),
col = "darkred", cex = 1.2)
