Introducción
La regresión lineal múltiple (RLM) es una técnica estadística básica empleada por los hidrólogos para transferir información de las características del escurrimiento, calculadas en cuencas con hidrometría hacia sitios o cuencas donde tal información es necesaria y no existen aforos. Por ejemplo, la RLM se emplea para completar series de volumen escurrido anual (Y) en sitios con pocos datos, con base en varios registros amplios y cercanos (Xp); también se ha utilizado para encontrar ecuaciones empíricas que relacionan la creciente media anual (Qma) o de cierto periodo de retorno (QTr ) con las características físicas de las cuencas con aforos, por lo común el tamaño de la cuenca y ciertas propiedades del cauce, o tormentas de la zona (Pandey & Nguyen, 1999; Griffis & Stedinger, 2007; Salas et al, 2008; Wilks, 2011).
Desde el inicio del siglo XXI se han comenzado a notar los efectos negativos del cambio climático a través de observar eventos extremos, como tormentas, crecientes y sequías más severos y más frecuentes. Ante tal evidencia se han sugerido varias acciones (Campos-Aranda, 2015), una de ellas, quizá la más importante, sea aumentar los sitios de registro de lluvias, gastos y niveles en lagos, embalses y acuíferos. Sin embargo, la realidad demuestra que el número de estaciones pluviográficas, pluviométricas e hidrométricas ha disminuido notablemente, hasta llegar a una situación crítica, como lo han destacado Lafragua-Contreras, González-Rojas y Solís-Alvarado (2006). En este escenario, es prioritario hacer un uso eficiente de la información hidrológica disponible y la RLM es la técnica estadística que lo permite, al trabajar los datos regionalmente tanto en el transporte de datos como en la estimación de ecuaciones empíricas.
Los parámetros de ajuste de un modelo definido por una RLM se estiman a través del método o técnica de mínimos cuadrados ordinarios (MCO), la cual acepta o asume que las varianzas de la variable dependiente (Yi) son las mismas para toda i, lo cual se conoce como condición de homocedasticidad. Lo anterior implica que todas las observaciones de Y sean "igualmente confiables". En cualquier aplicación hidrológica regional de la RLM, tal condición seguramente será violada debido a que la confiabilidad en la estimación de las características del escurrimiento depende de la amplitud de su registro, o bien de las condiciones de medición en las estaciones de aforos. Cuando la hipótesis de homocedasticidad no es satisfecha, los parámetros estimados no tienen varianza mínima y todas las estimaciones asociadas con la RLM no son exactas (Tasker, 1980; Stedinger & Tasker, 1985; Tasker & Stedinger, 1989; Pandey & Nguyen, 1999; Kottegoda & Rosso, 2008).
Como las varianzas de la variable dependiente σ2(Yi) no son iguales, se puede aplicar una función de ponderado (wi) que corrija tal variación, según se expone en la teoría del método de mínimos cuadrados ponderados (MCP). Tal función es óptima cuando wi = 1/σ2(Yi).
El objetivo de este trabajo consiste en exponer con detalle dos procedimientos para encontrar y aplicar la función de ponderado óptima. El primero utiliza el algoritmo desarrollado por Tasker (1980), que toma en cuenta la teoría de los residuales y los resultados del método de MCO. El segundo se basa en una propiedad de los valores cercanos del regresor (X1), que son considerados puntos de repetición. Ambos procedimientos se aplican en la obtención de dos ecuaciones de regresión potencial, que permiten estimar el gasto máximo medio anual (Qma) de la Región Hidrológica 10 (Sinaloa, México). Como los resultados del método de MCP mejoran el ajuste de MCO, según indicadores evaluados en el dominio real, se recomienda su aplicación sistemática al estimar ecuaciones empíricas por RLM.
Resumen de la teoría operativa
Regresión lineal múltiple (RLM)
Algunas veces se puede establecer una relación de tipo lineal entre la variable dependiente (Y) y varias (p) independientes X1, X2..., Xp o regresores, la cual es la generalización o extensión natural de la regresión lineal simple. Su expresión es (Ryan, 1998):
Entonces, los principios que rigen la regresión lineal se aplican a la RLM, por ejemplo, que tanto Y como las Xp estén normalmente distribuidas, y que los errores e sean independientes y tengan distribución normal de media cero y misma varianza (σ2) para cada X. La solución de mínimos cuadrados de los residuos de forma matricial para el caso general expuesto y con n observaciones o datos de Y, y de los regresores es la siguiente (Ryan, 1998):
siendo:
El planteamiento de esta solución implica que la sumatoria de 1 a n de los residuos al cuadrado debe ser minimizada, es decir que:
Entonces, diferenciando el lado derecho de la ecuación anterior con respecto a β0, β1, β2,.. βp por separado, se originan las ecuaciones llamadas normales, función de los parámetros desconocidos. En notación matricial, estas ecuaciones son:
cuya solución es:
en la cual X' es la matriz transpuesta de X y (X'-X)-1 indica la matriz inversa de X'-X.
Solución de mínimos cuadrados ponderados
Como ya se indicó, las suposiciones que por lo común se establecen en relación con RLM (ecuación (2)) son que E(ε) = 0 y que Var(e) = σ2-I, siendo I la matriz unitaria o identidad. Con frecuencia, tales premisas son irrazonables, pues se tiene que Var(ε) = a2-V, siendo V una matriz conocida de n x n. Si V es diagonal, con elementos diagonales distintos, las observaciones Y no están correlacionadas, pero tienen varianzas desiguales; en cambio, si existen algunos elementos fuera de la diagonal principal de V, las observaciones están correlacionadas y la solución de mínimos cuadrados de los residuos es (Montgomery, Peck, & Vining, 2002):
En la ecuación anterior,
en la cual
Aplicando MCO a los datos transformados, se obtiene el estimador de mínimos cuadrados ponderados, que será:
Para usar la técnica de mínimos cuadrados ponderados se deben conocer los pesos wi. Con frecuencia se puede recurrir a la experiencia o conocimiento previo, a la información de un modelo teórico, o bien el análisis de los residuos puede indicar que la varianza de los errores puede ser una función de uno de los regresores; por ejemplo, si Var(ε) = σ2.XLi , entonces wi = 1/ X1i r Incluso, en aplicaciones prácticas, se pueden suponer los pesos y hacer iteraciones para mejorar la regresión y/o minimizar algunos residuos (Montgomery et al., 2002).
Indicadores de calidad del ajuste
Cuando se obtienen ecuaciones de RLM con mínimos cuadrados ponderados es necesario compararlas para escoger la de mejor ajuste, lo cual implica evaluar indicadores basados en los residuos. El más común de tales indicadores es el coeficiente de determinación (R2), que indica la proporción de la varianza de la variable dependiente que es explicada por la ecuación de regresión. Por ello su expresión es (Ryan, 1998):
en la cual
siendo npa el número de parámetro de ajuste de la ecuación de RLM. EEM es el error estándar medio también con las unidades de Yi:
Por último, EREM es el error relativo estándar medio, que es adimensional:
Aplicación hidrológica
Antecedentes
Clarke (1994) presenta los datos de 23 estaciones hidrométricas ubicadas dentro del sistema del río Itajaí-Acú en Brasil, relativos al gasto máximo medio anual (Qma, m3/s), áreas de cuenca (A, km2) y número de años de registro (NA). Los valores de Qma variaron de 31 a 3194 m3/s; los de A de 105 a 11 719 km2, y los de NA de 3 a 118 años. Realiza una regresión lineal del tipo Qma = b0 + b1A, cuyo coeficiente de determinación (R2 ) resultó de 0.828. Para mejorar la estimación, se emplea la regresión potencial Qma = b0·-Abl , cuyo R2 resultó de 0.909. El análisis de los residuos de esta última ecuación indica que los mayores errores corresponden a los datos que tienen las menores amplitudes de registro, por ello aplica la técnica de MCP utilizando wi = NA, obteniendo un R2 de 0.935 y observando que los residuos disminuyeron. La función de ponderado aplicada por Clarke (1994) es en realidad una versión bastante simplificada de la función óptima, como se deduce a continuación en su primer procedimiento de búsqueda. Este enfoque simple para la aplicación del ajuste de MCP también ha sido aplicado por Vogel, Wilson y Daly (1999).
Primera aproximación
Procedimiento que divide la varianza residual
Tasker (1980) desarrolló un procedimiento para estimar la función de ponderado óptima wi = 1/σ2 (Yi ), a utilizar en el ajuste por MCP. Partió de los resultados de Matalas y Gilroy (1968), que establecieron que la varianza de la variable dependiente (Y i ) se puede dividir en dos componentes: la primera originada por el error del modelo σ2(δi ) y la segunda debida al error de muestreo σ2(ei); es decir:
El subíndice i varía de uno al número de datos o valores de Yi, y también es igual al número de estaciones hidrométricas (NE) o registros procesados en los análisis regionales. La condición de homocedasticidad requiere que tanto σ2(δ i) como σ2(ei) sean independientes de i. En el procedimiento desarrollado por Tasker (1980), sólo se acepta la independencia del error del modelo. Para estimar la varianza del error de muestreo, se considera que las crecientes asociadas con un cierto periodo de retorno (Tr) siguen una distribución Pearson tipo III y entonces de acuerdo con Bobée (1973) se tiene:
en la cual a es la desviación estándar de los gastos máximos anuales; n, el número de gastos anuales observados en la estación hidrométrica i; γ, el coeficiente de asimetría de los gastos máximos anuales, y Kp es la desviación estandarizada con distribución Pearson tipo III asociada con el valor de γ y de la probabilidad de no excedencia p. Algunas veces Kp se designa por KTr, pues Tr = 1/(1 - p). En estudios hidrológicos regionales se puede aceptar que y y a2 son aproximadamente constantes en todos los sitios de tal zona o región debido precisamente a la homogeneidad regional verificada previamente (Hosking & Wallis, 1997) y entonces la varianza de Yi (ecuación 13) se expresa como (Tasker, 1980):
siendo co = σ2(δi) una constante y c1 otra, que se logra estimar con base en la información regional disponible, según la expresión:
en la cual
en donde
siendo ni el número de datos de cada registro procesado. Tasker (1980) indica que cuando el error debido al modelo es grande,
Cuando se disponga de los datos anuales de cada variable dependiente (Yi), el procedimiento anterior puede ser mejorado, calculando la magnitud de la constante
Aspectos operativos previos
Para aplicar el procedimiento de Tasker (1980), descrito en las ecuaciones (13) a (18), primero se deben definir las expresiones de las estimaciones regionales
siendo:
Para la estimación del valor de
La ecuación anterior está limitada
Análisis de resultados
Los datos que serán procesados se muestran en el Cuadro 2 en sus primeras seis columnas; proceden de Campos-Aranda (2013), y corresponden a los valores del gasto máximo medio anual (Qma) o creciente media anual, así como de varias propiedades fisiográficas de 22 cuencas de las estaciones hidrométricas de la Región Hidrológica 10 (Sinaloa), que no presentan régimen hidrológico modificado, las cuales fueron tomadas de Escalante-Sandoval y Reyes-Chávez (2002). La regresión del tipo Qma = b0•Abl se obtuvo con el método de MCO, b0 = 11.6751 y b1 = 0.5258, y las estimaciones y residuales mostrados en las columnas 2 y 3 del Cuadro 3. Este método condujo a un valor de R2 de 0.813, con el resto de indicadores de ajuste mostrados al final de la citada columna 3 y un error estándar de ajuste de 309.8 m3/s.
Simbología:
1 mínimos cuadrados ordinarios.
2 mínimos cuadrados ponderados, con función de ponderado con procedimiento que divide la varianza residual.
3 mínimos cuadrados ponderados, con función de ponderado con procedimiento basado en datos cercanos.
Para los datos de la columna 3 del Cuadro 2 se obtiene que su coeficiente de asimetría es 2.0303, lo cual conduce a un valor de K = -0.3104. Con base en estos valores y los citados anteriormente, se obtuvo que c0 = 88 468.7 y c1 = 27 7787.1, ambos con unidades de varianza (m6/s2). Los correspondientes factores de ponderación (ecuación (18)) se muestran en la columna 7 del Cuadro 2.
El método de MCP aporta b0 = 11.7339 y b1 = 0.5251, con las estimaciones de Qma, los residuos y sus indicadores de desempeño (ID) que se tienen en las columnas 4 y 5 del Cuadro 3. Se observa en la porción final de la columna 5 que los ID del método de MCP son casi iguales a los del ajuste de MCO. Este resultado se considera congruente, pues en esta aplicación numérica, las amplitudes de los registros procesados son semejantes (ver Cuadro 2), variando de 21 a 56, con una media de 37 años.
También se revisó la regresión potencial del tipo
Segunda aproximación
Procedimiento basado en datos cercanos
Sugerido por Draper y Smith (1998), y por Montgomery et al. (2002), comienza por definir conjuntos de valores del regresor X1 que son "vecinos cercanos", por tener observaciones con magnitudes semejantes de X1. El procedimiento supone que tales conjuntos pueden considerarse "puntos de repetición" y por lo tanto se puede usar la varianza promedio de sus respuestas (Yi ) para estimar la forma en que, de manera aproximada, cambia Var(Y) en función de X1.
En la columna 9 del Cuadro 2 se indica el renglón donde comienza cada conjunto de datos repetidos, el número de elementos que incluye entre paréntesis y el valor promedio de X1, es decir, del área de cuenca (A). En la columna 10 se indica la varianza muestral (promedio, aproximadamente) de las Yi de cada conjunto, estimada con la expresión siguiente:
Con base en las seis parejas de valores calculados de
con a0 = -2 025.336, a1 = 26.55365, a2 = -1.122622E-02, a3 = 1.481098E-06, a4 = -4.118003E-11 y error estándar de la estimación de 9 315 m3/s. Con base en tal ecuación se obtuvieron los factores de ponderación (wi ) mostrados en la columna 11 del Cuadro 2.
Cuando se tienen dos o más regresores, resulta muy difícil la identificación visual de los datos que son vecinos cercanos y por ello se debe aplicar una técnica analítica para buscar pares de puntos cercanos entre sí en el espacio de Xp (Montgomery et al., 2002), o bien aplicar la ecuación (8), con la función de ponderado estimada con el primer regresor (X1), que es el más importante.
Análisis de resultados
Al aplicar la técnica de MCP, con la función de ponderado mostrada en la columna final del Cuadro 2 y la regresión del tipo Qma = b0·Ab1, se obtienen b0 = 7.7896 y b1 = 0.5784, con una ligera mejoría en el ajuste, pues ahora se tiene un R2 de 0.824 y los valores del error disminuyen (ver EEM y EREM) al final de la columna 7 del Cuadro 3. Las estimaciones de Qma y sus residuos se tienen en las columnas 6 y 7 del citado Cuadro 3. Este ajuste reduce notablemente los residuos positivos de las estaciones hidrométricas Huites y Guatenipa II, lo cual se ve reflejado en el valor del EEM. Lo anterior se puede observar al comparar las Figuras 1 y 2.
Al aplicar la técnica de MCP a la segunda regresión potencial, se obtuvieron estos resultados: b0 = 8.0767, b1 = 0.6563 y b2 = -0.1422. Las estimaciones y los residuos de este método se tienen en las dos columnas finales del Cuadro 3. Los ID relacionados con el error muestran una mejoría de ajuste, como se observa al comparar los tres últimos renglones de las columnas 9 y 13. Al igual que la ecuación potencial anterior, los residuos de las estaciones Huites y Guate-nipa II se reducen de modo sustancial y ello se aprecia en los ID del error.
Conclusiones
En los experimentos numéricos realizados por Tasker (1980) variaron: (1) el error del modelo de 0 a 100%; (2) la correlación entre estaciones tomó valores de 0.0, 0.4 y 0.8, y (3) la amplitud de los registros fluctuó de 10 a 50 años, con tres formas aleatorias de variación. Concluye, con base en la simulación numérica, que cuando se aplica la función de ponderado (wi) definida por la ecuación (18), en el método de mínimos cuadrados ponderados (MCP), la ecuación de regresión resultante siempre tiene mejores indicadores de desempeño (ecuaciones (9) a (12)), que la obtenida con mínimos cuadrados ordinarios (MCO); excepto cuando ni no varía y/o existe correlación cruzada importante entre los eventos anuales de las variables dependientes. En este último caso, habrá que aplicar la técnica de mínimos cuadrados generalizados (Griffis & Stedinger, 2007).
Con base en la aplicación hidrológica descrita, se pudo verificar que siempre alguno o varios de los indicadores de desempeño (ID) mostraron un mejor ajuste, es decir, se redujeron, al aplicar la técnica de MCP, en comparación con los ID obtenidos por MCO. Para el caso mostrado, en general los errores residuales se reducen más (se obtienen valores menores de los ID) con la segunda función de ponderado, la cual se obtiene con base en los datos cercanos.
Por lo anterior, se recomienda revisar las ecuaciones potenciales obtenidas para estimar Qma, o bien QTr , a través del ajuste de MCP, empleando al menos la primera función de ponderado descrita, pues la segunda requiere la ocurrencia de datos cercanos. Incluso en su versión simplificada, como la utilizó Clarke (1994) y Vogel et al. (1999), es decir, con w. = NA, podrá aportar una mejoría estadística, esto es, reducir los ID o cuando menos verificar la similitud numérica de resultados del ajuste.