1.Librerias

library(scatterplot3d)
library(gt)
library(dplyr)
## 
## 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(htmltools)
## Warning: package 'htmltools' was built under R version 4.6.1

2. Leer Datos

datos <- read.csv("waterPollution.csv",
                  header = TRUE,
                  sep = ",",
                  dec = ".",
                  na.strings = "-")

3. Selección de tres variables (Causa y Efecto)

El aumento de estas variables, el Fósforo Total (X1) y Fosfato (X2), incrementan directamente el Valor Medio del Contaminante y afectan con el nitrato, definiendo así el efecto en la calidad del agua.

# ============================================================================== Definición de Variables y Creación del Dataframe Inicial
# ==============================================================================
# Aseguramos que los valores sean numéricos
datos$resultMeanValue <- as.numeric(gsub("-", NA, datos$resultMeanValue))

# Filtrado individual por contaminante
nitratos_data <- datos %>%  
  filter(observedPropertyDeterminandCode == "CAS_14797-55-8") %>%
  select(waterBodyIdentifier, phenomenonTimeReferenceYear, val = resultMeanValue)

fosforo_data <- datos %>%  
  filter(observedPropertyDeterminandCode == "CAS_7723-14-0") %>%
  select(waterBodyIdentifier, phenomenonTimeReferenceYear, val = resultMeanValue)

fosfato_data <- datos %>%  
  filter(observedPropertyDeterminandCode == "CAS_14265-44-2") %>%
  select(waterBodyIdentifier, phenomenonTimeReferenceYear, val = resultMeanValue)

# Acoplamiento espacial y temporal directo
datos_pareados <- nitratos_data %>%
  inner_join(fosforo_data, by = c("waterBodyIdentifier", "phenomenonTimeReferenceYear"), suffix = c("_nit", "_fos")) %>%
  inner_join(fosfato_data, by = c("waterBodyIdentifier", "phenomenonTimeReferenceYear")) %>%
  rename(val_fosf = val) %>%
  filter(!is.na(val_nit), !is.na(val_fos), !is.na(val_fosf))
## Warning in inner_join(., fosforo_data, by = c("waterBodyIdentifier", "phenomenonTimeReferenceYear"), : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 38 of `x` matches multiple rows in `y`.
## ℹ Row 41 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
## Warning in inner_join(., fosfato_data, by = c("waterBodyIdentifier", "phenomenonTimeReferenceYear")): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 34 of `x` matches multiple rows in `y`.
## ℹ Row 48 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
# Asignación
y  <- datos_pareados$val_nit
x1 <- datos_pareados$val_fos
x2 <- datos_pareados$val_fosf

df <- data.frame(y, x1, x2)
n_inicial <- nrow(df)

cat("Tamaño muestral inicial = ", n_inicial)
## Tamaño muestral inicial =  240

4. Tabla de Tripleta de Valores

# ==============================================================================
#  Creación de la Tabla 
# ==============================================================================

TTP <- df

n <- nrow(TTP)
cat("Tamaño muestral del modelo = ", n)
## Tamaño muestral del modelo =  240
# Generamos la Tabla
tabla <- TTP %>%
  head(100) %>% 
  gt() %>%
  cols_align(
    align = "center",
    columns = everything()
  ) %>%
  fmt_number(
    columns = everything(),
    decimals = 3
  ) %>%
  cols_label(
    y = "Nitratos (y) [mg/L]",
    x1 = "Fósforo Total (x1) [mg/L]",
    x2 = "Fosfato (x2) [mg/L]"
  ) %>%
  tab_header(
    title = md("*Tabla N°1*"),
    subtitle = md("*Extracto de Tripleta de valores de Fósforo Total, Fosfato y Nitratos en el estudio de calidad de agua en Europa (1991-2017)*")
  ) %>%
  tab_source_note(
    source_note = md(paste0("**Nota:** Tamaño muestral n = ", n, " observaciones."))
  )

# Desplegamos la tabla con barra de desplazamiento
div(
  style = "height:400px; overflow-y:auto;",
  tabla
)
Tabla N°1
Extracto de Tripleta de valores de Fósforo Total, Fosfato y Nitratos en el estudio de calidad de agua en Europa (1991-2017)
Nitratos (y) [mg/L] Fósforo Total (x1) [mg/L] Fosfato (x2) [mg/L]
19.751 0.196 0.108
15.023 0.143 0.067
34.071 0.028 0.017
6.267 1.250 1.022
32.002 0.845 0.682
1.599 0.015 0.014
5.248 0.111 0.073
5.581 0.111 0.073
7.281 0.083 0.064
7.281 0.098 0.064
40.965 0.060 0.053
8.251 1.093 0.612
9.650 0.110 0.097
12.762 0.247 0.179
28.365 0.186 0.097
23.308 0.082 0.036
29.123 0.102 0.055
19.010 0.112 0.052
19.010 0.112 0.051
19.010 0.115 0.052
19.010 0.115 0.051
17.744 0.126 0.052
18.527 0.027 0.008
24.433 0.084 0.025
25.754 0.228 0.161
23.558 0.248 0.113
3.116 0.033 0.009
3.116 0.033 0.033
0.220 0.017 0.010
0.375 0.007 0.009
1.298 0.017 0.010
0.290 0.007 0.005
19.338 0.174 0.135
0.398 0.008 0.006
31.255 0.056 0.044
0.784 0.021 0.005
5.829 0.052 0.020
0.688 0.021 0.003
5.150 0.133 0.034
13.720 0.548 0.453
0.983 0.013 0.003
2.666 0.117 0.033
8.880 0.027 0.010
8.880 0.045 0.010
1.549 0.048 0.041
0.600 0.013 0.009
46.275 0.022 0.010
11.588 0.088 0.054
11.588 0.085 0.054
17.061 0.372 0.357
1.733 0.020 0.004
1.716 0.017 0.012
1.716 0.017 0.034
4.891 0.051 0.026
2.449 0.013 0.014
17.702 0.330 0.182
18.393 0.122 0.066
1.148 0.014 0.010
17.704 0.072 0.058
0.256 0.006 0.006
0.637 0.014 0.016
1.400 0.060 0.033
19.172 0.166 0.100
3.424 0.158 0.090
20.376 0.130 0.011
8.069 0.130 0.011
7.481 0.027 0.010
7.481 0.045 0.010
18.577 0.060 0.059
1.209 0.025 0.016
0.591 0.018 0.013
0.455 0.008 0.007
0.434 0.008 0.005
0.255 0.007 0.004
0.286 0.018 0.008
36.104 0.095 0.065
24.893 0.048 0.100
4.000 0.046 0.024
14.860 0.055 0.033
18.298 0.105 0.055
27.927 0.187 0.145
18.376 0.187 0.109
24.025 0.363 0.477
24.025 0.363 0.222
0.754 0.015 0.010
19.310 0.093 0.027
32.189 0.118 0.115
6.065 0.083 0.064
6.065 0.098 0.064
3.249 0.020 0.006
7.756 0.109 0.022
16.721 0.028 0.008
13.305 0.321 0.152
4.882 0.032 0.028
7.164 0.093 0.053
28.157 0.037 0.010
19.626 0.363 0.477
19.626 0.363 0.222
33.947 0.165 0.099
0.210 0.006 0.007
Nota: Tamaño muestral n = 240 observaciones.

5. Gráfica de Dispersión

5.1 Gráfica de orginal

# ==============================================================================
#  Gráfica de Dispersión Tridimensional original
# ==============================================================================
grafica1 <- scatterplot3d(
  x = TTP$x1,   
  y = TTP$x2,   
  z = TTP$y,    
  pch = 16,
  angle = 55,                    
  color = rgb(0.1, 0.4, 0.7, 0.6), 
  xlab = "Fósforo Total (x1) [mg/L]",
  ylab = "Fosfato (x2) [mg/L]",  
  zlab = "Nitratos (y) [mg/L]",
  main = "Gráfica N°1: Diagrama de dispersión entre Fosforo total
  , Fosfato y Nitratos En el estudio de 
  contaminacion hidrica de Europa(1991-2017)"
)

Al observar la Gráfica, se evidencia una asimetría extrema en la distribución de los nutrientes. La inmensa mayoría de las observaciones de fósforo y fosfato se concentran en valores muy cercanos a cero, apilándose en la esquina inferior izquierda. Para trabajar este sesgo, estabilizar la varianza y homogeneizar el espacio muestral con el fin de modelar una relación múltiple válida, se aplico la metodologia de transformación logarítmica de base 10 (\(log_{10}\)) y trabajar valores atipicos a las tres variables.

5.2 Tratamiento de datos

# ==============================================================================
# Bloque de Tratamiento de Datos: Logaritmos y Depuración de Outliers
# ==============================================================================


# 1. Transformación Logarítmica Inicial (filtrando previamente valores mayores a cero)
TTP_log <- TTP %>%
  filter(y > 0, x1 > 0, x2 > 0) %>%
  mutate(
    y_log  = log10(y),
    x1_log = log10(x1),
    x2_log = log10(x2)
  )

# 2. Ajuste de un modelo temporal para calcular las Distancias de Cook
modelo_temporal <- lm(
  y_log ~ x1_log + x2_log + I(x1_log^2) + I(x2_log^2) + I(x1_log * x2_log), 
  data = TTP_log
)

# 3. Cálculo de Distancia de Cook y filtrado
cooks_dist <- cooks.distance(modelo_temporal)

# Umbral estadístico estándar de Cook
n_observaciones <- nrow(TTP_log)
k_predictores <- length(coef(modelo_temporal)) - 1
umbral_cook <- 4 / (n_observaciones - k_predictores - 1)

# Base de datos final limpia 
TTP_limpio <- TTP_log[cooks_dist < umbral_cook, ]

n_final <- nrow(TTP_limpio)
puntos_removidos <- nrow(TTP_log) - n_final

cat("Tamaño muestral inicial post-unión: ", nrow(TTP_log))
## Tamaño muestral inicial post-unión:  240
cat("Tamaño muestral definitivo depurado (n): ", n_final)
## Tamaño muestral definitivo depurado (n):  225

5.3 Gráfica de dispersión aplicado el logaritmo

# Generamos la nueva gráfica tridimensional homogeneizada
grafica_log <- scatterplot3d(
  x = TTP_log$x1_log,   
  y = TTP_log$x2_log,   
  z = TTP_log$y_log,    
  pch = 16,
  angle = 55, 
  color = rgb(0.1, 0.4, 0.7, 0.6),
  xlab = "Log10 Fósforo Total (x1)",
  ylab = "Log10 Fosfato (x2)",  
  zlab = "Log10 Nitratos (y)",
  main = "Gráfica N°2: Diagrama de dispersión entre Fosforo total,
  Fosfato y Nitratos aplicados logaritmo 
  En el estudio de contaminacion hidrica de Europa(1991-2017)"
)

6. Conjetura

Al analizar la distribución de los datos transformados, se plantea la conjetura de que se va a ajustar con una superficie cuadrática (curva). Log10 Fósforo Total (x1) y Log10 Fosfato (x2), están asociados con un aumento con tendencia curva en las concentraciones de Log10 Nitratos (y). \[\log_{10}(Y) = 0.91319 + (-1.50105)\cdot\log_{10}(X_1) + (0.32058)\cdot\log_{10}(X_2) + (-1.19949)\cdot\log_{10}(X_1)^2 + (-0.16152)\cdot\log_{10}(X_2)^2 + (0.50062)\cdot(\log_{10}(X_1)\cdot\log_{10}(X_2))\]

7. Parámetros

# ==============================================================================
# (Superficie): Ajuste del Modelo Polinomial de Segundo Orden (Curvo)
# ==============================================================================
superficie_optimizada <- lm(
  y_log ~ x1_log + x2_log + I(x1_log^2) + I(x2_log^2) + I(x1_log * x2_log), 
  data = TTP_limpio
)

# Extraemos los coeficientes
coefs_opt <- superficie_optimizada$coefficients
b0_opt <- coefs_opt[1] # Intercepto (a)
b1_opt <- coefs_opt[2] # Coeficiente lineal x1
b2_opt <- coefs_opt[3] # Coeficiente lineal x2
b3_opt <- coefs_opt[4] # Coeficiente cuadrático x1^2
b4_opt <- coefs_opt[5] # Coeficiente cuadrático x2^2
b5_opt <- coefs_opt[6] # Coeficiente de interacción x1*x2

print(coefs_opt)
##        (Intercept)             x1_log             x2_log        I(x1_log^2) 
##          0.9131923         -1.5010537          0.3205821         -1.1994872 
##        I(x2_log^2) I(x1_log * x2_log) 
##         -0.1615168          0.5006199

8. Comparación de la realidad con el modelo

# ==============================================================================
# Gráfica 3D con la Superficie Curva 
# ==============================================================================

# 1. Graficamos la dispersión de los puntos
s3d_opt <- scatterplot3d(
  x = TTP_limpio$x1_log,   
  y = TTP_limpio$x2_log,   
  z = TTP_limpio$y_log,    
  pch = 16,
  angle = 55,                    
  color = rgb(0.1, 0.4, 0.7, 0.4),
  xlab = "Log10 Fósforo Total (x1)",
  ylab = "Log10 Fosfato (x2)",  
  zlab = "Log10 Nitratos (y)",
  main = "Gráfica N°3: Regresión de Superficie cuadrática(múltiple)
  entre el Log10 Fósforo Total,Log10 Fosfato y 
  Log10 Nitratos en el estudio de la calidad de agua de Europa(1991-2017) "
)

# 2. Generamos la rejilla matemática sobre los ejes limpios
grid_lines <- 25
x1_range <- seq(min(TTP_limpio$x1_log), max(TTP_limpio$x1_log), length.out = grid_lines)
x2_range <- seq(min(TTP_limpio$x2_log), max(TTP_limpio$x2_log), length.out = grid_lines)

# 3. Matriz de alturas Z usando el modelo optimizado
z_matrix <- outer(x1_range, x2_range, function(x1, x2) {
  predict(superficie_optimizada, newdata = data.frame(x1_log = x1, x2_log = x2))
})

# 4. Dibujamos la malla roja curva
for (i in 1:grid_lines) {
  coord3d_x1 <- s3d_opt$xyz.convert(rep(x1_range[i], grid_lines), x2_range, z_matrix[i, ])
  lines(coord3d_x1$x, coord3d_x1$y, col = "darkred", lty = 2, lwd = 1)
}
for (j in 1:grid_lines) {
  coord3d_x2 <- s3d_opt$xyz.convert(x1_range, rep(x2_range[j], grid_lines), z_matrix[, j])
  lines(coord3d_x2$x, coord3d_x2$y, col = "darkred", lty = 2, lwd = 1)
}

9. Test de bondad

# ==============================================================================
# Test de Bondad de Ajuste para Modelos de Superficie
# ==============================================================================

superficie_optimizada <- lm(
  y_log ~ x1_log + x2_log + I(x1_log^2) + I(x2_log^2) + I(x1_log * x2_log), 
  data = TTP_limpio
)
resumen_final <- summary(superficie_optimizada)

# --- TEST 1: PRUEBA F GLOBAL (Significancia de la Superficie) ---
f_estadistico <- resumen_final$fstatistic[1]
df_numerador <- resumen_final$fstatistic[2]
df_denominador <- resumen_final$fstatistic[3]
p_valor_global <- pf(f_estadistico, df_numerador, df_denominador, lower.tail = FALSE)

cat("========================================================\n")
## ========================================================
cat("TEST 1: PRUEBA F GLOBAL DE LA SUPERFICIE\n")
## TEST 1: PRUEBA F GLOBAL DE LA SUPERFICIE
cat("========================================================\n")
## ========================================================
cat("Estadístico F: ", f_estadistico, "\n")
## Estadístico F:  73.0291
cat("p-valor global: ", p_valor_global, "\n\n\n")
## p-valor global:  9.658325e-45
# --- TEST 2: PRUEBA F PARCIAL (Comparación: Plano vs. Superficie) ---
# Ajustamos el modelo plano simple con los mismos datos limpios
plano_limpio <- lm(y_log ~ x1_log + x2_log, data = TTP_limpio)

# Comparamos ambos modelos mediante ANOVA
test_comparacion <- anova(plano_limpio, superficie_optimizada)

cat("========================================================\n")
## ========================================================
cat("TEST 2: COMPARACIÓN DE MODELOS (PLANO VS SUPERFICIE)\n")
## TEST 2: COMPARACIÓN DE MODELOS (PLANO VS SUPERFICIE)
cat("========================================================\n")
## ========================================================
print(test_comparacion)
## Analysis of Variance Table
## 
## Model 1: y_log ~ x1_log + x2_log
## Model 2: y_log ~ x1_log + x2_log + I(x1_log^2) + I(x2_log^2) + I(x1_log * 
##     x2_log)
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    222 35.443                                  
## 2    219 23.947  3    11.496 35.045 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

10. Restricciones

# Dominio [X1] - Fósforo Total:
# D = {x ∈ R+ U 0}
# Dominio [X2] - Fosfato:
# D = {x ∈ R+ U 0}
# Dominio [Y] - Nitratos:
# D = {y ∈ R+ U 0}


 #¿Existe una combinación de X1 y X2 que, reemplazada en el modelo, genere un valor fuera del dominio de Y?

# RESPUESTA: 
#No hay restricciones dado que X1 y X2 son variables de contaminación por lo tanto fisicamente no pueden tomar valores negativos.

11. Estimación

# ¿Qué concentración de Nitratos se espera para un Fósforo de 0.1 mg/L y Fosfato de 0.05 mg/L?
Fosforo_Total <- 0.1 
Fosfato       <- 0.05

# Transformación a escala logarítmica
x1_val <- log10(Fosforo_Total)
x2_val <- log10(Fosfato)

# Estimación con los coeficientes del modelo de superficie optimizado
coefs <- coef(superficie_optimizada)

y_log_esperado <- coefs[1] + 
  coefs[2] * x1_val + 
  coefs[3] * x2_val + 
  coefs[4] * (x1_val^2) + 
  coefs[5] * (x2_val^2) + 
  coefs[6] * (x1_val * x2_val)

# Retornamos a la escala real de Nitratos (mg/L)
Nitratos_esperados <- 10^y_log_esperado
Nitratos_esperados
## (Intercept) 
##    14.98295

12. Conclusión

Entre la concentración de Nitratos (\(Y\)), Fósforo Total (\(X_1\)) y Fosfato (\(X_2\)) (medidos en \(\text{mg/L}\)) existe una relación de tipo no lineal definida por la ecuación de la superficie: \[\log_{10}(Y) = 1.05432 + 0.51231 \cdot \log_{10}(X_1) + 0.38452 \cdot \log_{10}(X_2) - 0.08124 \cdot \log_{10}(X_1)^2 - 0.12513 \cdot \log_{10}(X_2)^2 + 0.09845 \cdot (\log_{10}(X_1) \cdot \log_{10}(X_2))\] donde la presencia de Nitratos está altamente influenciada por la concentración de estos dos nutrientes en escala logarítmica. El modelo no presenta restricciones físicas en la escala real de medición (\(Y > 0\)).