#=====================================================
# REGRESIÓN POLINOMIAL SIMPLE CON RELACIÓN 1 A 1
# 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 vacíos y valores físicamente imposibles
datos_limpios <- datos_seleccionados %>%
filter(!is.na(PHI_10), !is.na(CLAY_PCT)) %>%
filter(CLAY_PCT >= 0, CLAY_PCT <= 100) %>%
distinct(PHI_10, CLAY_PCT, .keep_all = TRUE)
cat("Datos después de eliminar vacíos, valores inválidos y duplicados exactos:",
nrow(datos_limpios), "\n")
## Datos después de eliminar vacíos, valores inválidos y duplicados exactos: 18930
# 3.2 Relación estricta 1 a 1
# Para cada X se obtiene un solo Y usando promedio.
# Luego, para cada Y se obtiene un solo X usando promedio.
# Se repite hasta que no existan valores repetidos ni en X ni en Y.
TPV <- datos_limpios
iteracion <- 1
repeat {
n_antes <- nrow(TPV)
# Un solo Y para cada X
TPV <- TPV %>%
group_by(PHI_10) %>%
summarise(
CLAY_PCT = mean(CLAY_PCT),
n_y_por_x = n(),
.groups = "drop"
)
# Redondear para evitar duplicados por decimales muy largos
TPV$CLAY_PCT <- round(TPV$CLAY_PCT, 4)
# Un solo X para cada Y
TPV <- TPV %>%
group_by(CLAY_PCT) %>%
summarise(
PHI_10 = mean(PHI_10),
n_x_por_y = n(),
.groups = "drop"
)
TPV$PHI_10 <- round(TPV$PHI_10, 4)
# Ordenar columnas
TPV <- TPV %>%
select(PHI_10, CLAY_PCT, n_x_por_y)
n_despues <- nrow(TPV)
cat("Iteración", iteracion,
"- observaciones:", n_antes, "→", n_despues, "\n")
repetidos_x <- sum(duplicated(TPV$PHI_10))
repetidos_y <- sum(duplicated(TPV$CLAY_PCT))
if(repetidos_x == 0 & repetidos_y == 0){
break
}
if(iteracion >= 10){
break
}
iteracion <- iteracion + 1
}
## Iteración 1 - observaciones: 18930 → 3857
## Iteración 2 - observaciones: 3857 → 3807
cat("\nDatos finales con relación 1 a 1:", nrow(TPV), "\n")
##
## Datos finales con relación 1 a 1: 3807
cat("Valores repetidos en X:", sum(duplicated(TPV$PHI_10)), "\n")
## Valores repetidos en X: 0
cat("Valores repetidos en Y:", sum(duplicated(TPV$CLAY_PCT)), "\n")
## Valores repetidos en Y: 0
cat("\nTRATAMIENTO DE DATOS REALIZADO:\n")
##
## TRATAMIENTO DE DATOS REALIZADO:
cat("1. Se eliminaron valores vacíos en PHI_10 y CLAY_PCT.\n")
## 1. Se eliminaron valores 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 registros duplicados exactos.\n")
## 3. Se eliminaron registros duplicados exactos.
cat("4. Para cada valor repetido de PHI_10 se calculó el promedio de CLAY_PCT.\n")
## 4. Para cada valor repetido de PHI_10 se calculó el promedio de CLAY_PCT.
cat("5. Para cada valor repetido de CLAY_PCT se calculó el promedio de PHI_10.\n")
## 5. Para cada valor repetido de CLAY_PCT se calculó el promedio de PHI_10.
cat("6. El proceso se repitió hasta obtener una relación estricta 1 a 1.\n\n")
## 6. El proceso se repitió hasta obtener una relación estricta 1 a 1.
# 4. TABLA DE PARES DE VALORES ------------------------
View(TPV)
head(TPV, 20)
## # A tibble: 20 × 3
## PHI_10 CLAY_PCT n_x_por_y
## <dbl> <dbl> <int>
## 1 7.87 0 1
## 2 0.01 0.0533 1
## 3 0.02 0.069 1
## 4 7.18 0.105 1
## 5 10.3 0.112 1
## 6 0.03 0.116 1
## 7 4.78 0.135 1
## 8 9.69 0.15 1
## 9 0.06 0.186 1
## 10 0.0754 0.197 1
## 11 0.328 0.204 1
## 12 0.07 0.256 1
## 13 0.08 0.279 1
## 14 8.22 0.3 1
## 15 0.11 0.306 1
## 16 3.34 0.329 1
## 17 0.12 0.352 1
## 18 3.34 0.352 1
## 19 0.14 0.369 1
## 20 3.99 0.410 1
# 5. GRÁFICA: 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("\nCONJETURA DEL MODELO\n")
##
## CONJETURA DEL MODELO
cat("Se plantea que existe una relación polinomial simple entre PHI_10 y CLAY_PCT.\n")
## Se plantea que existe una relación polinomial simple entre PHI_10 y CLAY_PCT.
cat("La relación observada no es completamente lineal.\n")
## La relación observada no es completamente lineal.
cat("Por ello se propone un modelo polinomial de cuarto grado.\n\n")
## Por ello se propone un modelo polinomial de cuarto grado.
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 DEL MODELO -----------------
regresion_polinomica <- lm(
y ~ poly(x, 4, raw = TRUE)
)
summary(regresion_polinomica)
##
## Call:
## lm(formula = y ~ poly(x, 4, raw = TRUE))
##
## Residuals:
## Min 1Q Median 3Q Max
## -49.442 -5.480 0.347 3.604 72.207
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.259e+01 7.180e-01 17.534 < 2e-16 ***
## poly(x, 4, raw = TRUE)1 8.791e-01 2.245e-01 3.916 9.15e-05 ***
## poly(x, 4, raw = TRUE)2 6.781e-02 2.039e-02 3.326 0.000891 ***
## poly(x, 4, raw = TRUE)3 -1.262e-03 6.645e-04 -1.898 0.057706 .
## poly(x, 4, raw = TRUE)4 2.620e-06 6.930e-06 0.378 0.705338
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.67 on 3802 degrees of freedom
## Multiple R-squared: 0.6684, Adjusted R-squared: 0.6681
## F-statistic: 1916 on 4 and 3802 DF, p-value: < 2.2e-16
coeficientes <- coef(regresion_polinomica)
a <- coeficientes[1]
b <- coeficientes[2]
c <- coeficientes[3]
d <- coeficientes[4]
e <- coeficientes[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 = 12.5903 + 0.8791 X + 0.0678 X^2 + -0.0013 X^3 + 0 X^4
# 8. SUPERPONER EL 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 = quantile(x, 0.01),
to = quantile(x, 0.99),
add = TRUE,
col = "red",
lwd = 2
)
legend(
"topleft",
legend = c("Datos tratados 1 a 1", "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 = 85.936, df = 3805, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8012937 0.8229122
## sample estimates:
## cor
## 0.8123819
r <- pearson$estimate
porcentaje_pearson <- r * 100
cat("\nCoeficiente de Pearson:", round(r, 4), "\n")
##
## Coeficiente de Pearson: 0.8124
cat("Porcentaje de Pearson:", round(porcentaje_pearson, 2), "%\n")
## Porcentaje de Pearson: 81.24 %
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 <- predict(
regresion_polinomica,
newdata = data.frame(x = x0)
)
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: 23.3276 %
# 11. CONCLUSIONES ------------------------------------
r2 <- summary(regresion_polinomica)$r.squared
r2_ajustado <- summary(regresion_polinomica)$adj.r.squared
cat("\nCONCLUSIONES\n")
##
## CONCLUSIONES
cat("Se seleccionó PHI_10 como variable independiente y CLAY_PCT como variable dependiente.\n")
## Se seleccionó PHI_10 como variable independiente y CLAY_PCT como variable dependiente.
cat("El tratamiento de datos permitió obtener una relación estricta 1 a 1.\n")
## El tratamiento de datos permitió obtener una relación estricta 1 a 1.
cat("Es decir, cada valor de PHI_10 quedó asociado con un único valor de CLAY_PCT,\n")
## Es decir, cada valor de PHI_10 quedó asociado con un único valor de CLAY_PCT,
cat("y cada valor de CLAY_PCT quedó asociado con un único valor de PHI_10.\n\n")
## y cada valor de CLAY_PCT quedó asociado con un único valor de PHI_10.
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 = 12.5903 + 0.8791 X + 0.0678 X^2 + -0.0013 X^3 + 0 X^4
cat("Coeficiente de Pearson:", round(r, 4), "\n")
## Coeficiente de Pearson: 0.8124
cat("Porcentaje de Pearson:", round(porcentaje_pearson, 2), "%\n")
## Porcentaje de Pearson: 81.24 %
cat("R²:", round(r2, 4), "\n")
## R²: 0.6684
cat("R² ajustado:", round(r2_ajustado, 4), "\n\n")
## R² ajustado: 0.6681
cat("El modelo explica aproximadamente", round(r2 * 100, 2),
"% de la variabilidad de CLAY_PCT.\n")
## El modelo explica aproximadamente 66.84 % de la variabilidad de CLAY_PCT.
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("Las estimaciones son válidas únicamente dentro del rango observado del dataset tratado.\n")
## Las estimaciones son válidas únicamente dentro del rango observado del dataset tratado.