#=====================================================
# REGRESIÓN POLINOMIAL SIMPLE CON TRATAMIENTO DE DATOS
# X = PHI_10
# Y = CLAY_PCT
#=====================================================
# 1. CARGA DE DATOS Y LIBRERÍAS -----------------------
datos <- read.csv(
"~/ESTADISTICA/dataset_geologico_limpio_80.csv",
header = TRUE,
sep = ",",
dec = "."
)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.3
##
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
# 2. SELECCIÓN DE VARIABLES ---------------------------
datos_seleccionados <- data.frame(
PHI_10 = as.numeric(datos$PHI_10),
CLAY_PCT = as.numeric(datos$CLAY_PCT)
)
cat("Datos originales seleccionados:", nrow(datos_seleccionados), "\n")
## Datos originales seleccionados: 27784
# 3. TRATAMIENTO DE DATOS -----------------------------
# 3.1 Eliminar datos vacíos
datos_tratados <- datos_seleccionados %>%
filter(!is.na(PHI_10), !is.na(CLAY_PCT))
cat("Después de eliminar vacíos:", nrow(datos_tratados), "\n")
## Después de eliminar vacíos: 26484
# 3.2 Filtrar valores físicamente posibles
# CLAY_PCT debe estar entre 0 y 100 %
datos_tratados <- datos_tratados %>%
filter(CLAY_PCT >= 0, CLAY_PCT <= 100)
cat("Después de filtrar CLAY_PCT entre 0 y 100:", nrow(datos_tratados), "\n")
## Después de filtrar CLAY_PCT entre 0 y 100: 26484
# 3.3 Eliminar valores atípicos en PHI_10 usando IQR
Q1_x <- quantile(datos_tratados$PHI_10, 0.25)
Q3_x <- quantile(datos_tratados$PHI_10, 0.75)
IQR_x <- Q3_x - Q1_x
LI_x <- Q1_x - 1.5 * IQR_x
LS_x <- Q3_x + 1.5 * IQR_x
datos_tratados <- datos_tratados %>%
filter(PHI_10 >= LI_x, PHI_10 <= LS_x)
cat("Después de eliminar atípicos en PHI_10:", nrow(datos_tratados), "\n")
## Después de eliminar atípicos en PHI_10: 24948
# 3.4 Eliminar valores atípicos en CLAY_PCT usando IQR
Q1_y <- quantile(datos_tratados$CLAY_PCT, 0.25)
Q3_y <- quantile(datos_tratados$CLAY_PCT, 0.75)
IQR_y <- Q3_y - Q1_y
LI_y <- Q1_y - 1.5 * IQR_y
LS_y <- Q3_y + 1.5 * IQR_y
datos_tratados <- datos_tratados %>%
filter(CLAY_PCT >= LI_y, CLAY_PCT <= LS_y)
cat("Después de eliminar atípicos en CLAY_PCT:", nrow(datos_tratados), "\n")
## Después de eliminar atípicos en CLAY_PCT: 24255
# 3.5 Agrupar PHI_10 por intervalos
# Se redondea PHI_10 cada 0.5 unidades para obtener un solo Y por cada X
datos_tratados <- datos_tratados %>%
mutate(
PHI_10_INTERVALO = round(PHI_10 / 0.5) * 0.5
)
# 3.6 Para cada X se calcula un único valor de Y
# Se usa el promedio de CLAY_PCT dentro de cada intervalo
TPV <- datos_tratados %>%
group_by(PHI_10_INTERVALO) %>%
summarise(
CLAY_PCT = mean(CLAY_PCT),
n_repeticiones = n()
) %>%
ungroup() %>%
filter(n_repeticiones >= 3)
colnames(TPV) <- c("PHI_10", "CLAY_PCT", "n_repeticiones")
cat("Datos finales tratados:", nrow(TPV), "\n")
## Datos finales tratados: 38
cat("\nTRATAMIENTO DE DATOS REALIZADO:\n")
##
## TRATAMIENTO DE DATOS REALIZADO:
cat("1. Se eliminaron registros vacíos en PHI_10 y CLAY_PCT.\n")
## 1. Se eliminaron registros vacíos en PHI_10 y CLAY_PCT.
cat("2. Se conservaron valores de CLAY_PCT entre 0 y 100%.\n")
## 2. Se conservaron valores de CLAY_PCT entre 0 y 100%.
cat("3. Se eliminaron valores atípicos usando el método IQR.\n")
## 3. Se eliminaron valores atípicos usando el método IQR.
cat("4. PHI_10 fue agrupado en intervalos de 0.5 unidades.\n")
## 4. PHI_10 fue agrupado en intervalos de 0.5 unidades.
cat("5. Para cada intervalo de PHI_10 se calculó el promedio de CLAY_PCT.\n")
## 5. Para cada intervalo de PHI_10 se calculó el promedio de CLAY_PCT.
cat("6. Se eliminaron intervalos con menos de 3 datos.\n")
## 6. Se eliminaron intervalos con menos de 3 datos.
cat("Resultado: cada valor de X tiene un único valor de Y.\n\n")
## Resultado: cada valor de X tiene un único valor de Y.
# 4. TABLA DE PARES DE VALORES ------------------------
View(TPV)
head(TPV, 20)
## # A tibble: 20 × 3
## PHI_10 CLAY_PCT n_repeticiones
## <dbl> <dbl> <int>
## 1 0 0.349 6636
## 2 0.5 1.62 3010
## 3 1 3.48 1565
## 4 1.5 5.47 1135
## 5 2 6.80 917
## 6 2.5 8.73 765
## 7 3 10.0 723
## 8 3.5 11.4 586
## 9 4 12.8 604
## 10 4.5 13.8 525
## 11 5 13.5 1253
## 12 5.5 16.9 441
## 13 6 18.5 442
## 14 6.5 19.2 371
## 15 7 20.2 403
## 16 7.5 21.2 348
## 17 8 20.6 327
## 18 8.5 23.0 315
## 19 9 23.7 316
## 20 9.5 24.6 299
# 5. DIAGRAMA DE DISPERSIÓN ---------------------------
x <- TPV$PHI_10
y <- TPV$CLAY_PCT
plot(
x, y,
pch = 16,
col = "blue",
main = "Gráfica Nº1: Diagrama de dispersión entre PHI_10 y Arcilla",
xlab = "PHI_10",
ylab = "Arcilla (%)"
)
# 6. CONJETURA DEL MODELO -----------------------------
cat("CONJETURA DEL MODELO\n")
## CONJETURA DEL MODELO
cat("Se plantea que existe una relación polinomial entre PHI_10 y CLAY_PCT.\n")
## Se plantea que existe una relación polinomial entre PHI_10 y CLAY_PCT.
cat("La relación observada no es lineal, por lo que se propone un modelo polinomial.\n\n")
## La relación observada no es lineal, por lo que se propone un modelo polinomial.
cat("Modelo propuesto:\n")
## Modelo propuesto:
cat("Y = a + bX + cX^2 + dX^3 + eX^4\n")
## Y = a + bX + cX^2 + dX^3 + eX^4
cat("Donde:\n")
## Donde:
cat("X = PHI_10\n")
## X = PHI_10
cat("Y = CLAY_PCT\n\n")
## Y = CLAY_PCT
# 7. CÁLCULO DE PARÁMETROS ----------------------------
xcuad <- x^2
xcub <- x^3
xcta <- x^4
regresion_polinomica <- lm(y ~ x + xcuad + xcub + xcta)
summary(regresion_polinomica)
##
## Call:
## lm(formula = y ~ x + xcuad + xcub + xcta)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5848 -0.5928 -0.2105 0.7100 2.0083
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.7847452 0.7413417 -1.059 0.297493
## x 5.0891182 0.5710346 8.912 2.67e-10 ***
## xcuad -0.5703944 0.1279891 -4.457 9.06e-05 ***
## xcub 0.0445880 0.0104618 4.262 0.000159 ***
## xcta -0.0012441 0.0002804 -4.436 9.61e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.062 on 33 degrees of freedom
## Multiple R-squared: 0.9919, Adjusted R-squared: 0.9909
## F-statistic: 1012 on 4 and 33 DF, p-value: < 2.2e-16
a <- regresion_polinomica$coefficients[1]
b <- regresion_polinomica$coefficients[2]
c <- regresion_polinomica$coefficients[3]
d <- regresion_polinomica$coefficients[4]
e <- regresion_polinomica$coefficients[5]
cat("\nECUACIÓN DEL MODELO\n")
##
## ECUACIÓN DEL MODELO
cat(
"Y =",
round(a, 4), "+",
round(b, 4), "X +",
round(c, 4), "X^2 +",
round(d, 4), "X^3 +",
round(e, 4), "X^4\n"
)
## Y = -0.7847 + 5.0891 X + -0.5704 X^2 + 0.0446 X^3 + -0.0012 X^4
# 8. SUPERPONER MODELO CON LA REALIDAD ----------------
plot(
x, y,
pch = 16,
col = "blue",
main = "Gráfica Nº2: Realidad vs Modelo Polinomial",
xlab = "PHI_10",
ylab = "Arcilla (%)"
)
curve(
a + b*x + c*x^2 + d*x^3 + e*x^4,
from = min(x),
to = max(x),
add = TRUE,
col = "red",
lwd = 2
)
legend(
"topleft",
legend = c("Datos tratados", "Modelo polinomial"),
col = c("blue", "red"),
pch = c(16, NA),
lty = c(NA, 1)
)
# 9. TEST DE PEARSON ----------------------------------
pearson <- cor.test(x, y, method = "pearson")
pearson
##
## Pearson's product-moment correlation
##
## data: x and y
## t = 29.191, df = 36, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9606576 0.9893913
## sample estimates:
## cor
## 0.9795232
r <- pearson$estimate
porcentaje_pearson <- r * 100
cat("\nCoeficiente de Pearson =", round(r, 4), "\n")
##
## Coeficiente de Pearson = 0.9795
cat("Porcentaje de Pearson =", round(porcentaje_pearson, 2), "%\n")
## Porcentaje de Pearson = 97.95 %
if(abs(porcentaje_pearson) >= 70){
cat("El test de Pearson APRUEBA porque supera el 70%.\n")
} else {
cat("El test de Pearson NO APRUEBA porque no supera el 70%.\n")
}
## El test de Pearson APRUEBA porque supera el 70%.
# 10. CÁLCULO DE ESTIMACIONES -------------------------
x0 <- 8
arcilla_esp <- a + b*x0 + c*x0^2 + d*x0^3 + e*x0^4
cat("\nESTIMACIÓN\n")
##
## ESTIMACIÓN
cat("Si PHI_10 =", x0, "\n")
## Si PHI_10 = 8
cat("El porcentaje estimado de arcilla es:", round(arcilla_esp, 4), "%\n")
## El porcentaje estimado de arcilla es: 21.1562 %
# 11. CONCLUSIONES ------------------------------------
r2 <- summary(regresion_polinomica)$r.squared
r2_ajustado <- summary(regresion_polinomica)$adj.r.squared
cat("\nCONCLUSIONES\n")
##
## CONCLUSIONES
cat("Variable independiente: PHI_10\n")
## Variable independiente: PHI_10
cat("Variable dependiente: CLAY_PCT\n\n")
## Variable dependiente: CLAY_PCT
cat("Antes de la regresión se realizó tratamiento de datos.\n")
## Antes de la regresión se realizó tratamiento de datos.
cat("Se eliminaron vacíos, valores fuera de rango y valores atípicos.\n")
## Se eliminaron vacíos, valores fuera de rango y valores atípicos.
cat("Luego PHI_10 fue agrupado en intervalos de 0.5 unidades.\n")
## Luego PHI_10 fue agrupado en intervalos de 0.5 unidades.
cat("Para cada intervalo se calculó el promedio de CLAY_PCT.\n")
## Para cada intervalo se calculó el promedio de CLAY_PCT.
cat("De esta forma, cada valor de X quedó asociado con un único valor de Y.\n\n")
## De esta forma, cada valor de X quedó asociado con un único valor de Y.
cat("Modelo obtenido:\n")
## Modelo obtenido:
cat(
"Y =",
round(a, 4), "+",
round(b, 4), "X +",
round(c, 4), "X^2 +",
round(d, 4), "X^3 +",
round(e, 4), "X^4\n\n"
)
## Y = -0.7847 + 5.0891 X + -0.5704 X^2 + 0.0446 X^3 + -0.0012 X^4
cat("Coeficiente de Pearson:", round(r, 4), "\n")
## Coeficiente de Pearson: 0.9795
cat("Porcentaje de Pearson:", round(porcentaje_pearson, 2), "%\n")
## Porcentaje de Pearson: 97.95 %
cat("R²:", round(r2, 4), "\n")
## R²: 0.9919
cat("R² ajustado:", round(r2_ajustado, 4), "\n\n")
## R² ajustado: 0.9909
if(abs(porcentaje_pearson) >= 70){
cat("El modelo cumple el criterio de correlación mayor al 70%.\n")
} else {
cat("El modelo no cumple el criterio de correlación mayor al 70%.\n")
}
## El modelo cumple el criterio de correlación mayor al 70%.
cat("El modelo permite estimar el porcentaje de arcilla a partir de PHI_10.\n")
## El modelo permite estimar el porcentaje de arcilla a partir de PHI_10.
cat("Las estimaciones son válidas únicamente dentro del rango observado de los datos tratados.\n")
## Las estimaciones son válidas únicamente dentro del rango observado de los datos tratados.