#=====================================================
# REGRESIÓN POLINOMIAL SIMPLE
# X = MEDIAN
# Y = CLAY_PCT
#=====================================================

# 1. CARGAR 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
library(knitr)
library(gt)

# 2. SELECCIÓN DE VARIABLES ---------------------------

mediana_grano <- as.numeric(datos$MEDIAN)
arcilla <- as.numeric(datos$CLAY_PCT)

x <- mediana_grano      # Variable independiente
y <- arcilla            # Variable dependiente

# 3. TABLA DE PARES DE VALORES ------------------------

TPV <- data.frame(x, y)
TPV <- na.omit(TPV)

colnames(TPV) <- c("MEDIAN", "CLAY_PCT")

View(TPV)

kable(
  head(TPV, 20),
  caption = "Tabla de pares de valores: MEDIAN y CLAY_PCT"
)
Tabla de pares de valores: MEDIAN y CLAY_PCT
MEDIAN CLAY_PCT
1 -0.290000 0.00000
2 0.430000 0.00000
3 1.750000 0.00000
4 1.700000 0.00000
5 0.970000 0.00000
6 3.270000 12.70000
7 1.680000 19.60000
8 -1.740000 0.00000
9 7.644568 13.73566
10 1.170000 0.00000
11 -0.660000 0.00000
13 3.190000 16.20000
14 4.750000 19.70000
15 6.360000 29.70000
16 2.740000 0.00000
17 1.070000 0.00000
18 2.140000 0.00000
19 1.800000 0.00000
20 1.770000 0.00000
21 1.220000 0.00000 
# 4. GRÁFICA: DIAGRAMA DE DISPERSIÓN ------------------

plot(
  TPV$MEDIAN,
  TPV$CLAY_PCT,
  pch = 16,
  col = "blue",
  main = "Gráfica Nº1: Diagrama de dispersión entre MEDIAN y Arcilla",
  xlab = "MEDIAN",
  ylab = "Arcilla (%)"
)

# 5. CONJETURA DEL MODELO -----------------------------

cat("CONJETURA DEL MODELO\n")
## CONJETURA DEL MODELO
cat("Se plantea que existe una relación polinomial entre MEDIAN y CLAY_PCT.\n")
## Se plantea que existe una relación polinomial entre MEDIAN y CLAY_PCT.
cat("A medida que cambia el valor de MEDIAN, también cambia el porcentaje de arcilla.\n")
## A medida que cambia el valor de MEDIAN, también cambia el porcentaje de arcilla.
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 = MEDIAN\n")
## X = MEDIAN
cat("Y = CLAY_PCT\n\n")
## Y = CLAY_PCT
# 6. CÁLCULO DE PARÁMETROS DEL MODELO -----------------

x <- TPV$MEDIAN
y <- TPV$CLAY_PCT

xcuad <- x^2
xcub <- x^3
xcta <- x^4

regresion_polinomica <- lm(y ~ x + xcuad + xcub + xcta)

regresion_polinomica
## 
## Call:
## lm(formula = y ~ x + xcuad + xcub + xcta)
## 
## Coefficients:
## (Intercept)            x        xcuad         xcub         xcta  
##    2.436694    -0.044828     0.235973     0.066409    -0.001011
summary(regresion_polinomica)
## 
## Call:
## lm(formula = y ~ x + xcuad + xcub + xcta)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -154.751   -3.691   -2.252    1.892   89.088 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.4366945  0.1304421  18.680  < 2e-16 ***
## x           -0.0448279  0.0731116  -0.613  0.53979    
## xcuad        0.2359729  0.0131723  17.914  < 2e-16 ***
## xcub         0.0664087  0.0036743  18.074  < 2e-16 ***
## xcta        -0.0010105  0.0003175  -3.183  0.00146 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.075 on 26770 degrees of freedom
## Multiple R-squared:  0.7549, Adjusted R-squared:  0.7548 
## F-statistic: 2.061e+04 on 4 and 26770 DF,  p-value: < 2.2e-16
beta0 <- regresion_polinomica$coefficients[1]
beta1 <- regresion_polinomica$coefficients[2]
beta2 <- regresion_polinomica$coefficients[3]
beta3 <- regresion_polinomica$coefficients[4]
beta4 <- regresion_polinomica$coefficients[5]

a <- beta0
b <- beta1
c <- beta2
d <- beta3
e <- beta4

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 = 2.4367 + -0.0448 X + 0.236 X^2 + 0.0664 X^3 + -0.001 X^4
# 7. SUPERPONER EL MODELO CON LA REALIDAD --------------

plot(
  x,
  y,
  pch = 16,
  col = "blue",
  main = "Gráfica Nº2: Realidad vs Modelo Polinomial",
  xlab = "MEDIAN",
  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(
  "topright",
  legend = c("Datos reales", "Modelo polinomial"),
  col = c("blue", "red"),
  pch = c(16, NA),
  lty = c(NA, 1)
)

# Mostrar ecuación en una gráfica aparte

plot(
  1,
  type = "n",
  axes = FALSE,
  xlab = "",
  ylab = ""
)

eq <- paste0(
  "Ecuación polinómica de cuarto grado\n\n",
  "Y = a + bX + cX^2 + dX^3 + eX^4\n\n",
  "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",
  "Donde:\n",
  "X = MEDIAN\n",
  "Y = Arcilla (%)"
)

text(
  1,
  1,
  labels = eq,
  cex = 1.1,
  col = "blue",
  font = 2
)

# 8. TEST DE PEARSON ----------------------------------

pearson <- cor.test(x, y, method = "pearson")

pearson
## 
##  Pearson's product-moment correlation
## 
## data:  x and y
## t = 204.8, df = 26773, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7765614 0.7858958
## sample estimates:
##       cor 
## 0.7812723
r <- pearson$estimate
porcentaje_pearson <- r * 100

cat("\nCoeficiente de Pearson =", r, "\n")
## 
## Coeficiente de Pearson = 0.7812723
cat("Porcentaje de Pearson =", porcentaje_pearson, "%\n")
## Porcentaje de Pearson = 78.12723 %
# 9. CÁLCULO DE ESTIMACIONES --------------------------

# Ejemplo:
# ¿Cuál sería el porcentaje de arcilla
# si MEDIAN = 6?

x0 <- 6

arcilla_esp <- a +
  b*x0 +
  c*x0^2 +
  d*x0^3 +
  e*x0^4

arcilla_esp
## (Intercept) 
##    23.69741
plot(
  1,
  type = "n",
  axes = FALSE,
  xlab = "",
  ylab = ""
)

text(
  1,
  1,
  labels = paste0(
    "¿Cuál sería el porcentaje de arcilla\n",
    "si MEDIAN = 6?\n\n",
    "R = ",
    round(arcilla_esp, 4),
    " %"
  ),
  cex = 1.4,
  col = "blue",
  font = 2
)

# 10. CONCLUSIONES ------------------------------------

r2 <- summary(regresion_polinomica)$r.squared

cat("\nCONCLUSIONES\n")
## 
## CONCLUSIONES
cat("Entre MEDIAN y CLAY_PCT existe una relación polinomial simple.\n\n")
## Entre MEDIAN y CLAY_PCT existe una relación polinomial simple.
cat("El modelo obtenido fue:\n")
## El modelo obtenido fue:
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 = 2.4367 + -0.0448 X + 0.236 X^2 + 0.0664 X^3 + -0.001 X^4
cat("El coeficiente de Pearson fue:", round(r, 4), "\n")
## El coeficiente de Pearson fue: 0.7813
cat("El porcentaje de Pearson fue:", round(porcentaje_pearson, 2), "%\n")
## El porcentaje de Pearson fue: 78.13 %
cat("El coeficiente de determinación R² fue:", round(r2, 4), "\n\n")
## El coeficiente de determinación R² fue: 0.7549
if(abs(porcentaje_pearson) >= 70){
  cat("El modelo APRUEBA el criterio de correlación mayor al 70%.\n")
} else {
  cat("El modelo NO aprueba el criterio de correlación mayor al 70%.\n")
}
## El modelo APRUEBA el criterio de correlación mayor al 70%.
cat("El modelo permite estimar el porcentaje de arcilla a partir de MEDIAN.\n")
## El modelo permite estimar el porcentaje de arcilla a partir de MEDIAN.
cat("Sin embargo, las estimaciones solo son confiables dentro del rango observado del dataset.\n")
## Sin embargo, las estimaciones solo son confiables dentro del rango observado del dataset.