Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

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 de mínimos cuadrados no lineales normales son apropiados cuando todos los errores de medición tienen la misma varianza. Cuando esa suposición no es cierta, conviene utilizar un ajuste ponderado. En este ejemplo se muestra cómo utilizar ponderaciones con la función.fitnlm

Datos y modelo para el Fit

Utilizaremos los datos recopilados para estudiar la contaminación del agua causada por los desechos industriales y domésticos. Estos datos se describen detalladamente en Box, G.P., W.G. Hunter y J.S. Hunter, Statistics for Experimenters (Wiley, 1978, PP. 483-487). La variable de respuesta es la demanda bioquímica de oxígeno 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'); 

Supondremos que se sabe que las dos primeras observaciones se realizaron con menos precisión que las observaciones restantes. Por ejemplo, podrían haberse hecho con un instrumento diferente. Otra razón común para los datos de peso es que cada observación grabada es en realidad la media de varias mediciones tomadas en el mismo valor de x. En los datos aquí, supongamos que los dos primeros valores representan una sola medida RAW, mientras que los cuatro restantes son cada uno la media de 5 mediciones en bruto. Entonces sería apropiado ponderar por el número de mediciones que iban en cada observación.

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

El modelo que ajustaremos a estos datos es una curva exponencial a escala 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 trazada a través de los puntos podría nivelarse a un valor de alrededor de 240 en algún lugar en el vecindario de x = 15. Así que usaremos 240 como valor inicial para B1, y puesto que e ^ (-.5 * 15) es pequeño comparado con 1, usaremos 5 como el valor inicial para B2.

start = [240; .5]; 

Ajuste el modelo sin pesos

El peligro de ignorar el error de medición es que el ajuste puede estar excesivamente 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') 

Ajuste el modelo con pesos

Observe que la curva ajustada se tira hacia los dos primeros puntos, pero parece perderse la tendencia de los otros puntos. Intentemos repetir el ajuste usando pesos.

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 promedio 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 visualización del coeficiente muestra los errores estándar de 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 los límites de confianza de pointwise para el valor pronosticado, pero solicitamos 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 descendentes no son aptos también por la curva como los puntos restantes. Eso es lo que cabría esperar para un ajuste ponderado.

También es posible estimar los intervalos de predicción para futuras observaciones a valores especificados de x. Esos intervalos, en efecto, asumirán 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 no afecta realmente a las estimaciones de parámetros. La reescalado de los pesos por cualquier constante nos habría dado las mismas estimaciones. Pero sí afectan a 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.

Mientras que el método no nos permite cambiar las ponderaciones, es posible para nosotros hacer un poco de post-procesamiento e investigar cómo la curva se vería para 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 en esta trama. 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, vamos a trazar los residuos de un ajuste en contra de los predictores, para diagnosticar cualquier problema con el modelo. Los residuos deben aparecer independientes e idénticos distribuidos, pero con una varianza proporcional a la inversa de las ponderaciones. Podemos estandarizar esta varianza para que la trama 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 evidencias de patrones sistemáticos en esta trama residual. Observe cómo los últimos cuatro residuales tienen una tendencia lineal, lo que sugiere que el modelo podría no aumentar lo suficientemente rápido como x aumenta. Además, la magnitud de los residuos tiende a disminuir a medida que aumenta x, 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 de adjuntar significado a estos patrones aparentes.