Main Content

Esta página es para la versión anterior. La página correspondiente en inglés ha sido eliminada en la versión actual.

Regresión no lineal ponderada

En este ejemplo se muestra cómo ajustar un modelo de regresión no lineal para datos con varianza de error no constante.

Los algoritmos regulares de mínimos cuadrados no lineales son apropiados cuando los errores de medición tienen la misma varianza. Cuando esa suposición no es cierta, es apropiado utilizar un ajuste ponderado. En este ejemplo se muestra cómo utilizar ponderaciones con la función.fitnlm

Datos y modelo para el ajuste

Usaremos los datos recopilados para estudiar la contaminación del agua causada por los desechos industriales y domésticos. Estos datos se describen en detalle en Box, G.P., W.G. Hunter y J.S. Hunter, Statistics for Experimenters (Wiley, 1978, págs. 483-487). La variable de respuesta es la demanda de oxígeno bioquímico en mg/l, y la variable predictora es el tiempo de incubación en días.

x = [1 2 3 5 7 10]'; y = [109 149 149 191 213 224]';  plot(x,y,'ko'); xlabel('Incubation (days), x'); ylabel('Biochemical oxygen demand (mg/l), y'); 

Suponeremos que se sabe que las dos primeras observaciones se hicieron con menos precisión que las observaciones restantes. Podrían, por ejemplo, haberse fabricado con un instrumento diferente. Otra razón común para los datos de peso es que cada observación registrada es en realidad la media de varias mediciones tomadas al mismo valor de x. En los datos aquí, supongamos que los dos primeros valores representan una sola medición sin procesar, mientras que los cuatro restantes son cada uno la media de 5 mediciones sin procesar. Entonces sería apropiado ponderar por el número de mediciones que entraron en cada observación.

w = [1 1 5 5 5 5]'; 

El modelo que ajustaremos a estos datos es una curva exponencial escalada que se convierte en nivel a medida que x se vuelve grande.

modelFun = @(b,x) b(1).*(1-exp(-b(2).*x)); 

Sólo basado en un ajuste visual áspero, parece que una curva dibujada a través de los puntos podría nivelar se nivela en un valor de alrededor de 240 en algún lugar de la vecindad de x 15. Por lo tanto, usaremos 240 como el valor inicial para b1, y dado que e(-.5*15) es pequeño en comparación con 1, usaremos .5 como el valor inicial para b2.

start = [240; .5]; 

Ajustar el modelo sin pesos

El peligro de ignorar el error de medición es que el ajuste puede estar demasiado influenciado por mediciones imprecisas y, por lo tanto, puede no proporcionar un buen ajuste a las mediciones que se conocen con precisión. Vamos a ajustar los datos sin ponderaciones y compararlos con los puntos.

nlm = fitnlm(x,y,modelFun,start); xx = linspace(0,12)'; line(xx,predict(nlm,xx),'linestyle','--','color','k') 

Ajustar el modelo con pesos

Observe que la curva ajustada se tira hacia los dos primeros puntos, pero parece perder se pierde la tendencia de los otros puntos. Vamos a tratar de repetir el ajuste usando pesas.

wnlm = fitnlm(x,y,modelFun,start,'Weight',w) line(xx,predict(wnlm,xx),'color','b') 
 wnlm =    Nonlinear regression model:     y ~ b1*(1 - exp( - b2*x))  Estimated Coefficients:           Estimate       SE       tStat       pValue             ________    ________    ______    __________      b1     225.17         10.7    21.045    3.0134e-05     b2    0.40078     0.064296    6.2333     0.0033745   Number of observations: 6, Error degrees of freedom: 4 Root Mean Squared Error: 24 R-Squared: 0.908,  Adjusted R-Squared 0.885 F-statistic vs. zero model: 696, p-value = 8.2e-06 

La desviación estándar de población estimada en este caso describe la variación media para una observación "estándar" con un peso, o precisión de medición, de 1.

wnlm.RMSE 
 ans =     24.0096  

Una parte importante de cualquier análisis es una estimación de la precisión del ajuste del modelo. La pantalla del coeficiente muestra errores estándar para los parámetros, pero también podemos calcular los intervalos de confianza para ellos.

coefCI(wnlm) 
 ans =    195.4650  254.8788     0.2223    0.5793  

Estimar la curva de respuesta

A continuación, calcularemos los valores de respuesta ajustados y los intervalos de confianza para ellos. De forma predeterminada, esos anchos son para límites de confianza puntuales para el valor predicho, pero solicitaremos intervalos simultáneos para toda la curva.

[ypred,ypredci] = predict(wnlm,xx,'Simultaneous',true); plot(x,y,'ko', xx,ypred,'b-', xx,ypredci,'b:'); xlabel('x'); ylabel('y'); legend({'Data', 'Weighted fit', '95% Confidence Limits'},'location','SouthEast'); 

Observe que los dos puntos ponderados hacia abajo no se ajustan tan bien por la curva como los puntos restantes. Eso es lo que cabría esperar de un ajuste ponderado.

También es posible estimar los intervalos de predicción para futuras observaciones en valores especificados de x. Esos intervalos asumirán en efecto un peso, o precisión de medición, de 1.

[ypred,ypredci] = predict(wnlm,xx,'Simultaneous',true,'Prediction','observation'); plot(x,y,'ko', xx,ypred,'b-', xx,ypredci,'b:'); xlabel('x'); ylabel('y'); legend({'Data', 'Weighted fit', '95% Prediction Limits'},'location','SouthEast'); 

La escala absoluta de las ponderaciones en realidad no afecta a las estimaciones de parámetros. Reescalar los pesos por cualquier constante nos habría dado las mismas estimaciones. Pero sí afectan los límites de confianza, ya que los límites representan una observación con el peso 1. Aquí se puede ver que los puntos con mayor peso parecen demasiado cerca de la línea ajustada, en comparación con los límites de confianza.

Si bien el método no nos permite cambiar los pesos, es posible que hagamos un post-procesamiento e investiguemos cómo la curva buscaría una estimación más precisa.predict Supongamos que estamos interesados en una nueva observación que se basa en el promedio de cinco mediciones, al igual que los últimos cuatro puntos de esta gráfica. Podríamos reducir el ancho de los intervalos por un factor de sqrt(5).

halfwidth = ypredci(:,2)-ypred; newwidth = halfwidth/sqrt(5); newci = [ypred-newwidth, ypred+newwidth]; plot(x,y,'ko', xx,ypred,'b-', xx,newci,'r:'); xlabel('x'); ylabel('y'); legend({'Data', 'Weighted fit', 'Limits for weight=5'},'location','SouthEast'); 

Análisis residual

Además de trazar los datos y el ajuste, trazaremos los residuos de un ajuste contra los predictores, para diagnosticar cualquier problema con el modelo. Los residuos deben aparecer independientes e distribuidos de forma idéntica, pero con una varianza proporcional a la inversa de los pesos. Podemos estandarizar esta varianza para que la gráfica sea más fácil de interpretar.

r = wnlm.Residuals.Raw; plot(x,r.*sqrt(w),'b^'); xlabel('x'); ylabel('Residuals, yFit - y'); 

Hay algunas pruebas de patrones sistemáticos en esta trama residual. Observe cómo los últimos cuatro residuos tienen una tendencia lineal, lo que sugiere que el modelo podría no aumentar lo suficientemente rápido a medida que x aumenta. Además, la magnitud de los residuos tiende a disminuir a medida que x aumenta, lo que sugiere que el error de medición puede depender de x. Estos merecen una investigación, sin embargo, hay tan pocos puntos de datos, que es difícil dar importancia a estos patrones aparentes.