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 Conteo de Locales Comerciales (\(X_1\)) y el Factor de Sequías, Inundaciones y Temperatura (\(X_2\)), incrementan directamente el Valor Medio del Contaminante (\(Y\)), definiendo así el efecto en la calidad del agua.

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

# 2. 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)

# 3. 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.
# 4. Asignación directa según el formato de tu rúbrica
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

# ==============================================================================
# Paso 3: Depuración de Datos y Creación de la Tabla con gt()
# ==============================================================================
# Mantenemos los datos originales en su escala de medición real
TTP <- df

n <- nrow(TTP)
cat("Tamaño muestral del modelo = ", n)
## Tamaño muestral del modelo =  240
# Generamos la Tabla con los valores reales
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("*Valores crudos de Fósforo Total, Fosfato y Nitratos en cuerpos de agua de Europa*")
  ) %>%
  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
Valores crudos de Fósforo Total, Fosfato y Nitratos en cuerpos de agua de Europa
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

# ==============================================================================
# Paso 4: Gráfica de Dispersión Tridimensional Cruda
# ==============================================================================
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: Distribución tridimensional cruda de nutrientes hídricos"
)

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 corregir 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}\)) 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 de puntos influyentes
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 (225 observaciones listas)
TTP_limpio <- TTP_log[cooks_dist < umbral_cook, ]

# Confirmamos el tamaño en la consola de compilación
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: Distribución tridimensional aplicada logaritmo de nutrientes hídricos"
)

Al analizar visualmente la distribución de los datos transformados en la Gráfica N°2, se plantea la conjetura de que la relación de codependencia entre los nutrientes hídricos se puede modelar mediante un Plano de Regresión Lineal Múltiple con pendiente positiva de la forma

6. Parámetros

# ==============================================================================
# Paso 6 (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 exactos del modelo final
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

7. Comparación de la realidad con el modelo

# ==============================================================================
# Paso 7: Gráfica 3D con la Superficie Curva Ajustada (225 puntos)
# ==============================================================================

# 1. Graficamos la dispersión de los 225 puntos limpios
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), # Puntos con transparencia para ver la malla
  xlab = "Log10 Fósforo Total (x1)",
  ylab = "Log10 Fosfato (x2)",  
  zlab = "Log10 Nitratos (y)",
  main = "Gráfica N°2: Superficie cuadrática de regresión múltiple\najustada sobre datos depurados de nutrientes"
)

# 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)
}

8. Comparación de la realidad con el modelo

# ==============================================================================
# Paso 8: Test de Bondad de Ajuste Definitivo
# ==============================================================================

resumen_final <- summary(superficie_optimizada)

r2_final <- resumen_final$r.squared
r_final <- sqrt(r2_final)


cat("Coeficiente de Correlación (R): ", r_final)
## Coeficiente de Correlación (R):  0.7906285

9. 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 $X_1$ y $X_2$ 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.

10. 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

11. 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á influenciada en un 62.51% por la concentración de estos dos nutrientes en escala logarítmica, mientras que el 37.49% restante se debe a factores externos. El modelo no presenta restricciones físicas en la escala real de medición (\(Y > 0\)).